mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-21 06:25:30 -04:00
Implement DecaySpectrum distribution type and utilize in R2S (#3930)
Co-authored-by: Copilot <copilot@github.com>
This commit is contained in:
Paul Romano
GitHub
parent
f3e1066d46
commit
d56cda2544
13 changed files with 908 additions and 134 deletions
|
|
@ -1058,17 +1058,19 @@ variable and whose sub-elements/attributes are as follows:
|
|||
|
||||
:type:
|
||||
The type of the distribution. Valid options are "uniform", "discrete",
|
||||
"tabular", "maxwell", "watt", and "mixture". The "uniform" option produces
|
||||
variates sampled from a uniform distribution over a finite interval. The
|
||||
"discrete" option produces random variates that can assume a finite number
|
||||
of values (i.e., a distribution characterized by a probability mass function).
|
||||
The "tabular" option produces random variates sampled from a tabulated
|
||||
distribution where the density function is either a histogram or
|
||||
"tabular", "maxwell", "watt", "mixture", and "decay_spectrum". The "uniform"
|
||||
option produces variates sampled from a uniform distribution over a finite
|
||||
interval. The "discrete" option produces random variates that can assume a
|
||||
finite number of values (i.e., a distribution characterized by a probability
|
||||
mass function). The "tabular" option produces random variates sampled from a
|
||||
tabulated distribution where the density function is either a histogram or
|
||||
linearly-interpolated between tabulated points. The "watt" option produces
|
||||
random variates is sampled from a Watt fission spectrum (only used for
|
||||
energies). The "maxwell" option produce variates sampled from a Maxwell
|
||||
fission spectrum (only used for energies). The "mixture" option produces samples
|
||||
from univariate sub-distributions with given probabilities.
|
||||
fission spectrum (only used for energies). The "mixture" option produces
|
||||
samples from univariate sub-distributions with given probabilities. The
|
||||
"decay_spectrum" option produces photon energies sampled from decay photon
|
||||
spectra in a depletion chain (only used for energies).
|
||||
|
||||
*Default*: None
|
||||
|
||||
|
|
@ -1086,6 +1088,10 @@ variable and whose sub-elements/attributes are as follows:
|
|||
:math:`(x,p)` pairs defining the discrete/tabular distribution. All :math:`x`
|
||||
points are given first followed by corresponding :math:`p` points.
|
||||
|
||||
For a "decay_spectrum" distribution, ``parameters`` gives the atom densities
|
||||
in [atom/b-cm] for the nuclides listed in the ``nuclides`` element, in the
|
||||
same order.
|
||||
|
||||
For a "watt" distribution, ``parameters`` should be given as two real numbers
|
||||
:math:`a` and :math:`b` that parameterize the distribution :math:`p(x) dx = c
|
||||
e^{-x/a} \sinh \sqrt{b \, x} dx`.
|
||||
|
|
@ -1115,6 +1121,21 @@ variable and whose sub-elements/attributes are as follows:
|
|||
This sub-element of a ``pair`` element provides information on the
|
||||
corresponding univariate distribution.
|
||||
|
||||
:volume:
|
||||
For a "decay_spectrum" distribution, this attribute specifies the source
|
||||
region volume in cm\ :sup:`3`. It is used together with atom densities to
|
||||
determine the absolute photon emission rate. When a source uses a
|
||||
"decay_spectrum" energy distribution, the source strength is set from this
|
||||
emission rate.
|
||||
|
||||
:nuclides:
|
||||
For a "decay_spectrum" distribution, this element specifies a
|
||||
whitespace-separated list of nuclide names contributing to the decay photon
|
||||
source. The atom densities for these nuclides are given by the ``parameters``
|
||||
element in the same order. Nuclides are resolved against the depletion chain,
|
||||
and nuclides without decay photon spectra do not contribute to the
|
||||
distribution.
|
||||
|
||||
:bias:
|
||||
This optional element specifies a biased distribution for importance sampling.
|
||||
For continuous distributions, the ``bias`` element should contain another
|
||||
|
|
@ -1300,7 +1321,7 @@ The ``<surface_grazing_cutoff>`` element specifies the surface flux cosine cutof
|
|||
``<surface_grazing_ratio>`` Element
|
||||
-----------------------------------
|
||||
|
||||
The ``<surface_grazing_ratio>`` element specifies the surface flux cosine
|
||||
The ``<surface_grazing_ratio>`` element specifies the surface flux cosine
|
||||
substitution ratio.
|
||||
|
||||
*Default*: 0.5
|
||||
|
|
|
|||
|
|
@ -22,6 +22,7 @@ Univariate Probability Distributions
|
|||
openmc.stats.Legendre
|
||||
openmc.stats.Mixture
|
||||
openmc.stats.Normal
|
||||
openmc.stats.DecaySpectrum
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated
|
||||
|
|
|
|||
|
|
@ -101,6 +101,8 @@ extern vector<unique_ptr<ChainNuclide>> chain_nuclides;
|
|||
|
||||
void read_chain_file_xml();
|
||||
|
||||
void free_memory_chain();
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_CHAIN_H
|
||||
|
|
|
|||
|
|
@ -407,6 +407,71 @@ private:
|
|||
double integral_; //!< Integral of distribution
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
// DecaySpectrum — non-owning mixture of decay photon distributions
|
||||
//==============================================================================
|
||||
|
||||
//! Energy distribution formed by mixing multiple decay photon spectra.
|
||||
//!
|
||||
//! Unlike the general Mixture distribution, this class holds non-owning
|
||||
//! pointers to the component distributions (which live in
|
||||
//! data::chain_nuclides). Each component is weighted by the activity
|
||||
//! (atoms * decay_constant) of the corresponding nuclide.
|
||||
|
||||
class DecaySpectrum : public Distribution {
|
||||
public:
|
||||
//============================================================================
|
||||
// Types, aliases
|
||||
|
||||
struct Sample {
|
||||
double energy;
|
||||
double weight;
|
||||
int parent_nuclide;
|
||||
};
|
||||
|
||||
//============================================================================
|
||||
// Constructors
|
||||
|
||||
//! Construct from an XML node containing nuclide names and atom densities.
|
||||
//!
|
||||
//! Reads child ``<nuclide>`` elements with ``name`` and ``density``
|
||||
//! attributes, resolves them against the loaded depletion chain, and
|
||||
//! constructs the mixed distribution.
|
||||
explicit DecaySpectrum(pugi::xml_node node);
|
||||
|
||||
//============================================================================
|
||||
// Methods
|
||||
|
||||
//! Sample a value from the distribution and return the parent nuclide index
|
||||
//! \param seed Pseudorandom number seed pointer
|
||||
//! \return (Sampled energy, sample weight, chain nuclide index)
|
||||
Sample sample_with_parent(uint64_t* seed) const;
|
||||
|
||||
//! Sample a value from the distribution
|
||||
//! \param seed Pseudorandom number seed pointer
|
||||
//! \return (sampled value, sample weight)
|
||||
std::pair<double, double> sample(uint64_t* seed) const override;
|
||||
|
||||
double integral() const override;
|
||||
|
||||
protected:
|
||||
//! Sample a value (unbiased) from the distribution
|
||||
//! \param seed Pseudorandom number seed pointer
|
||||
//! \return Sampled value
|
||||
double sample_unbiased(uint64_t* seed) const override;
|
||||
|
||||
private:
|
||||
//! Initialize decay spectrum sampling data
|
||||
//! \param nuclide_indices Indices of decay photon emitters in
|
||||
//! data::chain_nuclides
|
||||
//! \param atoms Number of atoms for each component.
|
||||
void init(vector<int> nuclide_indices, const vector<double>& atoms);
|
||||
|
||||
vector<int> nuclide_indices_; //!< Indices of emitting nuclides in the chain
|
||||
DiscreteIndex di_; //!< Discrete index for component selection
|
||||
double integral_; //!< Total photon emission rate
|
||||
};
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_DISTRIBUTION_H
|
||||
|
|
|
|||
|
|
@ -1,5 +1,6 @@
|
|||
from __future__ import annotations
|
||||
from collections.abc import Sequence
|
||||
from contextlib import nullcontext
|
||||
import copy
|
||||
from datetime import datetime
|
||||
import json
|
||||
|
|
@ -8,6 +9,7 @@ from pathlib import Path
|
|||
import numpy as np
|
||||
import openmc
|
||||
from . import IndependentOperator, PredictorIntegrator
|
||||
from .chain import Chain
|
||||
from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5
|
||||
from .results import Results
|
||||
from ..checkvalue import PathLike
|
||||
|
|
@ -149,7 +151,7 @@ class R2SManager:
|
|||
photon_time_indices: Sequence[int] | None = None,
|
||||
output_dir: PathLike | None = None,
|
||||
bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
|
||||
chain_file: PathLike | None = None,
|
||||
chain_file: PathLike | Chain | None = None,
|
||||
micro_kwargs: dict | None = None,
|
||||
mat_vol_kwargs: dict | None = None,
|
||||
run_kwargs: dict | None = None,
|
||||
|
|
@ -187,9 +189,10 @@ class R2SManager:
|
|||
Dictionary mapping cell IDs to bounding boxes used for spatial
|
||||
source sampling in cell-based R2S calculations. Required if method
|
||||
is 'cell-based'.
|
||||
chain_file : PathLike, optional
|
||||
Path to the depletion chain XML file to use during activation. If
|
||||
not provided, the default configured chain file will be used.
|
||||
chain_file : PathLike or openmc.deplete.Chain, optional
|
||||
Path to the depletion chain XML file or depletion chain object to
|
||||
use during activation. If not provided, the default configured
|
||||
chain file will be used.
|
||||
micro_kwargs : dict, optional
|
||||
Additional keyword arguments passed to
|
||||
:func:`openmc.deplete.get_microxs_and_flux` during the neutron
|
||||
|
|
@ -216,6 +219,8 @@ class R2SManager:
|
|||
# consistency (different ranks may have slightly different times)
|
||||
stamp = datetime.now().strftime('%Y-%m-%dT%H-%M-%S')
|
||||
output_dir = Path(comm.bcast(f'r2s_{stamp}'))
|
||||
else:
|
||||
output_dir = Path(output_dir)
|
||||
|
||||
# Set run_kwargs for the neutron transport step
|
||||
if micro_kwargs is None:
|
||||
|
|
@ -226,22 +231,35 @@ class R2SManager:
|
|||
operator_kwargs = {}
|
||||
run_kwargs.setdefault('output', False)
|
||||
micro_kwargs.setdefault('run_kwargs', run_kwargs)
|
||||
# If a chain file is provided, prefer it for steps 1 and 2
|
||||
if chain_file is not None:
|
||||
micro_kwargs.setdefault('chain_file', chain_file)
|
||||
operator_kwargs.setdefault('chain_file', chain_file)
|
||||
|
||||
self.step1_neutron_transport(
|
||||
output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs
|
||||
)
|
||||
self.step2_activation(
|
||||
timesteps, source_rates, timestep_units, output_dir / 'activation',
|
||||
operator_kwargs=operator_kwargs
|
||||
)
|
||||
self.step3_photon_transport(
|
||||
photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
|
||||
mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
|
||||
# DecaySpectrum distributions are resolved in the C++ solver using
|
||||
# OPENMC_CHAIN_FILE. If a Chain object was passed, write an XML
|
||||
# representation alongside the R2S outputs.
|
||||
if isinstance(chain_file, Chain):
|
||||
output_dir.mkdir(parents=True, exist_ok=True)
|
||||
chain_path = output_dir / 'chain.xml'
|
||||
if comm.rank == 0:
|
||||
chain_file.export_to_xml(chain_path)
|
||||
comm.barrier()
|
||||
else:
|
||||
chain_path = chain_file
|
||||
|
||||
chain_context = (
|
||||
openmc.config.patch('chain_file', chain_path)
|
||||
if chain_path is not None else nullcontext()
|
||||
)
|
||||
with chain_context:
|
||||
self.step1_neutron_transport(
|
||||
output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs
|
||||
)
|
||||
self.step2_activation(
|
||||
timesteps, source_rates, timestep_units,
|
||||
output_dir / 'activation', operator_kwargs=operator_kwargs
|
||||
)
|
||||
self.step3_photon_transport(
|
||||
photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
|
||||
mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
|
||||
)
|
||||
|
||||
return output_dir
|
||||
|
||||
|
|
@ -516,45 +534,30 @@ class R2SManager:
|
|||
if different_photon_model:
|
||||
photon_cells = self.photon_model.geometry.get_all_cells()
|
||||
|
||||
for time_index in time_indices:
|
||||
# Create decay photon source
|
||||
if self.method == 'mesh-based':
|
||||
self.photon_model.settings.source = \
|
||||
self.get_decay_photon_source_mesh(time_index)
|
||||
else:
|
||||
sources = []
|
||||
results = self.results['depletion_results']
|
||||
for cell, original_mat in zip(self.domains, self.results['activation_materials']):
|
||||
# Skip if the cell is not in the photon model or the
|
||||
# material has changed
|
||||
if different_photon_model:
|
||||
if cell.id not in photon_cells or \
|
||||
# Determine eligible work items upfront (independent of time index).
|
||||
if self.method == 'mesh-based':
|
||||
work_items = self._get_mesh_work_items()
|
||||
else:
|
||||
work_items = []
|
||||
for cell, original_mat in zip(
|
||||
self.domains, self.results['activation_materials']):
|
||||
if different_photon_model:
|
||||
if cell.id not in photon_cells or \
|
||||
cell.fill.id != photon_cells[cell.id].fill.id:
|
||||
continue
|
||||
continue
|
||||
work_items.append((cell, original_mat, bounding_boxes[cell.id]))
|
||||
|
||||
# Get bounding box for the cell
|
||||
bounding_box = bounding_boxes[cell.id]
|
||||
|
||||
# Get activated material composition
|
||||
activated_mat = results[time_index].get_material(str(original_mat.id))
|
||||
|
||||
# Create decay photon source source
|
||||
space = openmc.stats.Box(*bounding_box)
|
||||
energy = activated_mat.get_decay_photon_energy()
|
||||
strength = energy.integral() if energy is not None else 0.0
|
||||
source = openmc.IndependentSource(
|
||||
space=space,
|
||||
energy=energy,
|
||||
particle='photon',
|
||||
strength=strength,
|
||||
constraints={'domains': [cell]}
|
||||
)
|
||||
sources.append(source)
|
||||
self.photon_model.settings.source = sources
|
||||
# Ensure photon transport is enabled in settings
|
||||
self.photon_model.settings.photon_transport = True
|
||||
|
||||
for time_index in time_indices:
|
||||
# Convert time_index (which may be negative) to a normal index
|
||||
if time_index < 0:
|
||||
time_index = len(self.results['depletion_results']) + time_index
|
||||
time_index += len(self.results['depletion_results'])
|
||||
|
||||
# Build decay photon sources and assign to the photon model
|
||||
sources = self._create_photon_sources(time_index, work_items)
|
||||
self.photon_model.settings.source = sources
|
||||
|
||||
# Run photon transport calculation
|
||||
photon_dir = Path(output_dir) / f'time_{time_index}'
|
||||
|
|
@ -567,58 +570,30 @@ class R2SManager:
|
|||
sp.tallies[tally.id] for tally in self.photon_model.tallies
|
||||
]
|
||||
|
||||
def get_decay_photon_source_mesh(
|
||||
self,
|
||||
time_index: int = -1
|
||||
) -> list[openmc.IndependentSource]:
|
||||
"""Create decay photon source for a mesh-based calculation.
|
||||
def _get_mesh_work_items(self):
|
||||
"""Enumerate mesh-based work items across all meshes.
|
||||
|
||||
For each mesh element-material combination across all meshes, an
|
||||
:class:`~openmc.IndependentSource` is created with a
|
||||
:class:`~openmc.stats.Box` spatial distribution based on the bounding
|
||||
box of the material within the mesh element. A material constraint is
|
||||
also applied so that sampled source sites are limited to the correct
|
||||
region.
|
||||
|
||||
When the photon transport model is different from the neutron model, the
|
||||
photon MeshMaterialVolumes is used to determine whether an (element,
|
||||
material) combination exists in the photon model.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
time_index : int, optional
|
||||
Time index for the decay photon source. Default is -1 (last time).
|
||||
Returns a list of (index_mat, mat_id, bbox) tuples for each eligible
|
||||
mesh element--material combination, where index_mat is the index into
|
||||
the activation materials list, mat_id is the material ID, and bbox is
|
||||
the bounding box for that mesh element--material combination.
|
||||
|
||||
Returns
|
||||
-------
|
||||
list of openmc.IndependentSource
|
||||
A list of IndependentSource objects for the decay photons, one for
|
||||
each mesh element-material combination with non-zero source strength.
|
||||
|
||||
list of tuple
|
||||
Each tuple is (index_mat, mat_id, bbox).
|
||||
"""
|
||||
mat_dict = self.neutron_model._get_all_materials()
|
||||
|
||||
# List to hold all sources
|
||||
sources = []
|
||||
|
||||
# Index in the overall list of activated materials
|
||||
index_mat = 0
|
||||
|
||||
# Get various results from previous steps
|
||||
mmv_list = self.results['mesh_material_volumes']
|
||||
materials = self.results['activation_materials']
|
||||
results = self.results['depletion_results']
|
||||
photon_mmv_list = self.results.get('mesh_material_volumes_photon')
|
||||
|
||||
work_items = []
|
||||
index_mat = 0
|
||||
for mesh_idx, mat_vols in enumerate(mmv_list):
|
||||
photon_mat_vols = photon_mmv_list[mesh_idx] \
|
||||
if photon_mmv_list is not None else None
|
||||
|
||||
# Total number of mesh elements for this mesh
|
||||
n_elements = mat_vols.num_elements
|
||||
|
||||
for index_elem in range(n_elements):
|
||||
# Determine which materials exist in the photon model for this element
|
||||
if photon_mat_vols is not None:
|
||||
photon_materials = {
|
||||
mat_id
|
||||
|
|
@ -626,36 +601,77 @@ class R2SManager:
|
|||
if mat_id is not None
|
||||
}
|
||||
|
||||
for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
|
||||
# Skip void volume
|
||||
for mat_id, _, bbox in mat_vols.by_element(
|
||||
index_elem, include_bboxes=True):
|
||||
if mat_id is None:
|
||||
continue
|
||||
|
||||
# Skip if this material doesn't exist in photon model
|
||||
if photon_mat_vols is not None and mat_id not in photon_materials:
|
||||
if photon_mat_vols is not None \
|
||||
and mat_id not in photon_materials:
|
||||
index_mat += 1
|
||||
continue
|
||||
|
||||
# Get activated material composition
|
||||
original_mat = materials[index_mat]
|
||||
activated_mat = results[time_index].get_material(str(original_mat.id))
|
||||
|
||||
# Create decay photon source
|
||||
energy = activated_mat.get_decay_photon_energy()
|
||||
if energy is not None:
|
||||
strength = energy.integral()
|
||||
space = openmc.stats.Box(*bbox)
|
||||
sources.append(openmc.IndependentSource(
|
||||
space=space,
|
||||
energy=energy,
|
||||
particle='photon',
|
||||
strength=strength,
|
||||
constraints={'domains': [mat_dict[mat_id]]}
|
||||
))
|
||||
|
||||
# Increment index of activated material
|
||||
work_items.append((index_mat, mat_id, bbox))
|
||||
index_mat += 1
|
||||
|
||||
return work_items
|
||||
|
||||
def _create_photon_sources(self, time_index, work_items):
|
||||
"""Create decay photon sources for a set of regions.
|
||||
|
||||
Builds :class:`openmc.IndependentSource` objects with
|
||||
:class:`openmc.stats.DecaySpectrum` energy distributions that will be
|
||||
serialized to XML and resolved against the depletion chain by the C++
|
||||
solver.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
time_index : int
|
||||
Index into depletion results.
|
||||
work_items : list of tuple
|
||||
For mesh-based: list of (index_mat, mat_id, bbox).
|
||||
For cell-based: list of (cell, original_mat, bbox).
|
||||
|
||||
Returns
|
||||
-------
|
||||
list of openmc.IndependentSource
|
||||
Photon sources for each activated region.
|
||||
"""
|
||||
step_result = self.results['depletion_results'][time_index]
|
||||
materials = self.results['activation_materials']
|
||||
mesh_based = self.method == 'mesh-based'
|
||||
if mesh_based:
|
||||
mat_dict = self.neutron_model._get_all_materials()
|
||||
|
||||
sources = []
|
||||
for item in work_items:
|
||||
if mesh_based:
|
||||
index_mat, domain_id, bbox = item
|
||||
original_mat = materials[index_mat]
|
||||
domain = mat_dict[domain_id]
|
||||
else:
|
||||
cell, original_mat, bbox = item
|
||||
domain = cell
|
||||
|
||||
activated_mat = step_result.get_material(str(original_mat.id))
|
||||
nuclides = activated_mat.get_nuclide_atom_densities()
|
||||
if not nuclides:
|
||||
continue
|
||||
|
||||
# Eliminate nuclides with zero density
|
||||
nuclides = {nuclide: density for nuclide, density in nuclides.items()
|
||||
if density > 0}
|
||||
|
||||
energy = openmc.stats.DecaySpectrum(nuclides, activated_mat.volume)
|
||||
energy.clip(inplace=True)
|
||||
if not energy.nuclides:
|
||||
continue
|
||||
|
||||
sources.append(openmc.IndependentSource(
|
||||
space=openmc.stats.Box(bbox.lower_left, bbox.upper_right),
|
||||
energy=energy,
|
||||
particle='photon',
|
||||
constraints={'domains': [domain]},
|
||||
))
|
||||
|
||||
return sources
|
||||
|
||||
def load_results(self, path: PathLike):
|
||||
|
|
|
|||
|
|
@ -2,8 +2,10 @@ from __future__ import annotations
|
|||
from abc import ABC, abstractmethod
|
||||
from collections import defaultdict
|
||||
from collections.abc import Iterable, Sequence
|
||||
from functools import cache
|
||||
from math import sqrt, pi, exp, log
|
||||
from numbers import Real
|
||||
from pathlib import Path
|
||||
from warnings import warn
|
||||
|
||||
import lxml.etree as ET
|
||||
|
|
@ -14,6 +16,7 @@ import scipy
|
|||
|
||||
import openmc.checkvalue as cv
|
||||
from openmc.data import atomic_mass, NEUTRON_MASS
|
||||
import openmc.data
|
||||
from .._xml import get_elem_list, get_text
|
||||
from ..mixin import EqualityMixin
|
||||
|
||||
|
|
@ -124,6 +127,8 @@ class Univariate(EqualityMixin, ABC):
|
|||
return Legendre.from_xml_element(elem)
|
||||
elif distribution == 'mixture':
|
||||
return Mixture.from_xml_element(elem)
|
||||
elif distribution == 'decay_spectrum':
|
||||
return DecaySpectrum.from_xml_element(elem)
|
||||
|
||||
@abstractmethod
|
||||
def _sample_unbiased(self, n_samples: int = 1, seed: int | None = None):
|
||||
|
|
@ -2196,6 +2201,310 @@ class Mixture(Univariate):
|
|||
return new_dist
|
||||
|
||||
|
||||
class DecaySpectrum(Univariate):
|
||||
"""Energy distribution from decay photon spectra of a mixture of nuclides.
|
||||
|
||||
This distribution stores nuclide names, their atom densities, and the volume
|
||||
of the region. When written to XML and read by the C++ solver, the nuclide
|
||||
names are resolved against the depletion chain to obtain the decay photon
|
||||
energy spectra and decay constants. The resulting distribution is a mixture
|
||||
of per-nuclide photon spectra weighted by absolute activity. The volume is
|
||||
necessary so that the C++ solver can compute the total photon emission rate
|
||||
in [photons/s], which is used as the source strength.
|
||||
|
||||
.. versionadded:: 0.15.4
|
||||
|
||||
Parameters
|
||||
----------
|
||||
nuclides : dict
|
||||
Dictionary mapping nuclide name (str) to atom density (float) in units
|
||||
of [atom/b-cm].
|
||||
volume : float
|
||||
Volume of the source region in [cm³]. Used together with atom densities
|
||||
to compute the absolute photon emission rate.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
nuclides : dict
|
||||
Dictionary mapping nuclide name to atom density in [atom/b-cm].
|
||||
volume : float
|
||||
Volume of the source region in [cm³].
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, nuclides: dict[str, float], volume: float):
|
||||
super().__init__(bias=None)
|
||||
self._dist_cache = None
|
||||
self._dist_cache_key = None
|
||||
self.nuclides = nuclides
|
||||
self.volume = volume
|
||||
|
||||
def __len__(self):
|
||||
return len(self.nuclides)
|
||||
|
||||
@property
|
||||
def nuclides(self):
|
||||
return self._nuclides
|
||||
|
||||
@nuclides.setter
|
||||
def nuclides(self, nuclides):
|
||||
cv.check_type('nuclides', nuclides, dict)
|
||||
for name, density in nuclides.items():
|
||||
cv.check_type('nuclide name', name, str)
|
||||
cv.check_type(f'atom density for {name}', density, Real)
|
||||
cv.check_greater_than(f'atom density for {name}', density, 0.0, True)
|
||||
self._nuclides = dict(nuclides)
|
||||
self._dist_cache = None
|
||||
self._dist_cache_key = None
|
||||
|
||||
@property
|
||||
def volume(self):
|
||||
return self._volume
|
||||
|
||||
@volume.setter
|
||||
def volume(self, volume):
|
||||
cv.check_type('volume', volume, Real)
|
||||
cv.check_greater_than('volume', volume, 0.0)
|
||||
self._volume = float(volume)
|
||||
self._dist_cache = None
|
||||
self._dist_cache_key = None
|
||||
|
||||
@staticmethod
|
||||
def _chain_file_cache_key():
|
||||
"""Return a hashable key for the active depletion chain."""
|
||||
chain_file = openmc.config.get('chain_file')
|
||||
if chain_file is None:
|
||||
return None
|
||||
|
||||
path = Path(chain_file).resolve()
|
||||
try:
|
||||
stat = path.stat()
|
||||
except OSError:
|
||||
return (path, None, None)
|
||||
return (path, stat.st_mtime, stat.st_size)
|
||||
|
||||
def to_distribution(self):
|
||||
"""Convert to a concrete distribution using decay chain data.
|
||||
|
||||
Builds a combined photon energy distribution by looking up each nuclide
|
||||
in the depletion chain via :func:`openmc.data.decay_photon_energy` and
|
||||
weighting by absolute atom count (``density * 1e24 * volume``). The
|
||||
result is cached on the object; the cache is invalidated automatically
|
||||
when :attr:`nuclides` or :attr:`volume` are reassigned.
|
||||
|
||||
Requires ``openmc.config['chain_file']`` to be set.
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.stats.Univariate or None
|
||||
Combined photon energy distribution, or ``None`` if no nuclide in
|
||||
:attr:`nuclides` has a photon source in the chain.
|
||||
|
||||
"""
|
||||
chain_key = self._chain_file_cache_key()
|
||||
if self._dist_cache is not None and self._dist_cache_key == chain_key:
|
||||
return self._dist_cache
|
||||
|
||||
dists = []
|
||||
weights = []
|
||||
for name, density in self.nuclides.items():
|
||||
dist = openmc.data.decay_photon_energy(name)
|
||||
if dist is not None:
|
||||
dists.append(dist)
|
||||
weights.append(density * 1e24 * self.volume)
|
||||
|
||||
if not dists:
|
||||
return None
|
||||
|
||||
self._dist_cache = combine_distributions(dists, weights)
|
||||
self._dist_cache_key = chain_key
|
||||
return self._dist_cache
|
||||
|
||||
def to_xml_element(self, element_name: str):
|
||||
"""Return XML representation of the decay photon distribution
|
||||
|
||||
Parameters
|
||||
----------
|
||||
element_name : str
|
||||
XML element name
|
||||
|
||||
Returns
|
||||
-------
|
||||
element : lxml.etree._Element
|
||||
XML element containing decay photon distribution data
|
||||
|
||||
"""
|
||||
element = ET.Element(element_name)
|
||||
element.set("type", "decay_spectrum")
|
||||
element.set("volume", str(self.volume))
|
||||
nuclides = ET.SubElement(element, "nuclides")
|
||||
nuclides.text = ' '.join(self.nuclides)
|
||||
parameters = ET.SubElement(element, "parameters")
|
||||
parameters.text = ' '.join(str(density) for density in self.nuclides.values())
|
||||
return element
|
||||
|
||||
@classmethod
|
||||
def from_xml_element(cls, elem: ET.Element):
|
||||
"""Generate decay photon distribution from an XML element
|
||||
|
||||
Parameters
|
||||
----------
|
||||
elem : lxml.etree._Element
|
||||
XML element
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.stats.DecaySpectrum
|
||||
Decay photon distribution generated from XML element
|
||||
|
||||
"""
|
||||
volume = float(elem.get('volume'))
|
||||
names = get_elem_list(elem, 'nuclides', str)
|
||||
densities = get_elem_list(elem, 'parameters', float)
|
||||
nuclides = dict(zip(names, densities))
|
||||
return cls(nuclides, volume)
|
||||
|
||||
def _sample_unbiased(self, n_samples=1, seed=None):
|
||||
dist = self.to_distribution()
|
||||
if dist is None:
|
||||
raise RuntimeError(
|
||||
"DecaySpectrum._sample_unbiased requires chain data but none "
|
||||
"was found. Ensure openmc.config['chain_file'] is set and the "
|
||||
"chain contains photon sources for the nuclides present."
|
||||
)
|
||||
return dist.sample(n_samples, seed)[0]
|
||||
|
||||
def integral(self):
|
||||
"""Return integral of the distribution
|
||||
|
||||
Returns the total photon emission rate in [photons/s] by delegating to
|
||||
:meth:`to_distribution`. Returns ``0.0`` when no chain data is
|
||||
available (e.g., ``openmc.config['chain_file']`` is not set).
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Total photon emission rate in [photons/s], or ``0.0`` if chain
|
||||
data is unavailable.
|
||||
"""
|
||||
try:
|
||||
dist = self.to_distribution()
|
||||
except Exception:
|
||||
return 0.0
|
||||
if dist is None:
|
||||
return 0.0
|
||||
return dist.integral()
|
||||
|
||||
@staticmethod
|
||||
@cache
|
||||
def _photon_integral(nuclide: str, chain_key) -> float | None:
|
||||
"""Return the per-atom photon emission integral for a nuclide"""
|
||||
dist = openmc.data.decay_photon_energy(nuclide)
|
||||
return dist.integral() if dist is not None else None
|
||||
|
||||
def clip(self, tolerance: float = 1e-9, inplace: bool = False):
|
||||
"""Remove nuclides with negligible contribution to photon emission.
|
||||
|
||||
Nuclides that are stable or have no photon source in the depletion
|
||||
chain are removed unconditionally. The remaining nuclides are ranked
|
||||
by their photon emission rate (proportional to
|
||||
``atom_density * decay_constant * photon_yield``) and the least
|
||||
important are discarded until the cumulative discarded fraction of the
|
||||
total emission rate exceeds *tolerance*.
|
||||
|
||||
Requires ``openmc.config['chain_file']`` to be set.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
tolerance : float
|
||||
Maximum fraction of total photon emission rate that may be
|
||||
discarded.
|
||||
inplace : bool
|
||||
Whether to modify the current object in-place or return a new one.
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.stats.DecaySpectrum
|
||||
Distribution with negligible nuclides removed.
|
||||
|
||||
"""
|
||||
# Compute per-nuclide emission rate; drop non-emitters
|
||||
emitting_names = []
|
||||
emitting_densities = []
|
||||
rates = []
|
||||
chain_key = self._chain_file_cache_key()
|
||||
for name, density in self.nuclides.items():
|
||||
integral = DecaySpectrum._photon_integral(name, chain_key)
|
||||
if integral is None:
|
||||
continue
|
||||
emitting_names.append(name)
|
||||
emitting_densities.append(density)
|
||||
rates.append(density * self.volume * integral)
|
||||
|
||||
if not emitting_names:
|
||||
new_nuclides = {}
|
||||
else:
|
||||
indices = _intensity_clip(rates, tolerance=tolerance)
|
||||
new_nuclides = {
|
||||
emitting_names[i]: emitting_densities[i] for i in indices
|
||||
}
|
||||
|
||||
if inplace:
|
||||
self._nuclides = new_nuclides
|
||||
self._dist_cache = None
|
||||
self._dist_cache_key = None
|
||||
return self
|
||||
return type(self)(new_nuclides, self.volume)
|
||||
|
||||
@property
|
||||
def support(self):
|
||||
return (0.0, np.inf)
|
||||
|
||||
def evaluate(self, x):
|
||||
"""Evaluate the probability density at a given value.
|
||||
|
||||
Delegates to the combined distribution built from chain data. Raises
|
||||
``NotImplementedError`` if the combined distribution is a
|
||||
:class:`~openmc.stats.Mixture` (which does not support
|
||||
``evaluate()``).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : float
|
||||
Value at which to evaluate the PDF.
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Probability density at *x*.
|
||||
"""
|
||||
dist = self.to_distribution()
|
||||
if dist is None:
|
||||
raise RuntimeError(
|
||||
"DecaySpectrum.evaluate requires chain data. Ensure "
|
||||
"openmc.config['chain_file'] is set."
|
||||
)
|
||||
return dist.evaluate(x)
|
||||
|
||||
def mean(self):
|
||||
"""Return the mean of the distribution.
|
||||
|
||||
Delegates to the combined distribution built from chain data.
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Mean photon energy in [eV].
|
||||
"""
|
||||
dist = self.to_distribution()
|
||||
if dist is None:
|
||||
raise RuntimeError(
|
||||
"DecaySpectrum.mean requires chain data. Ensure "
|
||||
"openmc.config['chain_file'] is set."
|
||||
)
|
||||
return dist.mean()
|
||||
|
||||
|
||||
def combine_distributions(
|
||||
dists: Sequence[Discrete | Tabular | Mixture],
|
||||
probs: Sequence[float]
|
||||
|
|
|
|||
|
|
@ -98,6 +98,8 @@ vector<unique_ptr<ChainNuclide>> chain_nuclides;
|
|||
|
||||
void read_chain_file_xml()
|
||||
{
|
||||
free_memory_chain();
|
||||
|
||||
char* chain_file_path = std::getenv("OPENMC_CHAIN_FILE");
|
||||
if (!chain_file_path) {
|
||||
return;
|
||||
|
|
@ -120,4 +122,10 @@ void read_chain_file_xml()
|
|||
}
|
||||
}
|
||||
|
||||
void free_memory_chain()
|
||||
{
|
||||
data::chain_nuclides.clear();
|
||||
data::chain_nuclide_map.clear();
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -7,7 +7,9 @@
|
|||
#include <numeric> // for accumulate
|
||||
#include <stdexcept> // for runtime_error
|
||||
#include <string> // for string, stod
|
||||
#include <unordered_set>
|
||||
|
||||
#include "openmc/chain.h"
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/error.h"
|
||||
#include "openmc/math_functions.h"
|
||||
|
|
@ -15,6 +17,10 @@
|
|||
#include "openmc/random_lcg.h"
|
||||
#include "openmc/xml_interface.h"
|
||||
|
||||
namespace {
|
||||
std::unordered_set<std::string> decay_spectrum_missing_chain_nuclides;
|
||||
}
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -758,6 +764,8 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
|
|||
dist = UPtrDist {new Tabular(node)};
|
||||
} else if (type == "mixture") {
|
||||
dist = UPtrDist {new Mixture(node)};
|
||||
} else if (type == "decay_spectrum") {
|
||||
dist = UPtrDist {new DecaySpectrum(node)};
|
||||
} else if (type == "muir") {
|
||||
openmc::fatal_error(
|
||||
"'muir' distributions are now specified using the openmc.stats.muir() "
|
||||
|
|
@ -768,4 +776,120 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
|
|||
return dist;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// DecaySpectrum implementation
|
||||
//==============================================================================
|
||||
|
||||
DecaySpectrum::DecaySpectrum(pugi::xml_node node)
|
||||
{
|
||||
// Read the region volume [cm^3] needed for absolute emission rate
|
||||
if (!check_for_node(node, "volume"))
|
||||
fatal_error("DecaySpectrum: 'volume' attribute is required.");
|
||||
double volume = std::stod(get_node_value(node, "volume"));
|
||||
|
||||
// Read nuclide names and atom densities from XML
|
||||
vector<int> nuclide_indices;
|
||||
vector<double> atoms;
|
||||
auto names = get_node_array<std::string>(node, "nuclides");
|
||||
auto densities = get_node_array<double>(node, "parameters");
|
||||
if (names.size() != densities.size()) {
|
||||
fatal_error("DecaySpectrum nuclides and parameters must have the same "
|
||||
"length.");
|
||||
}
|
||||
|
||||
for (size_t i = 0; i < names.size(); ++i) {
|
||||
const auto& name = names[i];
|
||||
double density = densities[i];
|
||||
|
||||
// Look up nuclide in the depletion chain
|
||||
auto it = data::chain_nuclide_map.find(name);
|
||||
if (it == data::chain_nuclide_map.end()) {
|
||||
if (decay_spectrum_missing_chain_nuclides.insert(name).second) {
|
||||
warning("Nuclide '" + name +
|
||||
"' appears in a DecaySpectrum source but is not present in "
|
||||
"the depletion chain; it will be ignored.");
|
||||
}
|
||||
continue;
|
||||
}
|
||||
|
||||
int nuclide_index = it->second;
|
||||
const auto& chain_nuc = data::chain_nuclides[nuclide_index];
|
||||
const Distribution* photon_dist = chain_nuc->photon_energy();
|
||||
if (!photon_dist)
|
||||
continue;
|
||||
|
||||
// Skip non-positive densities and warn if negative
|
||||
if (density <= 0.0) {
|
||||
if (density < 0.0) {
|
||||
warning("Nuclide '" + name +
|
||||
"' has a negative density in a DecaySpectrum source; it will "
|
||||
"be ignored.");
|
||||
}
|
||||
continue;
|
||||
}
|
||||
|
||||
// atoms = density [atom/b-cm] * 1e24 [b/cm^2] * volume [cm^3]
|
||||
double atoms_i = density * 1.0e24 * volume;
|
||||
|
||||
nuclide_indices.push_back(nuclide_index);
|
||||
atoms.push_back(atoms_i);
|
||||
}
|
||||
|
||||
init(std::move(nuclide_indices), atoms);
|
||||
}
|
||||
|
||||
void DecaySpectrum::init(
|
||||
vector<int> nuclide_indices, const vector<double>& atoms)
|
||||
{
|
||||
if (nuclide_indices.size() != atoms.size()) {
|
||||
fatal_error("DecaySpectrum nuclide index and atoms arrays must have "
|
||||
"the same length.");
|
||||
}
|
||||
|
||||
vector<double> probs;
|
||||
probs.reserve(nuclide_indices.size());
|
||||
for (size_t i = 0; i < nuclide_indices.size(); ++i) {
|
||||
// Distribution integral is in [photons/s/atom]; multiplying by atoms gives
|
||||
// the total emission rate [photons/s] for this nuclide.
|
||||
const auto* dist =
|
||||
data::chain_nuclides[nuclide_indices[i]]->photon_energy();
|
||||
probs.push_back(atoms[i] * dist->integral());
|
||||
}
|
||||
|
||||
nuclide_indices_ = std::move(nuclide_indices);
|
||||
integral_ = std::accumulate(probs.begin(), probs.end(), 0.0);
|
||||
if (nuclide_indices_.empty() || integral_ <= 0.0) {
|
||||
fatal_error("DecaySpectrum source did not resolve any nuclides with decay "
|
||||
"photon spectra and positive atom densities. Ensure "
|
||||
"OPENMC_CHAIN_FILE is set and matches the nuclides in the "
|
||||
"source definition.");
|
||||
}
|
||||
di_.assign(probs);
|
||||
}
|
||||
|
||||
DecaySpectrum::Sample DecaySpectrum::sample_with_parent(uint64_t* seed) const
|
||||
{
|
||||
size_t idx = di_.sample(seed);
|
||||
int parent_nuclide = nuclide_indices_[idx];
|
||||
const auto* dist = data::chain_nuclides[parent_nuclide]->photon_energy();
|
||||
auto [energy, weight] = dist->sample(seed);
|
||||
return {energy, weight, parent_nuclide};
|
||||
}
|
||||
|
||||
std::pair<double, double> DecaySpectrum::sample(uint64_t* seed) const
|
||||
{
|
||||
auto sample = sample_with_parent(seed);
|
||||
return {sample.energy, sample.weight};
|
||||
}
|
||||
|
||||
double DecaySpectrum::integral() const
|
||||
{
|
||||
return integral_;
|
||||
}
|
||||
|
||||
double DecaySpectrum::sample_unbiased(uint64_t* seed) const
|
||||
{
|
||||
return sample_with_parent(seed).energy;
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -2,6 +2,7 @@
|
|||
|
||||
#include "openmc/bank.h"
|
||||
#include "openmc/capi.h"
|
||||
#include "openmc/chain.h"
|
||||
#include "openmc/cmfd_solver.h"
|
||||
#include "openmc/collision_track.h"
|
||||
#include "openmc/constants.h"
|
||||
|
|
@ -44,6 +45,7 @@ void free_memory()
|
|||
free_memory_photon();
|
||||
free_memory_settings();
|
||||
free_memory_thermal();
|
||||
free_memory_chain();
|
||||
library_clear();
|
||||
nuclides_clear();
|
||||
free_memory_source();
|
||||
|
|
|
|||
|
|
@ -405,6 +405,10 @@ bool read_model_xml()
|
|||
write_message(
|
||||
fmt::format("Reading model XML file '{}' ...", model_filename), 5);
|
||||
|
||||
// Read chain data before settings so DecaySpectrum source distributions can
|
||||
// resolve nuclides while sources are constructed.
|
||||
read_chain_file_xml();
|
||||
|
||||
read_settings_xml(settings_root);
|
||||
|
||||
// If other XML files are present, display warning
|
||||
|
|
@ -420,9 +424,6 @@ bool read_model_xml()
|
|||
}
|
||||
}
|
||||
|
||||
// Read data from chain file
|
||||
read_chain_file_xml();
|
||||
|
||||
// Read materials and cross sections
|
||||
if (!check_for_node(root, "materials")) {
|
||||
fatal_error(fmt::format(
|
||||
|
|
@ -475,14 +476,15 @@ bool read_model_xml()
|
|||
|
||||
void read_separate_xml_files()
|
||||
{
|
||||
// Read chain data before settings so DecaySpectrum source distributions can
|
||||
// resolve nuclides while sources are constructed.
|
||||
read_chain_file_xml();
|
||||
|
||||
read_settings_xml();
|
||||
if (settings::run_mode != RunMode::PLOTTING) {
|
||||
read_cross_sections_xml();
|
||||
}
|
||||
|
||||
// Read data from chain file
|
||||
read_chain_file_xml();
|
||||
|
||||
read_materials_xml();
|
||||
read_geometry_xml();
|
||||
|
||||
|
|
|
|||
|
|
@ -354,6 +354,19 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node)
|
|||
if (check_for_node(node, "energy")) {
|
||||
pugi::xml_node node_dist = node.child("energy");
|
||||
energy_ = distribution_from_xml(node_dist);
|
||||
|
||||
// For decay photon sources, use the absolute photon emission rate in
|
||||
// [photons/s] as the source strength
|
||||
if (dynamic_cast<DecaySpectrum*>(energy_.get())) {
|
||||
if (strength_ != 1.0) {
|
||||
warning(fmt::format(
|
||||
"Source strength of {} is ignored because the source uses a "
|
||||
"DecaySpectrum energy distribution. The source strength will be "
|
||||
"set from the DecaySpectrum emission rate.",
|
||||
strength_));
|
||||
}
|
||||
strength_ = energy_->integral();
|
||||
}
|
||||
} else {
|
||||
// Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
|
||||
energy_ = UPtrDist {new Watt(0.988e6, 2.249e-6)};
|
||||
|
|
@ -414,6 +427,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
|
|||
// Check for monoenergetic source above maximum particle energy
|
||||
auto p = particle_.transport_index();
|
||||
auto energy_ptr = dynamic_cast<Discrete*>(energy_.get());
|
||||
auto decay_spectrum = dynamic_cast<DecaySpectrum*>(energy_.get());
|
||||
if (energy_ptr) {
|
||||
auto energies =
|
||||
tensor::Tensor<double>(energy_ptr->x().data(), energy_ptr->x().size());
|
||||
|
|
@ -424,10 +438,18 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
|
|||
}
|
||||
|
||||
while (true) {
|
||||
// Sample energy spectrum
|
||||
auto [E, E_wgt_temp] = energy_->sample(seed);
|
||||
site.E = E;
|
||||
E_wgt = E_wgt_temp;
|
||||
// Sample energy spectrum. For decay photon sources, also get the parent
|
||||
// nuclide index to store in the source site for tallying purposes.
|
||||
if (decay_spectrum) {
|
||||
auto sample = decay_spectrum->sample_with_parent(seed);
|
||||
site.E = sample.energy;
|
||||
E_wgt = sample.weight;
|
||||
site.parent_nuclide = sample.parent_nuclide;
|
||||
} else {
|
||||
auto [E, E_wgt_temp] = energy_->sample(seed);
|
||||
site.E = E;
|
||||
E_wgt = E_wgt_temp;
|
||||
}
|
||||
|
||||
// Resample if energy falls above maximum particle energy
|
||||
if (site.E < data::energy_max[p] &&
|
||||
|
|
|
|||
|
|
@ -1,5 +1,6 @@
|
|||
from collections import Counter
|
||||
from math import pi
|
||||
from pathlib import Path
|
||||
|
||||
import openmc
|
||||
import openmc.lib
|
||||
|
|
@ -91,6 +92,74 @@ def test_point_cloud_strengths(run_in_tmpdir, sphere_box_model):
|
|||
assert sampled_strength == expected_strength, f'Strength incorrect for {positions[i]}'
|
||||
|
||||
|
||||
def test_decay_spectrum_parent_nuclide(run_in_tmpdir):
|
||||
chain_file = Path('chain_decay_spectrum_parent.xml')
|
||||
chain_file.write_text("""<?xml version="1.0"?>
|
||||
<depletion_chain>
|
||||
<nuclide name="ParentA" decay_modes="0" reactions="0" half_life="1.0">
|
||||
<source type="discrete" particle="photon">
|
||||
<parameters>1000000.0 1.0</parameters>
|
||||
</source>
|
||||
</nuclide>
|
||||
<nuclide name="ParentB" decay_modes="0" reactions="0" half_life="1.0">
|
||||
<source type="discrete" particle="photon">
|
||||
<parameters>2000000.0 1.0</parameters>
|
||||
</source>
|
||||
</nuclide>
|
||||
</depletion_chain>
|
||||
""")
|
||||
|
||||
inner_sphere = openmc.Sphere(r=10.0)
|
||||
outer_sphere = openmc.Sphere(r=20.0, boundary_type='vacuum')
|
||||
|
||||
shell_mat = openmc.Material()
|
||||
shell_mat.add_nuclide('H1', 1.0)
|
||||
shell_mat.set_density('atom/b-cm', 1.0e-12)
|
||||
|
||||
void_cell = openmc.Cell(region=-inner_sphere)
|
||||
shell_cell = openmc.Cell(fill=shell_mat, region=+inner_sphere & -outer_sphere)
|
||||
|
||||
model = openmc.Model()
|
||||
model.geometry = openmc.Geometry([void_cell, shell_cell])
|
||||
model.materials = [shell_mat]
|
||||
model.settings.run_mode = 'fixed source'
|
||||
model.settings.photon_transport = True
|
||||
model.settings.particles = 1000
|
||||
model.settings.batches = 5
|
||||
model.settings.source = openmc.IndependentSource(
|
||||
particle='photon',
|
||||
space=openmc.stats.Point((0.0, 0.0, 0.0)),
|
||||
energy=openmc.stats.DecaySpectrum(
|
||||
{'ParentA': 1.0, 'ParentB': 1.0},
|
||||
volume=1.0
|
||||
)
|
||||
)
|
||||
|
||||
tally = openmc.Tally()
|
||||
tally.filters = [
|
||||
openmc.CellFilter([void_cell]),
|
||||
openmc.ParticleFilter(['photon']),
|
||||
openmc.EnergyFilter([0.0, 1.5e6, 2.5e6]),
|
||||
openmc.ParentNuclideFilter(['ParentA', 'ParentB'])
|
||||
]
|
||||
tally.scores = ['flux']
|
||||
model.tallies = [tally]
|
||||
|
||||
with openmc.config.patch('chain_file', chain_file):
|
||||
sp_filename = model.run()
|
||||
|
||||
with openmc.StatePoint(sp_filename) as sp:
|
||||
tally_out = sp.tallies[tally.id]
|
||||
mean = tally_out.get_reshaped_data('mean').squeeze()
|
||||
|
||||
assert mean.shape == (2, 2)
|
||||
assert mean[0, 0] > 0.0
|
||||
assert mean[1, 1] > 0.0
|
||||
assert mean[0, 1] == 0.0
|
||||
assert mean[1, 0] == 0.0
|
||||
assert np.count_nonzero(mean) == 2
|
||||
|
||||
|
||||
def test_source_file():
|
||||
filename = 'source.h5'
|
||||
src = openmc.FileSource(path=filename)
|
||||
|
|
|
|||
|
|
@ -1,10 +1,11 @@
|
|||
from math import pi
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import pytest
|
||||
import openmc
|
||||
import openmc.stats
|
||||
from openmc.stats.univariate import _INTERPOLATION_SCHEMES
|
||||
from openmc.stats.univariate import _INTERPOLATION_SCHEMES, DecaySpectrum
|
||||
from scipy.integrate import trapezoid
|
||||
|
||||
from tests.unit_tests import assert_sample_mean
|
||||
|
|
@ -1013,3 +1014,135 @@ def test_fusion_spectrum_invalid():
|
|||
# Temperature above 100 keV should raise an error
|
||||
with pytest.raises(ValueError):
|
||||
openmc.stats.fusion_neutron_spectrum(101e3, 'DT')
|
||||
|
||||
|
||||
@pytest.fixture(autouse=False)
|
||||
def decay_spectrum_chain():
|
||||
"""Set chain_file for the duration of a test and clear the _photon_integral
|
||||
cache so results from a different chain don't bleed across tests."""
|
||||
CHAIN_FILE = (Path(__file__).parents[1] / 'chain_simple.xml').resolve()
|
||||
DecaySpectrum._photon_integral.cache_clear()
|
||||
with openmc.config.patch('chain_file', CHAIN_FILE):
|
||||
yield
|
||||
DecaySpectrum._photon_integral.cache_clear()
|
||||
|
||||
|
||||
def test_decay_spectrum_construction():
|
||||
nuclides = {'I135': 1.5e-3, 'Xe135': 8.2e-4}
|
||||
d = openmc.stats.DecaySpectrum(nuclides, volume=100.0)
|
||||
assert d.nuclides == nuclides
|
||||
assert d.volume == pytest.approx(100.0)
|
||||
assert len(d) == 2
|
||||
|
||||
|
||||
def test_decay_spectrum_validation():
|
||||
# nuclides must be a dict
|
||||
with pytest.raises(TypeError):
|
||||
openmc.stats.DecaySpectrum(['I135'], volume=1.0)
|
||||
|
||||
# densities must be > 0
|
||||
with pytest.raises(ValueError):
|
||||
openmc.stats.DecaySpectrum({'I135': -1.0}, volume=1.0)
|
||||
|
||||
# volume must be > 0
|
||||
with pytest.raises(ValueError):
|
||||
openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=-1.0)
|
||||
|
||||
with pytest.raises(ValueError):
|
||||
openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=0.0)
|
||||
|
||||
|
||||
def test_decay_spectrum_xml_roundtrip():
|
||||
nuclides = {'I135': 1.5e-3, 'Xe135': 8.2e-4}
|
||||
d = openmc.stats.DecaySpectrum(nuclides, volume=100.0)
|
||||
|
||||
elem = d.to_xml_element('energy')
|
||||
assert elem.get('type') == 'decay_spectrum'
|
||||
assert float(elem.get('volume')) == pytest.approx(100.0)
|
||||
assert elem.findtext('nuclides').split() == list(nuclides)
|
||||
assert [float(x) for x in elem.findtext('parameters').split()] == pytest.approx(
|
||||
list(nuclides.values()))
|
||||
|
||||
# Round-trip via DecaySpectrum.from_xml_element
|
||||
d2 = openmc.stats.DecaySpectrum.from_xml_element(elem)
|
||||
assert d2.nuclides == nuclides
|
||||
assert d2.volume == pytest.approx(100.0)
|
||||
|
||||
# Round-trip via the Univariate dispatcher
|
||||
d3 = openmc.stats.Univariate.from_xml_element(elem)
|
||||
assert isinstance(d3, openmc.stats.DecaySpectrum)
|
||||
assert d3 == d
|
||||
|
||||
|
||||
def test_decay_spectrum_to_distribution(decay_spectrum_chain):
|
||||
# Single emitting nuclide -> concrete distribution, not None
|
||||
d = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=10.0)
|
||||
dist = d.to_distribution()
|
||||
assert dist is not None
|
||||
|
||||
# Result is cached on second call
|
||||
dist2 = d.to_distribution()
|
||||
assert dist2 is dist
|
||||
|
||||
# Nuclide with no photon source -> None
|
||||
d_stable = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0)
|
||||
assert d_stable.to_distribution() is None
|
||||
|
||||
# Mixture of emitters -> non-None combined distribution
|
||||
d_mix = openmc.stats.DecaySpectrum(
|
||||
{'I135': 1e-3, 'Xe135': 5e-4}, volume=10.0
|
||||
)
|
||||
dist_mix = d_mix.to_distribution()
|
||||
assert dist_mix is not None
|
||||
|
||||
|
||||
def test_decay_spectrum_integral(decay_spectrum_chain):
|
||||
# For an emitting nuclide, integral should be > 0
|
||||
d = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=10.0)
|
||||
assert d.integral() > 0.0
|
||||
|
||||
# Proportional to density: doubling density doubles integral
|
||||
d2 = openmc.stats.DecaySpectrum({'I135': 2e-3}, volume=10.0)
|
||||
assert d2.integral() == pytest.approx(2.0 * d.integral())
|
||||
|
||||
# Proportional to volume
|
||||
d3 = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=20.0)
|
||||
assert d3.integral() == pytest.approx(2.0 * d.integral())
|
||||
|
||||
# Pure non-emitter -> 0.0
|
||||
d_stable = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0)
|
||||
assert d_stable.integral() == pytest.approx(0.0)
|
||||
|
||||
|
||||
def test_decay_spectrum_clip(decay_spectrum_chain):
|
||||
# Stable / non-emitting nuclides are removed unconditionally
|
||||
d = openmc.stats.DecaySpectrum(
|
||||
{'I135': 1e-3, 'Xe135': 5e-4, 'Xe136': 1.0, 'Cs135': 1.0},
|
||||
volume=10.0,
|
||||
)
|
||||
d_clip = d.clip()
|
||||
assert 'Xe136' not in d_clip.nuclides
|
||||
assert 'Cs135' not in d_clip.nuclides
|
||||
assert 'I135' in d_clip.nuclides
|
||||
assert 'Xe135' in d_clip.nuclides
|
||||
# Original is unchanged
|
||||
assert 'Xe136' in d.nuclides
|
||||
|
||||
# inplace=True modifies and returns the same object
|
||||
d_same = d.clip(inplace=True)
|
||||
assert d_same is d
|
||||
assert 'Xe136' not in d.nuclides
|
||||
|
||||
# A nuclide with negligible emission rate is removed by tolerance clipping.
|
||||
# U235 has a very small integral (~4e-17 Bq/atom) compared with I135 (~4e-5)
|
||||
d_tight = openmc.stats.DecaySpectrum(
|
||||
{'I135': 1e-3, 'U235': 1e-3}, volume=10.0
|
||||
)
|
||||
d_tight_clip = d_tight.clip(tolerance=1e-9)
|
||||
assert 'U235' not in d_tight_clip.nuclides
|
||||
assert 'I135' in d_tight_clip.nuclides
|
||||
|
||||
# All non-emitters -> empty nuclides dict
|
||||
d_empty = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0)
|
||||
d_empty.clip(inplace=True)
|
||||
assert d_empty.nuclides == {}
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue