mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Merge branch 'mcpl2openmc' of github.com:ebknudsen/openmc into mcpl2openmc
This commit is contained in:
erkn
commit
2530d167db
22 changed files with 4471 additions and 96 deletions
|
|
@ -15,7 +15,7 @@
|
|||
|
||||
# sudo docker run image_name:tag_name or ID with no tag sudo docker run ID number
|
||||
|
||||
FROM debian:bullseye-slim
|
||||
FROM debian:bullseye-slim AS dependencies
|
||||
|
||||
# By default this Dockerfile builds OpenMC without DAGMC and LIBMESH support
|
||||
ARG build_dagmc=off
|
||||
|
|
@ -172,6 +172,8 @@ RUN if [ "$build_libmesh" = "on" ]; then \
|
|||
&& rm -rf ${LIBMESH_INSTALL_DIR}/build ${LIBMESH_INSTALL_DIR}/libmesh ; \
|
||||
fi
|
||||
|
||||
FROM dependencies AS build
|
||||
|
||||
# clone and install openmc
|
||||
RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \
|
||||
&& git clone --shallow-submodules --recurse-submodules --single-branch -b ${openmc_branch} --depth=1 ${OPENMC_REPO} \
|
||||
|
|
@ -207,5 +209,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \
|
|||
&& cd ../openmc && pip install .[test,depletion-mpi] \
|
||||
&& python -c "import openmc"
|
||||
|
||||
FROM build AS release
|
||||
|
||||
# Download cross sections (NNDC and WMP) and ENDF data needed by test suite
|
||||
RUN ${HOME}/OpenMC/openmc/tools/ci/download-xs.sh
|
||||
|
|
|
|||
|
|
@ -61,8 +61,10 @@ Core Functions
|
|||
|
||||
atomic_mass
|
||||
atomic_weight
|
||||
decay_constant
|
||||
dose_coefficients
|
||||
gnd_name
|
||||
half_life
|
||||
isotopes
|
||||
linearize
|
||||
thin
|
||||
|
|
|
|||
|
|
@ -36,7 +36,6 @@ public:
|
|||
int threshold;
|
||||
double n_electrons;
|
||||
double binding_energy;
|
||||
xt::xtensor<double, 1> cross_section;
|
||||
vector<Transition> transitions;
|
||||
};
|
||||
|
||||
|
|
@ -80,8 +79,11 @@ public:
|
|||
Tabulated1D coherent_anomalous_real_;
|
||||
Tabulated1D coherent_anomalous_imag_;
|
||||
|
||||
// Photoionization and atomic relaxation data
|
||||
// Photoionization and atomic relaxation data. Subshell cross sections are
|
||||
// stored separately to improve memory access pattern when calculating the
|
||||
// total cross section
|
||||
vector<ElectronSubshell> shells_;
|
||||
xt::xtensor<double, 2> cross_sections_;
|
||||
|
||||
// Compton profile data
|
||||
xt::xtensor<double, 2> profile_pdf_;
|
||||
|
|
|
|||
|
|
@ -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:
|
||||
|
|
|
|||
|
|
@ -1,10 +1,11 @@
|
|||
import itertools
|
||||
from math import sqrt
|
||||
import json
|
||||
import os
|
||||
import re
|
||||
from pathlib import Path
|
||||
from math import sqrt, log
|
||||
from warnings import warn
|
||||
|
||||
|
||||
# Isotopic abundances from Meija J, Coplen T B, et al, "Isotopic compositions
|
||||
# of the elements 2013 (IUPAC Technical Report)", Pure. Appl. Chem. 88 (3),
|
||||
# pp. 293-306 (2013). The "representative isotopic abundance" values from
|
||||
|
|
@ -199,6 +200,9 @@ _ATOMIC_MASS = {}
|
|||
# Regex for GND nuclide names (used in zam function)
|
||||
_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.
|
||||
|
|
@ -273,6 +277,62 @@ def atomic_weight(element):
|
|||
.format(element))
|
||||
|
||||
|
||||
def half_life(isotope):
|
||||
"""Return half-life of isotope in seconds or None if isotope is stable
|
||||
|
||||
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
|
||||
Name of isotope, e.g., 'Pu239'
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Half-life of isotope in [s]
|
||||
|
||||
"""
|
||||
global _HALF_LIFE
|
||||
if not _HALF_LIFE:
|
||||
# Load ENDF/B-VIII.0 data from JSON file
|
||||
half_life_path = Path(__file__).with_name('half_life.json')
|
||||
_HALF_LIFE = json.loads(half_life_path.read_text())
|
||||
|
||||
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.
|
||||
|
||||
|
|
|
|||
|
|
@ -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):
|
||||
|
|
|
|||
3563
openmc/data/half_life.json
Normal file
3563
openmc/data/half_life.json
Normal file
File diff suppressed because it is too large
Load diff
|
|
@ -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.
|
||||
|
|
|
|||
|
|
@ -85,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.
|
||||
|
||||
"""
|
||||
|
||||
|
|
@ -141,6 +144,11 @@ class Material(IDManagerMixin):
|
|||
|
||||
return string
|
||||
|
||||
@property
|
||||
def activity(self):
|
||||
"""Returns the total activity of the material in Becquerels."""
|
||||
return sum(self.get_nuclide_activity().values())
|
||||
|
||||
@property
|
||||
def name(self):
|
||||
return self._name
|
||||
|
|
@ -243,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
|
||||
|
||||
|
|
@ -407,6 +415,23 @@ class Material(IDManagerMixin):
|
|||
if nuclide == nuc.name:
|
||||
self.nuclides.remove(nuc)
|
||||
|
||||
def remove_element(self, element):
|
||||
"""Remove an element from the material
|
||||
|
||||
Parameters
|
||||
----------
|
||||
element : str
|
||||
Element to remove
|
||||
|
||||
"""
|
||||
cv.check_type('element', element, str)
|
||||
|
||||
# If the Material contains the element, delete it
|
||||
for nuc in reversed(self.nuclides):
|
||||
element_name = re.split(r'\d+', nuc.name)[0]
|
||||
if element_name == element:
|
||||
self.nuclides.remove(nuc)
|
||||
|
||||
def add_macroscopic(self, macroscopic):
|
||||
"""Add a macroscopic to the material. This will also set the
|
||||
density of the material to 1.0, unless it has been otherwise set,
|
||||
|
|
@ -753,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]
|
||||
|
||||
"""
|
||||
|
||||
|
|
@ -817,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
|
||||
|
||||
|
|
@ -837,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
|
||||
|
|
@ -1059,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) / \
|
||||
|
|
|
|||
242
openmc/mesh.py
242
openmc/mesh.py
|
|
@ -193,6 +193,81 @@ class StructuredMesh(MeshBase):
|
|||
s1 = (slice(1, None),)*ndim + (slice(None),)
|
||||
return (vertices[s0] + vertices[s1]) / 2
|
||||
|
||||
@property
|
||||
def num_mesh_cells(self):
|
||||
return np.prod(self.dimension)
|
||||
|
||||
def write_data_to_vtk(self, points, filename, datasets, volume_normalization=True):
|
||||
"""Creates a VTK object of the mesh
|
||||
|
||||
Parameters
|
||||
----------
|
||||
points : list or np.array
|
||||
List of (X,Y,Z) tuples.
|
||||
filename : str
|
||||
Name of the VTK file to write.
|
||||
datasets : dict
|
||||
Dictionary whose keys are the data labels
|
||||
and values are the data sets.
|
||||
volume_normalization : bool, optional
|
||||
Whether or not to normalize the data by
|
||||
the volume of the mesh elements.
|
||||
|
||||
Raises
|
||||
------
|
||||
RuntimeError
|
||||
When the size of a dataset doesn't match the number of cells
|
||||
|
||||
Returns
|
||||
-------
|
||||
vtk.vtkStructuredGrid
|
||||
the VTK object
|
||||
"""
|
||||
|
||||
import vtk
|
||||
from vtk.util import numpy_support as nps
|
||||
|
||||
# check that the data sets are appropriately sized
|
||||
errmsg = "The size of the dataset {} should be equal to the number of cells"
|
||||
for label, dataset in datasets.items():
|
||||
if isinstance(dataset, np.ndarray):
|
||||
if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]:
|
||||
raise RuntimeError(errmsg.format(label))
|
||||
else:
|
||||
if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]:
|
||||
raise RuntimeError(errmsg.format(label))
|
||||
cv.check_type('label', label, str)
|
||||
|
||||
vtk_grid = vtk.vtkStructuredGrid()
|
||||
|
||||
vtk_grid.SetDimensions(*self.dimension)
|
||||
|
||||
vtkPts = vtk.vtkPoints()
|
||||
vtkPts.SetData(nps.numpy_to_vtk(points, deep=True))
|
||||
vtk_grid.SetPoints(vtkPts)
|
||||
|
||||
# create VTK arrays for each of
|
||||
# the data sets
|
||||
for label, dataset in datasets.items():
|
||||
dataset = np.asarray(dataset).flatten()
|
||||
|
||||
if volume_normalization:
|
||||
dataset /= self.volumes.flatten()
|
||||
|
||||
dataset_array = vtk.vtkDoubleArray()
|
||||
dataset_array.SetName(label)
|
||||
dataset_array.SetArray(nps.numpy_to_vtk(dataset),
|
||||
dataset.size,
|
||||
True)
|
||||
vtk_grid.GetCellData().AddArray(dataset_array)
|
||||
|
||||
# write the .vtk file
|
||||
writer = vtk.vtkStructuredGridWriter()
|
||||
writer.SetFileName(str(filename))
|
||||
writer.SetInputData(vtk_grid)
|
||||
writer.Write()
|
||||
|
||||
return vtk_grid
|
||||
|
||||
class RegularMesh(StructuredMesh):
|
||||
"""A regular Cartesian mesh in one, two, or three dimensions
|
||||
|
|
@ -273,10 +348,6 @@ class RegularMesh(StructuredMesh):
|
|||
dims = self._dimension
|
||||
return [(u - l) / d for u, l, d in zip(us, ls, dims)]
|
||||
|
||||
@property
|
||||
def num_mesh_cells(self):
|
||||
return np.prod(self._dimension)
|
||||
|
||||
@property
|
||||
def volumes(self):
|
||||
"""Return Volumes for every mesh cell
|
||||
|
|
@ -621,7 +692,42 @@ class RegularMesh(StructuredMesh):
|
|||
root_cell.fill = lattice
|
||||
|
||||
return root_cell, cells
|
||||
|
||||
def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
|
||||
"""Creates a VTK object of the mesh
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : str or pathlib.Path
|
||||
Name of the VTK file to write.
|
||||
datasets : dict
|
||||
Dictionary whose keys are the data labels
|
||||
and values are the data sets.
|
||||
volume_normalization : bool, optional
|
||||
Whether or not to normalize the data by
|
||||
the volume of the mesh elements.
|
||||
Defaults to True.
|
||||
|
||||
Returns
|
||||
-------
|
||||
vtk.vtkStructuredGrid
|
||||
the VTK object
|
||||
"""
|
||||
|
||||
ll, ur = self.lower_left, self.upper_right
|
||||
x_vals = np.linspace(ll[0], ur[0], num=self.dimension[0] + 1)
|
||||
y_vals = np.linspace(ll[1], ur[1], num=self.dimension[1] + 1)
|
||||
z_vals = np.linspace(ll[2], ur[2], num=self.dimension[2] + 1)
|
||||
|
||||
# create points
|
||||
pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals])
|
||||
|
||||
return super().write_data_to_vtk(
|
||||
points=pts_cartesian,
|
||||
filename=filename,
|
||||
datasets=datasets,
|
||||
volume_normalization=volume_normalization
|
||||
)
|
||||
|
||||
def Mesh(*args, **kwargs):
|
||||
warnings.warn("Mesh has been renamed RegularMesh. Future versions of "
|
||||
|
|
@ -821,6 +927,36 @@ class RectilinearMesh(StructuredMesh):
|
|||
|
||||
return element
|
||||
|
||||
def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
|
||||
"""Creates a VTK object of the mesh
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : str or pathlib.Path
|
||||
Name of the VTK file to write.
|
||||
datasets : dict
|
||||
Dictionary whose keys are the data labels
|
||||
and values are the data sets.
|
||||
volume_normalization : bool, optional
|
||||
Whether or not to normalize the data by
|
||||
the volume of the mesh elements.
|
||||
Defaults to True.
|
||||
|
||||
Returns
|
||||
-------
|
||||
vtk.vtkStructuredGrid
|
||||
the VTK object
|
||||
"""
|
||||
# create points
|
||||
pts_cartesian = np.array([[x, y, z] for z in self.z_grid for y in self.y_grid for x in self.x_grid])
|
||||
|
||||
return super().write_data_to_vtk(
|
||||
points=pts_cartesian,
|
||||
filename=filename,
|
||||
datasets=datasets,
|
||||
volume_normalization=volume_normalization
|
||||
)
|
||||
|
||||
|
||||
class CylindricalMesh(StructuredMesh):
|
||||
"""A 3D cylindrical mesh
|
||||
|
|
@ -1012,6 +1148,46 @@ class CylindricalMesh(StructuredMesh):
|
|||
|
||||
return np.multiply.outer(np.outer(V_r, V_p), V_z)
|
||||
|
||||
def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
|
||||
"""Creates a VTK object of the mesh
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : str or pathlib.Path
|
||||
Name of the VTK file to write.
|
||||
datasets : dict
|
||||
Dictionary whose keys are the data labels
|
||||
and values are the data sets.
|
||||
volume_normalization : bool, optional
|
||||
Whether or not to normalize the data by
|
||||
the volume of the mesh elements.
|
||||
Defaults to True.
|
||||
|
||||
Returns
|
||||
-------
|
||||
vtk.vtkStructuredGrid
|
||||
the VTK object
|
||||
"""
|
||||
# create points
|
||||
pts_cylindrical = np.array(
|
||||
[
|
||||
[r, phi, z]
|
||||
for z in self.z_grid
|
||||
for phi in self.phi_grid
|
||||
for r in self.r_grid
|
||||
]
|
||||
)
|
||||
pts_cartesian = np.copy(pts_cylindrical)
|
||||
r, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1]
|
||||
pts_cartesian[:, 0] = r * np.cos(phi)
|
||||
pts_cartesian[:, 1] = r * np.sin(phi)
|
||||
|
||||
return super().write_data_to_vtk(
|
||||
points=pts_cartesian,
|
||||
filename=filename,
|
||||
datasets=datasets,
|
||||
volume_normalization=volume_normalization
|
||||
)
|
||||
|
||||
class SphericalMesh(StructuredMesh):
|
||||
"""A 3D spherical mesh
|
||||
|
|
@ -1204,6 +1380,49 @@ class SphericalMesh(StructuredMesh):
|
|||
|
||||
return np.multiply.outer(np.outer(V_r, V_t), V_p)
|
||||
|
||||
def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
|
||||
"""Creates a VTK object of the mesh
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : str or pathlib.Path
|
||||
Name of the VTK file to write.
|
||||
datasets : dict
|
||||
Dictionary whose keys are the data labels
|
||||
and values are the data sets.
|
||||
volume_normalization : bool, optional
|
||||
Whether or not to normalize the data by
|
||||
the volume of the mesh elements.
|
||||
Defaults to True.
|
||||
|
||||
Returns
|
||||
-------
|
||||
vtk.vtkStructuredGrid
|
||||
the VTK object
|
||||
"""
|
||||
|
||||
# create points
|
||||
pts_spherical = np.array(
|
||||
[
|
||||
[r, theta, phi]
|
||||
for phi in self.phi_grid
|
||||
for theta in self.theta_grid
|
||||
for r in self.r_grid
|
||||
]
|
||||
)
|
||||
pts_cartesian = np.copy(pts_spherical)
|
||||
r, theta, phi = pts_spherical[:, 0], pts_spherical[:, 1], pts_spherical[:, 2]
|
||||
pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta)
|
||||
pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta)
|
||||
pts_cartesian[:, 2] = r * np.cos(phi)
|
||||
|
||||
return super().write_data_to_vtk(
|
||||
points=pts_cartesian,
|
||||
filename=filename,
|
||||
datasets=datasets,
|
||||
volume_normalization=volume_normalization
|
||||
)
|
||||
|
||||
|
||||
class UnstructuredMesh(MeshBase):
|
||||
"""A 3D unstructured mesh
|
||||
|
|
@ -1328,6 +1547,7 @@ class UnstructuredMesh(MeshBase):
|
|||
"been loaded from a statepoint file.")
|
||||
return len(self._centroids)
|
||||
|
||||
|
||||
@centroids.setter
|
||||
def centroids(self, centroids):
|
||||
cv.check_type("Unstructured mesh centroids", centroids,
|
||||
|
|
@ -1373,7 +1593,7 @@ class UnstructuredMesh(MeshBase):
|
|||
|
||||
Parameters
|
||||
----------
|
||||
filename : str
|
||||
filename : str or pathlib.Path
|
||||
Name of the VTK file to write.
|
||||
datasets : dict
|
||||
Dictionary whose keys are the data labels
|
||||
|
|
@ -1381,6 +1601,11 @@ class UnstructuredMesh(MeshBase):
|
|||
volume_normalization : bool
|
||||
Whether or not to normalize the data by the
|
||||
volume of the mesh elements
|
||||
|
||||
Raises
|
||||
------
|
||||
RuntimeError
|
||||
when the size of a dataset doesn't match the number of cells
|
||||
"""
|
||||
|
||||
import vtk
|
||||
|
|
@ -1397,11 +1622,14 @@ class UnstructuredMesh(MeshBase):
|
|||
" mesh information from a statepoint file.")
|
||||
|
||||
# check that the data sets are appropriately sized
|
||||
errmsg = "The size of the dataset {} should be equal to the number of cells"
|
||||
for label, dataset in datasets.items():
|
||||
if isinstance(dataset, np.ndarray):
|
||||
assert dataset.size == self.n_elements
|
||||
if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]:
|
||||
raise RuntimeError(errmsg.format(label))
|
||||
else:
|
||||
assert len(dataset) == self.n_elements
|
||||
if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]:
|
||||
raise RuntimeError(errmsg.format(label))
|
||||
cv.check_type('label', label, str)
|
||||
|
||||
# create data arrays for the cells/points
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -60,6 +60,24 @@ class Univariate(EqualityMixin, ABC):
|
|||
elif distribution == 'mixture':
|
||||
return Mixture.from_xml_element(elem)
|
||||
|
||||
@abstractmethod
|
||||
def sample(n_samples=1, seed=None):
|
||||
"""Sample the univariate distribution
|
||||
|
||||
Parameters
|
||||
----------
|
||||
n_samples : int
|
||||
Number of sampled values to generate
|
||||
seed : int or None
|
||||
Initial random number seed.
|
||||
|
||||
Returns
|
||||
-------
|
||||
numpy.ndarray
|
||||
A 1-D array of sampled values
|
||||
"""
|
||||
pass
|
||||
|
||||
|
||||
class Discrete(Univariate):
|
||||
"""Distribution characterized by a probability mass function.
|
||||
|
|
@ -115,6 +133,18 @@ class Discrete(Univariate):
|
|||
cv.check_greater_than('discrete probability', pk, 0.0, True)
|
||||
self._p = p
|
||||
|
||||
def cdf(self):
|
||||
return np.insert(np.cumsum(self.p), 0, 0.0)
|
||||
|
||||
def sample(self, n_samples=1, seed=None):
|
||||
np.random.seed(seed)
|
||||
return np.random.choice(self.x, n_samples, p=self.p)
|
||||
|
||||
def normalize(self):
|
||||
"""Normalize the probabilities stored on the distribution"""
|
||||
norm = sum(self.p)
|
||||
self.p = [val / norm for val in self.p]
|
||||
|
||||
def to_xml_element(self, element_name):
|
||||
"""Return XML representation of the discrete distribution
|
||||
|
||||
|
|
@ -240,6 +270,10 @@ class Uniform(Univariate):
|
|||
t.c = [0., 1.]
|
||||
return t
|
||||
|
||||
def sample(self, n_samples=1, seed=None):
|
||||
np.random.seed(seed)
|
||||
return np.random.uniform(self.a, self.b, n_samples)
|
||||
|
||||
def to_xml_element(self, element_name):
|
||||
"""Return XML representation of the uniform distribution
|
||||
|
||||
|
|
@ -341,6 +375,14 @@ class PowerLaw(Univariate):
|
|||
cv.check_type('power law exponent', n, Real)
|
||||
self._n = n
|
||||
|
||||
def sample(self, n_samples=1, seed=None):
|
||||
np.random.seed(seed)
|
||||
xi = np.random.rand(n_samples)
|
||||
pwr = self.n + 1
|
||||
offset = self.a**pwr
|
||||
span = self.b**pwr - offset
|
||||
return np.power(offset + xi * span, 1/pwr)
|
||||
|
||||
def to_xml_element(self, element_name):
|
||||
"""Return XML representation of the power law distribution
|
||||
|
||||
|
|
@ -414,6 +456,16 @@ class Maxwell(Univariate):
|
|||
cv.check_greater_than('Maxwell temperature', theta, 0.0)
|
||||
self._theta = theta
|
||||
|
||||
def sample(self, n_samples=1, seed=None):
|
||||
np.random.seed(seed)
|
||||
return self.sample_maxwell(self.theta, n_samples)
|
||||
|
||||
@staticmethod
|
||||
def sample_maxwell(t, n_samples):
|
||||
r1, r2, r3 = np.random.rand(3, n_samples)
|
||||
c = np.cos(0.5 * np.pi * r3)
|
||||
return -t * (np.log(r1) + np.log(r2) * c * c)
|
||||
|
||||
def to_xml_element(self, element_name):
|
||||
"""Return XML representation of the Maxwellian distribution
|
||||
|
||||
|
|
@ -502,6 +554,13 @@ class Watt(Univariate):
|
|||
cv.check_greater_than('Watt b', b, 0.0)
|
||||
self._b = b
|
||||
|
||||
def sample(self, n_samples=1, seed=None):
|
||||
np.random.seed(seed)
|
||||
w = Maxwell.sample_maxwell(self.a, n_samples)
|
||||
u = np.random.uniform(-1., 1., n_samples)
|
||||
aab = self.a * self.a * self.b
|
||||
return w + 0.25*aab + u*np.sqrt(aab*w)
|
||||
|
||||
def to_xml_element(self, element_name):
|
||||
"""Return XML representation of the Watt distribution
|
||||
|
||||
|
|
@ -588,6 +647,10 @@ class Normal(Univariate):
|
|||
cv.check_greater_than('Normal std_dev', std_dev, 0.0)
|
||||
self._std_dev = std_dev
|
||||
|
||||
def sample(self, n_samples=1, seed=None):
|
||||
np.random.seed(seed)
|
||||
return np.random.normal(self.mean_value, self.std_dev, n_samples)
|
||||
|
||||
def to_xml_element(self, element_name):
|
||||
"""Return XML representation of the Normal distribution
|
||||
|
||||
|
|
@ -694,6 +757,15 @@ class Muir(Univariate):
|
|||
cv.check_greater_than('Muir kt', kt, 0.0)
|
||||
self._kt = kt
|
||||
|
||||
@property
|
||||
def std_dev(self):
|
||||
return np.sqrt(4.*self.e0*self.kt/self.m_rat)
|
||||
|
||||
def sample(self, n_samples=1, seed=None):
|
||||
# Based on LANL report LA-05411-MS
|
||||
np.random.seed(seed)
|
||||
return np.random.normal(self.e0, self.std_dev, n_samples)
|
||||
|
||||
def to_xml_element(self, element_name):
|
||||
"""Return XML representation of the Watt distribution
|
||||
|
||||
|
|
@ -803,6 +875,114 @@ class Tabular(Univariate):
|
|||
cv.check_value('interpolation', interpolation, _INTERPOLATION_SCHEMES)
|
||||
self._interpolation = interpolation
|
||||
|
||||
def cdf(self):
|
||||
if not self.interpolation in ('histogram', 'linear-linear'):
|
||||
raise NotImplementedError('Can only generate CDFs for tabular '
|
||||
'distributions using histogram or '
|
||||
'linear-linear interpolation')
|
||||
c = np.zeros_like(self.x)
|
||||
x = np.asarray(self.x)
|
||||
p = np.asarray(self.p)
|
||||
|
||||
if self.interpolation == 'histogram':
|
||||
c[1:] = p[:-1] * np.diff(x)
|
||||
elif self.interpolation == 'linear-linear':
|
||||
c[1:] = 0.5 * (p[:-1] + p[1:]) * np.diff(x)
|
||||
|
||||
return np.cumsum(c)
|
||||
|
||||
def mean(self):
|
||||
"""Compute the mean of the tabular distribution"""
|
||||
if not self.interpolation in ('histogram', 'linear-linear'):
|
||||
raise NotImplementedError('Can only compute mean for tabular '
|
||||
'distributions using histogram '
|
||||
'or linear-linear interpolation.')
|
||||
if self.interpolation == 'linear-linear':
|
||||
mean = 0.0
|
||||
self.normalize()
|
||||
for i in range(1, len(self.x)):
|
||||
y_min = self.p[i-1]
|
||||
y_max = self.p[i]
|
||||
x_min = self.x[i-1]
|
||||
x_max = self.x[i]
|
||||
|
||||
m = (y_max - y_min) / (x_max - x_min)
|
||||
|
||||
exp_val = (1./3.) * m * (x_max**3 - x_min**3)
|
||||
exp_val += 0.5 * m * x_min * (x_min**2 - x_max**2)
|
||||
exp_val += 0.5 * y_min * (x_max**2 - x_min**2)
|
||||
mean += exp_val
|
||||
|
||||
elif self.interpolation == 'histogram':
|
||||
mean = 0.5 * (self.x[:-1] + self.x[1:])
|
||||
mean *= np.diff(self.cdf())
|
||||
mean = sum(mean)
|
||||
|
||||
return mean
|
||||
|
||||
def normalize(self):
|
||||
"""Normalize the probabilities stored on the distribution"""
|
||||
self.p = np.asarray(self.p) / self.cdf().max()
|
||||
|
||||
def sample(self, n_samples=1, seed=None):
|
||||
if not self.interpolation in ('histogram', 'linear-linear'):
|
||||
raise NotImplementedError('Can only sample tabular distributions '
|
||||
'using histogram or '
|
||||
'linear-linear interpolation')
|
||||
np.random.seed(seed)
|
||||
xi = np.random.rand(n_samples)
|
||||
cdf = self.cdf()
|
||||
cdf /= cdf.max()
|
||||
# always use normalized probabilities when sampling
|
||||
p = self.p / cdf.max()
|
||||
|
||||
# get CDF bins that are above the
|
||||
# sampled values
|
||||
c_i = np.full(n_samples, cdf[0])
|
||||
cdf_idx = np.zeros(n_samples, dtype=int)
|
||||
for i, val in enumerate(cdf[:-1]):
|
||||
mask = xi > val
|
||||
c_i[mask] = val
|
||||
cdf_idx[mask] = i
|
||||
|
||||
# get table values at each index where
|
||||
# the random number is less than the next cdf
|
||||
# entry
|
||||
x_i = self.x[cdf_idx]
|
||||
p_i = self.p[cdf_idx]
|
||||
|
||||
if self.interpolation == 'histogram':
|
||||
# mask where probability is greater than zero
|
||||
pos_mask = p_i > 0.0
|
||||
# probabilities greater than zero are set proportional to the
|
||||
# position of the random numebers in relation to the cdf value
|
||||
p_i[pos_mask] = x_i[pos_mask] + (xi[pos_mask] - c_i[pos_mask]) \
|
||||
/ p_i[pos_mask]
|
||||
# probabilities smaller than zero are set to the random number value
|
||||
p_i[~pos_mask] = x_i[~pos_mask]
|
||||
|
||||
samples_out = p_i
|
||||
|
||||
elif self.interpolation == 'linear-linear':
|
||||
# get variable and probability values for the
|
||||
# next entry
|
||||
x_i1 = self.x[cdf_idx + 1]
|
||||
p_i1 = self.p[cdf_idx + 1]
|
||||
# compute slope between entries
|
||||
m = (p_i1 - p_i) / (x_i1 - x_i)
|
||||
# set values for zero slope
|
||||
zero = m == 0.0
|
||||
m[zero] = x_i[zero] + (xi[zero] - c_i[zero]) / p_i[zero]
|
||||
# set values for non-zero slope
|
||||
non_zero = ~zero
|
||||
quad = np.power(p_i[non_zero], 2) + 2.0 * m[non_zero] * (xi[non_zero] - c_i[non_zero])
|
||||
quad[quad < 0.0] = 0.0
|
||||
m[non_zero] = x_i[non_zero] + (np.sqrt(quad) - p_i[non_zero]) / m[non_zero]
|
||||
samples_out = m
|
||||
|
||||
assert all(samples_out < self.x[-1])
|
||||
return samples_out
|
||||
|
||||
def to_xml_element(self, element_name):
|
||||
"""Return XML representation of the tabular distribution
|
||||
|
||||
|
|
@ -890,6 +1070,9 @@ class Legendre(Univariate):
|
|||
def coefficients(self, coefficients):
|
||||
self._coefficients = np.asarray(coefficients)
|
||||
|
||||
def sample(self, n_samples=1, seed=None):
|
||||
raise NotImplementedError
|
||||
|
||||
def to_xml_element(self, element_name):
|
||||
raise NotImplementedError
|
||||
|
||||
|
|
@ -947,6 +1130,25 @@ class Mixture(Univariate):
|
|||
Iterable, Univariate)
|
||||
self._distribution = distribution
|
||||
|
||||
def cdf(self):
|
||||
return np.insert(np.cumsum(self.probability), 0, 0.0)
|
||||
|
||||
def sample(self, n_samples=1, seed=None):
|
||||
np.random.seed(seed)
|
||||
idx = np.random.choice(self.distribution, n_samples, p=self.probability)
|
||||
|
||||
out = np.zeros_like(idx)
|
||||
for i in np.unique(idx):
|
||||
n_dist_samples = np.count_nonzero(idx == i)
|
||||
samples = self.distribution[i].sample(n_dist_samples)
|
||||
out[idx == i] = samples
|
||||
return out
|
||||
|
||||
def normalize(self):
|
||||
"""Normalize the probabilities stored on the distribution"""
|
||||
norm = sum(self.probability)
|
||||
self.probability = [val / norm for val in self.probability]
|
||||
|
||||
def to_xml_element(self, element_name):
|
||||
"""Return XML representation of the mixture distribution
|
||||
|
||||
|
|
|
|||
2
setup.py
2
setup.py
|
|
@ -32,7 +32,7 @@ kwargs = {
|
|||
# Data files and libraries
|
||||
'package_data': {
|
||||
'openmc.lib': ['libopenmc.{}'.format(suffix)],
|
||||
'openmc.data': ['mass16.txt', 'BREMX.DAT', '*.h5'],
|
||||
'openmc.data': ['mass16.txt', 'BREMX.DAT', 'half_life.json', '*.h5'],
|
||||
'openmc.data.effective_dose': ['*.txt']
|
||||
},
|
||||
|
||||
|
|
|
|||
|
|
@ -14,7 +14,9 @@
|
|||
#include "openmc/settings.h"
|
||||
|
||||
#include "xtensor/xbuilder.hpp"
|
||||
#include "xtensor/xmath.hpp"
|
||||
#include "xtensor/xoperation.hpp"
|
||||
#include "xtensor/xslice.hpp"
|
||||
#include "xtensor/xview.hpp"
|
||||
|
||||
#include <cmath>
|
||||
|
|
@ -44,6 +46,8 @@ vector<unique_ptr<PhotonInteraction>> elements;
|
|||
|
||||
PhotonInteraction::PhotonInteraction(hid_t group)
|
||||
{
|
||||
using namespace xt::placeholders;
|
||||
|
||||
// Set index of element in global vector
|
||||
index_ = data::elements.size();
|
||||
|
||||
|
|
@ -125,6 +129,7 @@ PhotonInteraction::PhotonInteraction(hid_t group)
|
|||
}
|
||||
|
||||
shells_.resize(n_shell);
|
||||
cross_sections_ = xt::zeros<double>({energy_.size(), n_shell});
|
||||
|
||||
// Create mapping from designator to index
|
||||
std::unordered_map<int, int> shell_map;
|
||||
|
|
@ -155,13 +160,15 @@ PhotonInteraction::PhotonInteraction(hid_t group)
|
|||
read_attribute(tgroup, "num_electrons", shell.n_electrons);
|
||||
|
||||
// Read subshell cross section
|
||||
xt::xtensor<double, 1> xs;
|
||||
dset = open_dataset(tgroup, "xs");
|
||||
read_attribute(dset, "threshold_idx", shell.threshold);
|
||||
close_dataset(dset);
|
||||
read_dataset(tgroup, "xs", shell.cross_section);
|
||||
read_dataset(tgroup, "xs", xs);
|
||||
|
||||
auto& xs = shell.cross_section;
|
||||
xs = xt::where(xs > 0.0, xt::log(xs), -500.0);
|
||||
auto cross_section =
|
||||
xt::view(cross_sections_, xt::range(shell.threshold, _), i);
|
||||
cross_section = xt::where(xs > 0, xt::log(xs), 0);
|
||||
|
||||
if (object_exists(tgroup, "transitions")) {
|
||||
// Determine dimensions of transitions
|
||||
|
|
@ -566,18 +573,13 @@ void PhotonInteraction::calculate_xs(Particle& p) const
|
|||
|
||||
// Calculate microscopic photoelectric cross section
|
||||
xs.photoelectric = 0.0;
|
||||
for (const auto& shell : shells_) {
|
||||
// Check threshold of reaction
|
||||
int i_start = shell.threshold;
|
||||
if (i_grid < i_start)
|
||||
continue;
|
||||
const auto& xs_lower = xt::row(cross_sections_, i_grid);
|
||||
const auto& xs_upper = xt::row(cross_sections_, i_grid + 1);
|
||||
|
||||
// Evaluation subshell photoionization cross section
|
||||
xs.photoelectric +=
|
||||
std::exp(shell.cross_section(i_grid - i_start) +
|
||||
f * (shell.cross_section(i_grid + 1 - i_start) -
|
||||
shell.cross_section(i_grid - i_start)));
|
||||
}
|
||||
for (int i = 0; i < xs_upper.size(); ++i)
|
||||
if (xs_lower(i) != 0)
|
||||
xs.photoelectric +=
|
||||
std::exp(xs_lower(i) + f * (xs_upper(i) - xs_lower(i)));
|
||||
|
||||
// Calculate microscopic pair production cross section
|
||||
xs.pair_production = std::exp(
|
||||
|
|
|
|||
|
|
@ -29,6 +29,7 @@
|
|||
|
||||
#include <algorithm> // for max, min, max_element
|
||||
#include <cmath> // for sqrt, exp, log, abs, copysign
|
||||
#include <xtensor/xview.hpp>
|
||||
|
||||
namespace openmc {
|
||||
|
||||
|
|
@ -344,25 +345,26 @@ void sample_photon_reaction(Particle& p)
|
|||
|
||||
// Photoelectric effect
|
||||
double prob_after = prob + micro.photoelectric;
|
||||
|
||||
if (prob_after > cutoff) {
|
||||
// Get grid index, interpolation factor, and bounding subshell
|
||||
// cross sections
|
||||
int i_grid = micro.index_grid;
|
||||
double f = micro.interp_factor;
|
||||
const auto& xs_lower = xt::row(element.cross_sections_, i_grid);
|
||||
const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1);
|
||||
|
||||
for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) {
|
||||
const auto& shell {element.shells_[i_shell]};
|
||||
|
||||
// Get grid index and interpolation factor
|
||||
int i_grid = micro.index_grid;
|
||||
double f = micro.interp_factor;
|
||||
|
||||
// Check threshold of reaction
|
||||
int i_start = shell.threshold;
|
||||
if (i_grid < i_start)
|
||||
if (xs_lower(i_shell) == 0)
|
||||
continue;
|
||||
|
||||
// Evaluation subshell photoionization cross section
|
||||
double xs = std::exp(shell.cross_section(i_grid - i_start) +
|
||||
f * (shell.cross_section(i_grid + 1 - i_start) -
|
||||
shell.cross_section(i_grid - i_start)));
|
||||
// Evaluation subshell photoionization cross section
|
||||
prob += std::exp(
|
||||
xs_lower(i_shell) + f * (xs_upper(i_shell) - xs_lower(i_shell)));
|
||||
|
||||
prob += xs;
|
||||
if (prob > cutoff) {
|
||||
double E_electron = p.E() - shell.binding_energy;
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -38,6 +38,17 @@ def test_remove_nuclide():
|
|||
assert m.nuclides[1].percent == 2.0
|
||||
|
||||
|
||||
def test_remove_elements():
|
||||
"""Test removing elements."""
|
||||
m = openmc.Material()
|
||||
for elem, percent in [('Li', 1.0), ('Be', 1.0)]:
|
||||
m.add_element(elem, percent)
|
||||
m.remove_element('Li')
|
||||
assert len(m.nuclides) == 1
|
||||
assert m.nuclides[0].name == 'Be9'
|
||||
assert m.nuclides[0].percent == 1.0
|
||||
|
||||
|
||||
def test_elements():
|
||||
"""Test adding elements."""
|
||||
m = openmc.Material()
|
||||
|
|
@ -89,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"):
|
||||
|
|
@ -106,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()
|
||||
|
|
@ -118,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()
|
||||
|
|
@ -126,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"):
|
||||
|
|
@ -144,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()
|
||||
|
|
@ -278,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')
|
||||
|
|
@ -330,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.
|
||||
|
|
@ -404,3 +424,30 @@ def test_mix_materials():
|
|||
assert m3.density == pytest.approx(dens3)
|
||||
assert m4.density == pytest.approx(dens4)
|
||||
assert m5.density == pytest.approx(dens5)
|
||||
|
||||
|
||||
def test_activity_of_stable():
|
||||
"""Creates a material with stable isotopes to checks the activity is 0"""
|
||||
m1 = openmc.Material()
|
||||
m1.add_element("Fe", 1)
|
||||
m1.set_density('g/cm3', 1)
|
||||
m1.volume = 1
|
||||
assert m1.activity == 0
|
||||
|
||||
|
||||
def test_activity_of_tritium():
|
||||
"""Checks that 1g of tritium has the correct activity"""
|
||||
m1 = openmc.Material()
|
||||
m1.add_nuclide("H3", 1)
|
||||
m1.set_density('g/cm3', 1)
|
||||
m1.volume = 1
|
||||
assert pytest.approx(m1.activity) == 3.559778e14
|
||||
|
||||
|
||||
def test_activity_of_metastable():
|
||||
"""Checks that 1 mol of a Tc99_m1 nuclides has the correct activity"""
|
||||
m1 = openmc.Material()
|
||||
m1.add_nuclide("Tc99_m1", 1)
|
||||
m1.set_density('g/cm3', 1)
|
||||
m1.volume = 98.9
|
||||
assert pytest.approx(m1.activity, rel=0.001) == 1.93e19
|
||||
|
|
|
|||
76
tests/unit_tests/test_mesh_to_vtk.py
Normal file
76
tests/unit_tests/test_mesh_to_vtk.py
Normal file
|
|
@ -0,0 +1,76 @@
|
|||
import numpy as np
|
||||
from pathlib import Path
|
||||
import pytest
|
||||
|
||||
vtk = pytest.importorskip("vtk")
|
||||
from vtk.util import numpy_support as nps
|
||||
|
||||
import openmc
|
||||
|
||||
|
||||
regular_mesh = openmc.RegularMesh()
|
||||
regular_mesh.lower_left = (0, 0, 0)
|
||||
regular_mesh.upper_right = (1, 1, 1)
|
||||
regular_mesh.dimension = [30, 20, 10]
|
||||
|
||||
rectilinear_mesh = openmc.RectilinearMesh()
|
||||
rectilinear_mesh.x_grid = np.linspace(1, 2, num=30)
|
||||
rectilinear_mesh.y_grid = np.linspace(1, 2, num=30)
|
||||
rectilinear_mesh.z_grid = np.linspace(1, 2, num=30)
|
||||
|
||||
cylinder_mesh = openmc.CylindricalMesh()
|
||||
cylinder_mesh.r_grid = np.linspace(1, 2, num=30)
|
||||
cylinder_mesh.phi_grid = np.linspace(0, np.pi, num=50)
|
||||
cylinder_mesh.z_grid = np.linspace(0, 1, num=30)
|
||||
|
||||
spherical_mesh = openmc.SphericalMesh()
|
||||
spherical_mesh.r_grid = np.linspace(1, 2, num=30)
|
||||
spherical_mesh.phi_grid = np.linspace(0, np.pi, num=50)
|
||||
spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, num=30)
|
||||
|
||||
@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh])
|
||||
def test_write_data_to_vtk(mesh, tmpdir):
|
||||
# BUILD
|
||||
filename = Path(tmpdir) / "out.vtk"
|
||||
|
||||
data = np.random.random(mesh.num_mesh_cells)
|
||||
|
||||
# RUN
|
||||
mesh.write_data_to_vtk(filename=filename, datasets={"label1": data, "label2": data})
|
||||
|
||||
# TEST
|
||||
assert filename.is_file()
|
||||
|
||||
# read file
|
||||
reader = vtk.vtkStructuredGridReader()
|
||||
reader.SetFileName(filename)
|
||||
reader.Update()
|
||||
|
||||
# check name of datasets
|
||||
vtk_grid = reader.GetOutput()
|
||||
array1 = vtk_grid.GetCellData().GetArray(0)
|
||||
array2 = vtk_grid.GetCellData().GetArray(1)
|
||||
|
||||
assert array1.GetName() == "label1"
|
||||
assert array2.GetName() == "label2"
|
||||
|
||||
# check size of datasets
|
||||
assert nps.vtk_to_numpy(array1).size == data.size
|
||||
assert nps.vtk_to_numpy(array2).size == data.size
|
||||
|
||||
@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh])
|
||||
def test_write_data_to_vtk_size_mismatch(mesh):
|
||||
"""Checks that an error is raised when the size of the dataset
|
||||
doesn't match the mesh number of cells
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mesh : openmc.StructuredMesh
|
||||
The mesh to test
|
||||
"""
|
||||
right_size = mesh.num_mesh_cells
|
||||
data = np.random.random(right_size + 1)
|
||||
|
||||
expected_error_msg = "The size of the dataset label should be equal to the number of cells"
|
||||
with pytest.raises(RuntimeError, match=expected_error_msg):
|
||||
mesh.write_data_to_vtk(filename="out.vtk", datasets={"label": data})
|
||||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -26,6 +26,21 @@ def test_discrete():
|
|||
assert d2.p == [1.0]
|
||||
assert len(d2) == 1
|
||||
|
||||
vals = np.array([1.0, 2.0, 3.0])
|
||||
probs = np.array([0.1, 0.7, 0.2])
|
||||
|
||||
exp_mean = (vals * probs).sum()
|
||||
|
||||
d3 = openmc.stats.Discrete(vals, probs)
|
||||
|
||||
# sample discrete distribution
|
||||
n_samples = 1_000_000
|
||||
samples = d3.sample(n_samples, seed=100)
|
||||
# check that the mean of the samples is within 3 std. dev.
|
||||
# of the expected mean
|
||||
std_dev = samples.std() / np.sqrt(n_samples)
|
||||
assert np.abs(exp_mean - samples.mean()) < 3*std_dev
|
||||
|
||||
|
||||
def test_merge_discrete():
|
||||
x1 = [0.0, 1.0, 10.0]
|
||||
|
|
@ -65,8 +80,17 @@ def test_uniform():
|
|||
assert t.p == [1/(b-a), 1/(b-a)]
|
||||
assert t.interpolation == 'histogram'
|
||||
|
||||
exp_mean = 0.5 * (a + b)
|
||||
n_samples = 1_000_000
|
||||
samples = d.sample(n_samples, seed=100)
|
||||
# check that the mean of the samples is within 3 std. dev.
|
||||
# of the expected mean
|
||||
std_dev = samples.std() / np.sqrt(n_samples)
|
||||
assert np.abs(exp_mean - samples.mean()) < 3*std_dev
|
||||
|
||||
|
||||
def test_powerlaw():
|
||||
a, b, n = 10.0, 20.0, 2.0
|
||||
a, b, n = 10.0, 100.0, 2.0
|
||||
d = openmc.stats.PowerLaw(a, b, n)
|
||||
elem = d.to_xml_element('distribution')
|
||||
|
||||
|
|
@ -76,6 +100,17 @@ def test_powerlaw():
|
|||
assert d.n == n
|
||||
assert len(d) == 3
|
||||
|
||||
exp_mean = 100.0 * (n+1) / (n+2)
|
||||
|
||||
# sample power law distribution
|
||||
n_samples = 1_000_000
|
||||
samples = d.sample(n_samples, seed=100)
|
||||
# check that the mean of the samples is within 3 std. dev.
|
||||
# of the expected mean
|
||||
std_dev = samples.std() / np.sqrt(n_samples)
|
||||
assert np.abs(exp_mean - samples.mean()) < 3*std_dev
|
||||
|
||||
|
||||
def test_maxwell():
|
||||
theta = 1.2895e6
|
||||
d = openmc.stats.Maxwell(theta)
|
||||
|
|
@ -85,6 +120,25 @@ def test_maxwell():
|
|||
assert d.theta == theta
|
||||
assert len(d) == 1
|
||||
|
||||
exp_mean = 3/2 * theta
|
||||
|
||||
# sample maxwell distribution
|
||||
n_samples = 1_000_000
|
||||
samples = d.sample(n_samples, seed=100)
|
||||
# check that the mean of the samples is within 3 std. dev.
|
||||
# of the expected mean
|
||||
std_dev = samples.std() / np.sqrt(n_samples)
|
||||
assert np.abs(exp_mean - samples.mean()) < 3*std_dev
|
||||
|
||||
# A second sample with a different seed
|
||||
samples_2 = d.sample(n_samples, seed=200)
|
||||
# check that the mean of the samples is within 3 std. dev.
|
||||
# of the expected mean
|
||||
std_dev = samples_2.std() / np.sqrt(n_samples)
|
||||
assert np.abs(exp_mean - samples_2.mean()) < 3*std_dev
|
||||
|
||||
assert samples_2.mean() != samples.mean()
|
||||
|
||||
|
||||
def test_watt():
|
||||
a, b = 0.965e6, 2.29e-6
|
||||
|
|
@ -96,19 +150,59 @@ def test_watt():
|
|||
assert d.b == b
|
||||
assert len(d) == 2
|
||||
|
||||
# mean value form adapted from
|
||||
# "Prompt-fission-neutron average energy for 238U(n, f ) from
|
||||
# threshold to 200 MeV" Ethvignot et. al.
|
||||
# https://doi.org/10.1016/j.physletb.2003.09.048
|
||||
exp_mean = 3/2 * a + a**2 * b / 4
|
||||
|
||||
# sample Watt distribution
|
||||
n_samples = 1_000_000
|
||||
samples = d.sample(n_samples, seed=100)
|
||||
# check that the mean of the samples is within 3 std. dev.
|
||||
# of the expected mean
|
||||
std_dev = samples.std() / np.sqrt(n_samples)
|
||||
assert np.abs(exp_mean - samples.mean()) < 3*std_dev
|
||||
|
||||
|
||||
def test_tabular():
|
||||
x = [0.0, 5.0, 7.0]
|
||||
p = [0.1, 0.2, 0.05]
|
||||
x = np.array([0.0, 5.0, 7.0])
|
||||
p = np.array([0.1, 0.2, 0.05])
|
||||
d = openmc.stats.Tabular(x, p, 'linear-linear')
|
||||
elem = d.to_xml_element('distribution')
|
||||
|
||||
d = openmc.stats.Tabular.from_xml_element(elem)
|
||||
assert d.x == x
|
||||
assert d.p == p
|
||||
assert all(d.x == x)
|
||||
assert all(d.p == p)
|
||||
assert d.interpolation == 'linear-linear'
|
||||
assert len(d) == len(x)
|
||||
|
||||
# test linear-linear sampling
|
||||
d = openmc.stats.Tabular(x, p)
|
||||
|
||||
n_samples = 100_000
|
||||
samples = d.sample(n_samples, seed=100)
|
||||
diff = np.abs(samples - d.mean())
|
||||
# within_1_sigma = np.count_nonzero(diff < samples.std())
|
||||
# assert within_1_sigma / n_samples >= 0.68
|
||||
within_2_sigma = np.count_nonzero(diff < 2*samples.std())
|
||||
assert within_2_sigma / n_samples >= 0.95
|
||||
within_3_sigma = np.count_nonzero(diff < 3*samples.std())
|
||||
assert within_3_sigma / n_samples >= 0.99
|
||||
|
||||
# test histogram sampling
|
||||
d = openmc.stats.Tabular(x, p, interpolation='histogram')
|
||||
d.normalize()
|
||||
|
||||
samples = d.sample(n_samples, seed=100)
|
||||
diff = np.abs(samples - d.mean())
|
||||
# within_1_sigma = np.count_nonzero(diff < samples.std())
|
||||
# assert within_1_sigma / n_samples >= 0.68
|
||||
within_2_sigma = np.count_nonzero(diff < 2*samples.std())
|
||||
assert within_2_sigma / n_samples >= 0.95
|
||||
within_3_sigma = np.count_nonzero(diff < 3*samples.std())
|
||||
assert within_3_sigma / n_samples >= 0.99
|
||||
|
||||
|
||||
def test_legendre():
|
||||
# Pu239 elastic scattering at 100 keV
|
||||
|
|
@ -251,6 +345,7 @@ def test_point():
|
|||
d = openmc.stats.Point.from_xml_element(elem)
|
||||
assert d.xyz == pytest.approx(p)
|
||||
|
||||
|
||||
def test_normal():
|
||||
mean = 10.0
|
||||
std_dev = 2.0
|
||||
|
|
@ -264,6 +359,18 @@ def test_normal():
|
|||
assert d.std_dev == pytest.approx(std_dev)
|
||||
assert len(d) == 2
|
||||
|
||||
# sample normal distribution
|
||||
n_samples = 10000
|
||||
samples = d.sample(n_samples, seed=100)
|
||||
samples = np.abs(samples - mean)
|
||||
within_1_sigma = np.count_nonzero(samples < std_dev)
|
||||
assert within_1_sigma / n_samples >= 0.68
|
||||
within_2_sigma = np.count_nonzero(samples < 2*std_dev)
|
||||
assert within_2_sigma / n_samples >= 0.95
|
||||
within_3_sigma = np.count_nonzero(samples < 3*std_dev)
|
||||
assert within_3_sigma / n_samples >= 0.99
|
||||
|
||||
|
||||
def test_muir():
|
||||
mean = 10.0
|
||||
mass = 5.0
|
||||
|
|
@ -278,3 +385,14 @@ def test_muir():
|
|||
assert d.m_rat == pytest.approx(mass)
|
||||
assert d.kt == pytest.approx(temp)
|
||||
assert len(d) == 3
|
||||
|
||||
# sample muir distribution
|
||||
n_samples = 10000
|
||||
samples = d.sample(n_samples, seed=100)
|
||||
samples = np.abs(samples - mean)
|
||||
within_1_sigma = np.count_nonzero(samples < d.std_dev)
|
||||
assert within_1_sigma / n_samples >= 0.68
|
||||
within_2_sigma = np.count_nonzero(samples < 2*d.std_dev)
|
||||
assert within_2_sigma / n_samples >= 0.95
|
||||
within_3_sigma = np.count_nonzero(samples < 3*d.std_dev)
|
||||
assert within_3_sigma / n_samples >= 0.99
|
||||
|
|
|
|||
|
|
@ -1,4 +1,5 @@
|
|||
import numpy as np
|
||||
|
||||
import openmc
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue