crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #2095 from paulromano/get-nuclide-atoms

Browse source 
Add method on Material for getting number of atoms
This commit is contained in:
Jonathan Shimwell 2022-06-29 23:52:43 +01:00 • committed by GitHub
commit 2f8f8d06d6
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
11 changed files with 143 additions and 67 deletions
Download patch file Download diff file Expand all files Collapse all files

1
docs/source/pythonapi/data.rst
View file

@ -61,6 +61,7 @@ Core Functions
atomic_mass
atomic_weight
decay_constant
dose_coefficients
gnd_name
half_life

8
openmc/cell.py
View file

@ -214,8 +214,8 @@ class Cell(IDManagerMixin):
self._atoms = self._fill.get_nuclide_atom_densities()
# Convert to total number of atoms
for key, nuclide in self._atoms.items():
atom = nuclide[1] * self._volume * 1.0e+24
for key, atom_per_bcm in self._atoms.items():
atom = atom_per_bcm * self._volume * 1.0e+24
self._atoms[key] = atom
elif self.fill_type == 'distribmat':
@ -224,12 +224,12 @@ class Cell(IDManagerMixin):
partial_volume = self.volume / len(self.fill)
self._atoms = OrderedDict()
for mat in self.fill:
for key, nuclide in mat.get_nuclide_atom_densities().items():
for key, atom_per_bcm in mat.get_nuclide_atom_densities().items():
# To account for overlap of nuclides between distribmat
# we need to append new atoms to any existing value
# hence it is necessary to ask for default.
atom = self._atoms.setdefault(key, 0)
atom += nuclide[1] * partial_volume * 1.0e+24
atom += atom_per_bcm * partial_volume * 1.0e+24
self._atoms[key] = atom
else:

35
openmc/data/data.py
View file

@ -3,7 +3,7 @@ import json
import os
import re
from pathlib import Path
from math import sqrt
from math import sqrt, log
from warnings import warn
# Isotopic abundances from Meija J, Coplen T B, et al, "Isotopic compositions
@ -202,7 +202,7 @@ _GND_NAME_RE = re.compile(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)')
# Used in half_life function as a cache
_HALF_LIFE = {}
_LOG_TWO = log(2.0)
def atomic_mass(isotope):
"""Return atomic mass of isotope in atomic mass units.
@ -283,6 +283,8 @@ def half_life(isotope):
Half-life values are from the `ENDF/B-VIII.0 decay sublibrary
<https://www.nndc.bnl.gov/endf-b8.0/download.html>`_.
.. versionadded:: 0.13.1
Parameters
----------
isotope : str
@ -302,6 +304,35 @@ def half_life(isotope):
return _HALF_LIFE.get(isotope.lower())
def decay_constant(isotope):
"""Return decay constant of isotope in [s^-1]
Decay constants are based on half-life values from the
:func:`~openmc.data.half_life` function. When the isotope is stable, a decay
constant of zero is returned.
.. versionadded:: 0.13.1
Parameters
----------
isotope : str
Name of isotope, e.g., 'Pu239'
Returns
-------
float
Decay constant of isotope in [s^-1]
See also
--------
openmc.data.half_life
"""
t = half_life(isotope)
return _LOG_TWO / t if t else 0.0
def water_density(temperature, pressure=0.1013):
"""Return the density of liquid water at a given temperature and pressure.

9
openmc/data/decay.py
View file

@ -465,11 +465,10 @@ class Decay(EqualityMixin):
@property
def decay_constant(self):
if hasattr(self.half_life, 'n'):
return log(2.)/self.half_life
else:
mu, sigma = self.half_life
return ufloat(log(2.)/mu, log(2.)/mu**2*sigma)
if self.half_life.n == 0.0:
name = self.nuclide['name']
raise ValueError(f"{name} is listed as unstable but has a zero half-life.")
return log(2.)/self.half_life
@property
def decay_energy(self):

15
openmc/deplete/operator.py
View file

@ -528,9 +528,9 @@ class Operator(TransportOperator):
"""
mat_id = str(mat.id)
for nuclide, density in mat.get_nuclide_atom_densities().values():
number = density * 1.0e24
self.number.set_atom_density(mat_id, nuclide, number)
for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items():
atom_per_cc = atom_per_bcm * 1.0e24
self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
def _set_number_from_results(self, mat, prev_res):
"""Extracts material nuclides and number densities.
@ -556,16 +556,15 @@ class Operator(TransportOperator):
# Merge lists of nuclides, with the same order for every calculation
geom_nuc_densities.update(depl_nuc)
for nuclide in geom_nuc_densities.keys():
for nuclide, atom_per_bcm in geom_nuc_densities.items():
if nuclide in depl_nuc:
concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1]
volume = prev_res[-1].volume[mat_id]
number = concentration / volume
atom_per_cc = concentration / volume
else:
density = geom_nuc_densities[nuclide][1]
number = density * 1.0e24
atom_per_cc = atom_per_bcm * 1.0e24
self.number.set_atom_density(mat_id, nuclide, number)
self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
def initial_condition(self):
"""Performs final setup and returns initial condition.

72
openmc/material.py
View file

@ -3,7 +3,6 @@ from collections.abc import Iterable
from copy import deepcopy
from numbers import Real
from pathlib import Path
import math
import re
import warnings
from xml.etree import ElementTree as ET
@ -86,6 +85,9 @@ class Material(IDManagerMixin):
fissionable_mass : float
Mass of fissionable nuclides in the material in [g]. Requires that the
:attr:`volume` attribute is set.
activity : float
Activity of the material in [Bq]. Requires that the :attr:`volume`
attribute is set.
"""
@ -142,20 +144,10 @@ class Material(IDManagerMixin):
return string
@property
def activity(self):
"""Returns the total activity of the material in Becquerels."""
atoms_per_barn_cm = self.get_nuclide_atom_densities()
total_activity = 0
for key, value in atoms_per_barn_cm.items():
half_life = openmc.data.half_life(key)
if half_life:
total_activity += value[1] / half_life
total_activity *= math.log(2) * 1e24 * self.volume
return total_activity
return sum(self.get_nuclide_activity().values())
@property
def name(self):
@ -259,10 +251,10 @@ class Material(IDManagerMixin):
if self.volume is None:
raise ValueError("Volume must be set in order to determine mass.")
density = 0.0
for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values():
for nuc, atoms_per_bcm in self.get_nuclide_atom_densities().items():
Z = openmc.data.zam(nuc)[0]
if Z >= 90:
density += 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
density += 1e24 * atoms_per_bcm * openmc.data.atomic_mass(nuc) \
/ openmc.data.AVOGADRO
return density*self.volume
@ -786,11 +778,15 @@ class Material(IDManagerMixin):
"""Returns all nuclides in the material and their atomic densities in
units of atom/b-cm
.. versionchanged:: 0.13.1
The values in the dictionary were changed from a tuple containing
the nuclide name and the density to just the density.
Returns
-------
nuclides : dict
Dictionary whose keys are nuclide names and values are tuples of
(nuclide, density in atom/b-cm)
Dictionary whose keys are nuclide names and values are densities in
[atom/b-cm]
"""
@ -850,10 +846,46 @@ class Material(IDManagerMixin):
nuclides = OrderedDict()
for n, nuc in enumerate(nucs):
nuclides[nuc] = (nuc, nuc_densities[n])
nuclides[nuc] = nuc_densities[n]
return nuclides
def get_nuclide_activity(self):
"""Return activity in [Bq] for each nuclide in the material
.. versionadded:: 0.13.1
Returns
-------
dict
Dictionary whose keys are nuclide names and values are activity in
[Bq].
"""
activity = {}
for nuclide, atoms in self.get_nuclide_atoms().items():
inv_seconds = openmc.data.decay_constant(nuclide)
activity[nuclide] = inv_seconds * atoms
return activity
def get_nuclide_atoms(self):
"""Return number of atoms of each nuclide in the material
.. versionadded:: 0.13.1
Returns
-------
dict
Dictionary whose keys are nuclide names and values are number of
atoms present in the material.
"""
if self.volume is None:
raise ValueError("Volume must be set in order to determine atoms.")
atoms = {}
for nuclide, atom_per_bcm in self.get_nuclide_atom_densities().items():
atoms[nuclide] = 1.0e24 * atom_per_bcm * self.volume
return atoms
def get_mass_density(self, nuclide=None):
"""Return mass density of one or all nuclides
@ -870,9 +902,9 @@ class Material(IDManagerMixin):
"""
mass_density = 0.0
for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values():
for nuc, atoms_per_bcm in self.get_nuclide_atom_densities().items():
if nuclide is None or nuclide == nuc:
density_i = 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
density_i = 1e24 * atoms_per_bcm * openmc.data.atomic_mass(nuc) \
/ openmc.data.AVOGADRO
mass_density += density_i
return mass_density
@ -1092,7 +1124,7 @@ class Material(IDManagerMixin):
nuclides_per_cc = defaultdict(float)
mass_per_cc = defaultdict(float)
for mat, wgt in zip(materials, wgts):
for nuc, atoms_per_bcm in mat.get_nuclide_atom_densities().values():
for nuc, atoms_per_bcm in mat.get_nuclide_atom_densities().items():
nuc_per_cc = wgt*1.e24*atoms_per_bcm
nuclides_per_cc[nuc] += nuc_per_cc
mass_per_cc[nuc] += nuc_per_cc*openmc.data.atomic_mass(nuc) / \

12
openmc/plotter.py
View file

@ -542,15 +542,11 @@ def _calculate_cexs_elem_mat(this, types, temperature=294.,
if isinstance(this, openmc.Material):
# Expand elements in to nuclides with atomic densities
nuclides = this.get_nuclide_atom_densities()
# For ease of processing split out the nuclide and its fraction
nuc_fractions = {nuclide[1][0]: nuclide[1][1]
for nuclide in nuclides.items()}
nuc_fractions = this.get_nuclide_atom_densities()
# Create a dict of [nuclide name] = nuclide object to carry forward
# with a common nuclides format between openmc.Material and
# openmc.Element objects
nuclides = {nuclide[1][0]: nuclide[1][0]
for nuclide in nuclides.items()}
nuclides = {nuclide: nuclide for nuclide in nuc_fractions}
else:
# Expand elements in to nuclides with atomic densities
nuclides = this.expand(1., 'ao', enrichment=enrichment,
@ -885,13 +881,13 @@ def _calculate_mgxs_elem_mat(this, types, library, orders=None,
# Check to see if we have nuclides/elements or a macroscopic object
if this._macroscopic is not None:
# We have macroscopics
nuclides = {this._macroscopic: (this._macroscopic, this.density)}
nuclides = {this._macroscopic: this.density}
else:
# Expand elements in to nuclides with atomic densities
nuclides = this.get_nuclide_atom_densities()
# For ease of processing split out nuc and nuc_density
nuc_fraction = [nuclide[1][1] for nuclide in nuclides.items()]
nuc_fraction = list(nuclides.values())
else:
T = temperature
# Expand elements in to nuclides with atomic densities

13
tests/unit_tests/test_data_misc.py
View file

@ -1,6 +1,6 @@
#!/usr/bin/env python
from collections.abc import Mapping
from math import log
import os
from pathlib import Path
@ -126,3 +126,14 @@ def test_zam():
assert openmc.data.zam('Am242_m10') == (95, 242, 10)
with pytest.raises(ValueError):
openmc.data.zam('garbage')
def test_half_life():
assert openmc.data.half_life('H2') is None
assert openmc.data.half_life('U235') == pytest.approx(2.22102e16)
assert openmc.data.half_life('Am242') == pytest.approx(57672.0)
assert openmc.data.half_life('Am242_m1') == pytest.approx(4449622000.0)
assert openmc.data.decay_constant('H2') == 0.0
assert openmc.data.decay_constant('U235') == pytest.approx(log(2.0)/2.22102e16)
assert openmc.data.decay_constant('Am242') == pytest.approx(log(2.0)/57672.0)
assert openmc.data.decay_constant('Am242_m1') == pytest.approx(log(2.0)/4449622000.0)

2
tests/unit_tests/test_deplete_activation.py
View file

@ -92,7 +92,7 @@ def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts
w = model.geometry.get_materials_by_name('tungsten')[0]
atom_densities = w.get_nuclide_atom_densities()
atom_per_cc = 1e24 * atom_densities['W186'][1] # Density in atom/cm^3
atom_per_cc = 1e24 * atom_densities['W186'] # Density in atom/cm^3
n0 = atom_per_cc * w.volume # Absolute number of atoms
# Pick a random irradiation time and then determine necessary source rate to

25
tests/unit_tests/test_material.py
View file

@ -100,7 +100,7 @@ def test_add_elements_by_formula():
m.add_elements_from_formula('Li4SiO4')
# checking the ratio of elements is 4:1:4 for Li:Si:O
elem = defaultdict(float)
for nuclide, adens in m.get_nuclide_atom_densities().values():
for nuclide, adens in m.get_nuclide_atom_densities().items():
if nuclide.startswith("Li"):
elem["Li"] += adens
if nuclide.startswith("Si"):
@ -117,7 +117,7 @@ def test_add_elements_by_formula():
'O16': 0.443386, 'O17': 0.000168}
nuc_dens = m.get_nuclide_atom_densities()
for nuclide in ref_dens:
assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2)
assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2)
# testing the correct nuclides are added to the Material when enriched
m = openmc.Material()
@ -129,7 +129,7 @@ def test_add_elements_by_formula():
'O16': 0.443386, 'O17': 0.000168}
nuc_dens = m.get_nuclide_atom_densities()
for nuclide in ref_dens:
assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2)
assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2)
# testing the use of brackets
m = openmc.Material()
@ -137,7 +137,7 @@ def test_add_elements_by_formula():
# checking the ratio of elements is 2:2:6 for Mg:N:O
elem = defaultdict(float)
for nuclide, adens in m.get_nuclide_atom_densities().values():
for nuclide, adens in m.get_nuclide_atom_densities().items():
if nuclide.startswith("Mg"):
elem["Mg"] += adens
if nuclide.startswith("N"):
@ -155,7 +155,7 @@ def test_add_elements_by_formula():
'O16': 0.599772, 'O17': 0.000227}
nuc_dens = m.get_nuclide_atom_densities()
for nuclide in ref_dens:
assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2)
assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2)
# testing non integer multiplier results in a value error
m = openmc.Material()
@ -289,12 +289,21 @@ def test_get_nuclide_densities(uo2):
def test_get_nuclide_atom_densities(uo2):
nucs = uo2.get_nuclide_atom_densities()
for nuc, density in nucs.values():
for nuc, density in uo2.get_nuclide_atom_densities().items():
assert nuc in ('U235', 'O16')
assert density > 0
def test_get_nuclide_atoms():
mat = openmc.Material()
mat.add_nuclide('Li6', 1.0)
mat.set_density('atom/cm3', 3.26e20)
mat.volume = 100.0
atoms = mat.get_nuclide_atoms()
assert atoms['Li6'] == pytest.approx(mat.density * mat.volume)
def test_mass():
m = openmc.Material()
m.add_nuclide('Zr90', 1.0, 'wo')
@ -341,7 +350,7 @@ def test_borated_water():
'O16':2.4672e-02}
nuc_dens = m.get_nuclide_atom_densities()
for nuclide in ref_dens:
assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2)
assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2)
assert m.id == 50
# Test the Celsius conversion.

18
tests/unit_tests/test_model.py
View file

@ -391,16 +391,14 @@ def test_py_lib_attributes(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
assert openmc.lib.materials[1].get_density('atom/b-cm') == \
pytest.approx(0.06891296988603757, abs=1e-13)
mat_a_dens = np.sum(
[v[1] for v in test_model.materials[0].
get_nuclide_atom_densities().values()])
list(test_model.materials[0].get_nuclide_atom_densities().values()))
assert mat_a_dens == pytest.approx(0.06891296988603757, abs=1e-8)
# Change the density
test_model.update_densities(['UO2'], 2.)
assert openmc.lib.materials[1].get_density('atom/b-cm') == \
pytest.approx(2., abs=1e-13)
mat_a_dens = np.sum(
[v[1] for v in test_model.materials[0].
get_nuclide_atom_densities().values()])
list(test_model.materials[0].get_nuclide_atom_densities().values()))
assert mat_a_dens == pytest.approx(2., abs=1e-8)
# Now lets do the cell temperature updates.
@ -441,7 +439,7 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
test_model = openmc.Model(geom, mats, settings, tals, plots)
initial_mat = mats[0].clone()
initial_u = initial_mat.get_nuclide_atom_densities()['U235'][1]
initial_u = initial_mat.get_nuclide_atom_densities()['U235']
# Note that the chain file includes only U-235 fission to a stable Xe136 w/
# a yield of 100%. Thus all the U235 we lose becomes Xe136
@ -453,8 +451,8 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
operator_kwargs=op_kwargs,
power=1., output=False)
# Get the new Xe136 and U235 atom densities
after_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1]
after_u = mats[0].get_nuclide_atom_densities()['U235'][1]
after_xe = mats[0].get_nuclide_atom_densities()['Xe136']
after_u = mats[0].get_nuclide_atom_densities()['U235']
assert after_xe + after_u == pytest.approx(initial_u, abs=1e-15)
assert test_model.is_initialized is False
@ -462,7 +460,7 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
mats[0].nuclides.clear()
densities = initial_mat.get_nuclide_atom_densities()
tot_density = 0.
for nuc, density in densities.values():
for nuc, density in densities.items():
mats[0].add_nuclide(nuc, density)
tot_density += density
mats[0].set_density('atom/b-cm', tot_density)
@ -473,8 +471,8 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
operator_kwargs=op_kwargs,
power=1., output=False)
# Get the new Xe136 and U235 atom densities
after_lib_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1]
after_lib_u = mats[0].get_nuclide_atom_densities()['U235'][1]
after_lib_xe = mats[0].get_nuclide_atom_densities()['Xe136']
after_lib_u = mats[0].get_nuclide_atom_densities()['U235']
assert after_lib_xe + after_lib_u == pytest.approx(initial_u, abs=1e-15)
assert test_model.is_initialized is True