mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Merge remote-tracking branch 'upstream/develop' into bootsrap_mgxs
This commit is contained in:
John Tramm
commit
986bb765d2
68 changed files with 3265 additions and 348 deletions
1
.gitignore
vendored
1
.gitignore
vendored
|
|
@ -25,6 +25,7 @@ examples/**/*.xml
|
|||
|
||||
# Documentation builds
|
||||
docs/build
|
||||
docs/doxygen/xml
|
||||
docs/source/_images/*.pdf
|
||||
docs/source/_images/*.aux
|
||||
docs/source/pythonapi/generated/
|
||||
|
|
|
|||
|
|
@ -7,6 +7,8 @@ build:
|
|||
jobs:
|
||||
post_checkout:
|
||||
- git fetch --unshallow || true
|
||||
- cd docs/doxygen && doxygen && cd -
|
||||
|
||||
sphinx:
|
||||
configuration: docs/source/conf.py
|
||||
|
||||
|
|
|
|||
|
|
@ -45,6 +45,7 @@ help:
|
|||
clean:
|
||||
-rm -rf $(BUILDDIR)/*
|
||||
-rm -rf source/pythonapi/generated/
|
||||
-rm -rf doxygen/xml
|
||||
|
||||
html:
|
||||
$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
|
||||
|
|
|
|||
13
docs/doxygen/Doxyfile
Normal file
13
docs/doxygen/Doxyfile
Normal file
|
|
@ -0,0 +1,13 @@
|
|||
# Doxyfile 1.9.1
|
||||
|
||||
# This file describes the settings to be used by the documentation system
|
||||
# doxygen (www.doxygen.org) for a project.
|
||||
|
||||
# Difference with default Doxyfile 1.9.1
|
||||
PROJECT_NAME = OpenMC
|
||||
QUIET = YES
|
||||
WARN_IF_UNDOCUMENTED = NO
|
||||
INPUT = ../../include/openmc/capi.h
|
||||
GENERATE_HTML = NO
|
||||
GENERATE_LATEX = NO
|
||||
GENERATE_XML = YES
|
||||
|
|
@ -46,43 +46,20 @@ Type Definitions
|
|||
Functions
|
||||
---------
|
||||
|
||||
.. c:function:: int openmc_calculate_volumes()
|
||||
..
|
||||
Once documentation is complete in capi.h, use:
|
||||
.. doxygenfile:: capi.h
|
||||
to populate this documentation without using
|
||||
.. doxygenfunction::
|
||||
for every function.
|
||||
|
||||
Run a stochastic volume calculation
|
||||
.. doxygenfunction:: openmc_calculate_volumes
|
||||
|
||||
:return: Return status (negative if an error occurred)
|
||||
:rtype: int
|
||||
.. doxygenfunction:: openmc_cell_get_fill
|
||||
|
||||
.. c:function:: int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
|
||||
.. doxygenfunction:: openmc_cell_get_id
|
||||
|
||||
Get the fill for a cell
|
||||
|
||||
:param int32_t index: Index in the cells array
|
||||
:param int* type: Type of the fill
|
||||
:param int32_t** indices: Array of material indices for cell
|
||||
:param int32_t* n: Length of indices array
|
||||
:return: Return status (negative if an error occurred)
|
||||
:rtype: int
|
||||
|
||||
.. c:function:: int openmc_cell_get_id(int32_t index, int32_t* id)
|
||||
|
||||
Get the ID of a cell
|
||||
|
||||
:param int32_t index: Index in the cells array
|
||||
:param int32_t* id: ID of the cell
|
||||
:return: Return status (negative if an error occurred)
|
||||
:rtype: int
|
||||
|
||||
.. c:function:: int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T)
|
||||
|
||||
Get the temperature of a cell
|
||||
|
||||
:param int32_t index: Index in the cells array
|
||||
:param int32_t* instance: Which instance of the cell. If a null pointer is passed, the temperature
|
||||
of the first instance is returned.
|
||||
:param double* T: temperature of the cell
|
||||
:return: Return status (negative if an error occurred)
|
||||
:rtype: int
|
||||
.. doxygenfunction:: openmc_cell_get_temperature
|
||||
|
||||
.. c:function:: int openmc_cell_get_density(int32_t index, const int32_t* instance, double* density)
|
||||
|
||||
|
|
|
|||
|
|
@ -11,7 +11,10 @@
|
|||
# All configuration values have a default; values that are commented out
|
||||
# serve to show the default.
|
||||
|
||||
import sys, os
|
||||
import os
|
||||
from pathlib import Path
|
||||
import subprocess
|
||||
import sys
|
||||
|
||||
# Determine if we're on Read the Docs server
|
||||
on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
|
||||
|
|
@ -37,6 +40,7 @@ sys.path.insert(0, os.path.abspath('../..'))
|
|||
# Add any Sphinx extension module names here, as strings. They can be extensions
|
||||
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
|
||||
extensions = [
|
||||
'breathe',
|
||||
'sphinx.ext.autodoc',
|
||||
'sphinx.ext.napoleon',
|
||||
'sphinx.ext.autosummary',
|
||||
|
|
@ -47,6 +51,8 @@ extensions = [
|
|||
]
|
||||
if not on_rtd:
|
||||
extensions.append('sphinxcontrib.rsvgconverter')
|
||||
doxygen_dir = Path(__file__).parents[1] / 'doxygen'
|
||||
subprocess.run(['doxygen'], cwd=doxygen_dir, check=True)
|
||||
|
||||
# Add any paths that contain templates here, relative to this directory.
|
||||
templates_path = ['_templates']
|
||||
|
|
@ -117,6 +123,11 @@ pygments_style = 'tango'
|
|||
# A list of ignored prefixes for module index sorting.
|
||||
#modindex_common_prefix = []
|
||||
|
||||
# -- Options breathe + doxygen -------------------------------------------------
|
||||
|
||||
breathe_projects = {"OpenMC": "../doxygen/xml"}
|
||||
breathe_default_project = "OpenMC"
|
||||
breathe_domain_by_file_pattern = {"*capi.h": "c"}
|
||||
|
||||
# -- Options for HTML output ---------------------------------------------------
|
||||
|
||||
|
|
|
|||
|
|
@ -111,7 +111,8 @@ The TC consists of the following individuals:
|
|||
- `Paul Romano <https://github.com/paulromano>`_
|
||||
- `Patrick Shriwise <https://github.com/pshriwise>`_
|
||||
- `Adam Nelson <https://github.com/nelsonag>`_
|
||||
- `Benoit Forget <https://github.com/bforget>`_
|
||||
- `Jonathan Shimwell <https://github.com/shimwell>`_
|
||||
- `John Tramm <https://github.com/jtramm>`_
|
||||
|
||||
The Project Lead is Paul Romano.
|
||||
|
||||
|
|
|
|||
|
|
@ -14,6 +14,11 @@ Python API. That is, from the root directory of the OpenMC repository:
|
|||
|
||||
python -m pip install ".[docs]"
|
||||
|
||||
The OpenMC documentation also uses Doxygen to automatically generate its
|
||||
C/C++ API documentation directly from the docstrings available in the source
|
||||
code. You will need to have a working installation of Doxygen to generate the
|
||||
documentation locally.
|
||||
|
||||
-----------------------------------
|
||||
Building Documentation as a Webpage
|
||||
-----------------------------------
|
||||
|
|
|
|||
|
|
@ -425,13 +425,13 @@ Each ``<dagmc_universe>`` element can have the following attributes or sub-eleme
|
|||
material) assignment. Required.
|
||||
|
||||
:temperature:
|
||||
Temperature(s) in [K] to assign to the cell. Must be ≥ 0. Multiple
|
||||
space-separated values may be given.
|
||||
Temperature(s) in [K] to assign to the cell. Must be greater than or equal
|
||||
to 0. Multiple space-separated values may be given.
|
||||
|
||||
*Default*: None
|
||||
|
||||
:density:
|
||||
Density in [g/cm³] to assign to the cell. Must be > 0. Requires a non-void
|
||||
Density in [g/cm³] to assign to the cell. Must be greater than 0. Requires a non-void
|
||||
material fill. Multiple space-separated values may be given.
|
||||
|
||||
*Default*: None
|
||||
|
|
|
|||
|
|
@ -777,9 +777,9 @@ attributes/sub-elements:
|
|||
*Default*: 1.0
|
||||
|
||||
:type:
|
||||
Indicator of source type. One of ``independent``, ``file``, ``compiled``, or
|
||||
``mesh``. The type of the source will be determined by this attribute if it
|
||||
is present.
|
||||
Indicator of source type. One of ``independent``, ``file``, ``compiled``,
|
||||
``mesh``, or ``tokamak``. The type of the source will be determined by this
|
||||
attribute if it is present.
|
||||
|
||||
:particle:
|
||||
The source particle type, specified as a PDG number or a string alias (e.g.,
|
||||
|
|
@ -1015,6 +1015,80 @@ attributes/sub-elements:
|
|||
mesh element and follows the format for :ref:`source_element`. The number of
|
||||
``<source>`` sub-elements should correspond to the number of mesh elements.
|
||||
|
||||
For a source with ``type="tokamak"``, the spatial distribution is described by
|
||||
a Miller-style flux-surface parameterization and the following sub-elements
|
||||
are used instead of the ``space`` element:
|
||||
|
||||
:major_radius:
|
||||
The major radius :math:`R_0` of the plasma in [cm].
|
||||
|
||||
:minor_radius:
|
||||
The minor radius :math:`a` of the plasma in [cm]. Must be smaller than
|
||||
``major_radius``.
|
||||
|
||||
:elongation:
|
||||
The plasma elongation :math:`\kappa` (must be > 0).
|
||||
|
||||
:triangularity:
|
||||
The plasma triangularity :math:`\delta` (must be in [-1, 1]). Negative
|
||||
values describe negative-triangularity plasmas.
|
||||
|
||||
:shafranov_shift:
|
||||
The Shafranov shift :math:`\Delta` in [cm] (must be >= 0 and less than
|
||||
``minor_radius``/2).
|
||||
|
||||
:r_over_a:
|
||||
A list of normalized minor-radius grid points :math:`r/a`. Must be strictly
|
||||
increasing, start at 0, and end at 1.
|
||||
|
||||
:emission_density:
|
||||
A list of neutron emission densities :math:`S(r)` evaluated at each
|
||||
``r_over_a`` grid point (arbitrary units, must be non-negative). Only the
|
||||
shape matters, since the profile is normalized internally. Values are
|
||||
interpolated linearly between grid points and the profile is refined on an
|
||||
internal grid for radial sampling. Must have the same length as
|
||||
``r_over_a`` and contain at least one positive value.
|
||||
|
||||
:phi_start:
|
||||
The starting toroidal angle in [rad].
|
||||
|
||||
*Default*: 0.0
|
||||
|
||||
:phi_extent:
|
||||
The toroidal angle extent in [rad]. The source is sampled uniformly in
|
||||
:math:`[\phi_\text{start},\ \phi_\text{start} + \phi_\text{extent}]`.
|
||||
|
||||
*Default*: :math:`2\pi`
|
||||
|
||||
:n_alpha:
|
||||
The number of poloidal-angle grid points used to build the sampling CDFs
|
||||
(must be > 2). Larger values reduce discretization bias; values below 51
|
||||
produce a warning.
|
||||
|
||||
*Default*: 101
|
||||
|
||||
:vertical_shift:
|
||||
A vertical shift of the plasma center in [cm].
|
||||
|
||||
*Default*: 0.0
|
||||
|
||||
:energy:
|
||||
For a tokamak source, one or more ``energy`` sub-elements specify the
|
||||
neutron energy distribution(s). Either a single distribution is given (used
|
||||
at all radii) or exactly one distribution per ``r_over_a`` grid point is
|
||||
given, in which case the energy is sampled from one of the two
|
||||
distributions bracketing the sampled radius, selected stochastically with
|
||||
probability proportional to the proximity of the radius to each grid point
|
||||
(stochastic interpolation). Each follows the format of a univariate
|
||||
probability distribution (see :ref:`univariate`).
|
||||
|
||||
:time:
|
||||
An optional ``time`` sub-element specifying the time distribution of source
|
||||
particles, following the format of a univariate probability distribution
|
||||
(see :ref:`univariate`).
|
||||
|
||||
*Default*: particles are born at :math:`t=0`
|
||||
|
||||
.. note:: Biased sampling can be applied to the spatial and energy distributions
|
||||
of a source by using the ``<bias>`` sub-element (see
|
||||
:ref:`univariate` for details on how to specify bias distributions).
|
||||
|
|
|
|||
|
|
@ -26,6 +26,7 @@ Simulation Settings
|
|||
openmc.FileSource
|
||||
openmc.CompiledSource
|
||||
openmc.MeshSource
|
||||
openmc.TokamakSource
|
||||
openmc.SourceParticle
|
||||
openmc.VolumeCalculation
|
||||
openmc.Settings
|
||||
|
|
|
|||
|
|
@ -132,8 +132,9 @@ can be run::
|
|||
r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes)
|
||||
|
||||
If not specified otherwise, a photon transport calculation is run at each time
|
||||
in the depletion schedule. That means in the case above, we would see three
|
||||
photon transport calculations. To specify specific times at which photon
|
||||
in the depletion schedule for which a decay photon source exists. Times without
|
||||
a decay photon source, such as the initial state of a model containing only
|
||||
stable nuclides, are omitted. To specify particular times at which photon
|
||||
transport calculations should be run, pass the ``photon_time_indices`` argument.
|
||||
For example, if we wanted to run a photon transport calculation only on the last
|
||||
time (after the 5 hour decay), we would run::
|
||||
|
|
@ -141,6 +142,19 @@ time (after the 5 hour decay), we would run::
|
|||
r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes,
|
||||
photon_time_indices=[2])
|
||||
|
||||
To attribute photon tally results to their parent radionuclides, set
|
||||
``by_parent_nuclide=True``. This automatically adds a
|
||||
:class:`openmc.ParentNuclideFilter` to every photon tally that does not already
|
||||
have one. The filter bins are the union of radionuclides contributing to the
|
||||
prepared decay photon sources. The resulting bins can be used directly when
|
||||
inspecting the tally results::
|
||||
|
||||
r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes,
|
||||
photon_time_indices=[2], by_parent_nuclide=True)
|
||||
|
||||
photon_tally = r2s.results['photon_tallies'][2][0]
|
||||
tally_by_parent = photon_tally.get_pandas_dataframe()
|
||||
|
||||
After an R2S calculation has been run, the :class:`~openmc.deplete.R2SManager`
|
||||
instance will have a ``results`` dictionary that allows you to directly access
|
||||
results from each of the steps. It will also write out all the output files into
|
||||
|
|
@ -148,12 +162,13 @@ a directory that is named "r2s_<timestamp>/". The ``output_dir`` argument to the
|
|||
:meth:`~openmc.deplete.R2SManager.run` method enables you to override the
|
||||
default output directory name if desired.
|
||||
|
||||
The :meth:`~openmc.deplete.R2SManager.run` method actually runs three
|
||||
The :meth:`~openmc.deplete.R2SManager.run` method actually runs four
|
||||
lower-level methods under the hood::
|
||||
|
||||
r2s.step1_neutron_transport(...)
|
||||
r2s.step2_activation(...)
|
||||
r2s.step3_photon_transport(...)
|
||||
r2s.step3_photon_source(...)
|
||||
r2s.step4_photon_transport(...)
|
||||
|
||||
For users looking for more control over the calculation, these lower-level
|
||||
methods can be used in lieu of the :meth:`openmc.deplete.R2SManager.run` method.
|
||||
|
|
@ -255,4 +270,3 @@ relevant tallies. This can be done with the aid of the
|
|||
|
||||
# Apply time correction factors
|
||||
tally = d1s.apply_time_correction(dose_tally, factors, time_index)
|
||||
|
||||
|
|
|
|||
|
|
@ -291,6 +291,53 @@ example, the following would generate a photon source::
|
|||
For a full list of all classes related to statistical distributions, see
|
||||
:ref:`pythonapi_stats`.
|
||||
|
||||
Tokamak Plasma Sources
|
||||
----------------------
|
||||
|
||||
For fusion applications, the :class:`openmc.TokamakSource` class provides a
|
||||
native parametric neutron source for tokamak plasmas. Rather than specifying
|
||||
spatial, angular, and energy distributions separately, the source is defined by
|
||||
the plasma geometry (using a `Miller-style flux-surface parameterization
|
||||
<https://doi.org/10.1063/1.872666>`_) and a radial emission profile. Source
|
||||
sites are sampled directly from the plasma volume without rejection.
|
||||
|
||||
The plasma shape is described by the major radius :math:`R_0`, minor radius
|
||||
:math:`a`, elongation :math:`\kappa`, triangularity :math:`\delta`, and
|
||||
Shafranov shift :math:`\Delta`. The neutron birth profile is given as an
|
||||
emission density :math:`S(r/a)` tabulated on a normalized minor-radius grid that
|
||||
runs from 0 (magnetic axis) to 1 (last closed flux surface); only the shape of
|
||||
the profile matters, since it is normalized internally. The emission density is
|
||||
linearly interpolated between the supplied points and refined internally for
|
||||
radial sampling. For example::
|
||||
|
||||
import numpy as np
|
||||
|
||||
r_over_a = np.linspace(0.0, 1.0, 50)
|
||||
emission = (1.0 - r_over_a**2)**2 # peaked on-axis profile
|
||||
|
||||
source = openmc.TokamakSource(
|
||||
major_radius=620.0, # cm
|
||||
minor_radius=200.0, # cm
|
||||
elongation=1.8,
|
||||
triangularity=0.45,
|
||||
shafranov_shift=10.0, # cm
|
||||
r_over_a=r_over_a,
|
||||
emission_density=emission,
|
||||
energy=openmc.stats.muir(e0=14.08e6, m_rat=5.0, kt=20000.0),
|
||||
)
|
||||
|
||||
settings.source = source
|
||||
|
||||
The ``energy`` argument accepts either a single
|
||||
:class:`~openmc.stats.Univariate` distribution applied at all radii, or a
|
||||
sequence with one distribution per ``r_over_a`` grid point to model a
|
||||
radially-varying neutron spectrum (energies are then sampled by stochastic
|
||||
interpolation between the two distributions bracketing the sampled radius). A
|
||||
time distribution can be given with the ``time`` argument; by default, particles
|
||||
are born at :math:`t=0`. The toroidal extent can be restricted with
|
||||
``phi_start`` and ``phi_extent`` to model a sector of the plasma, and
|
||||
``vertical_shift`` translates the plasma center along the z-axis.
|
||||
|
||||
File-based Sources
|
||||
------------------
|
||||
|
||||
|
|
|
|||
|
|
@ -9,14 +9,41 @@
|
|||
extern "C" {
|
||||
#endif
|
||||
|
||||
//! Run a stochastic volume calculation
|
||||
//
|
||||
//! \return Status (negative if an error occurred)
|
||||
int openmc_calculate_volumes();
|
||||
|
||||
int openmc_cell_filter_get_bins(
|
||||
int32_t index, const int32_t** cells, int32_t* n);
|
||||
|
||||
//! Get the fill for a cell
|
||||
//
|
||||
//! \param index Index in the cells array
|
||||
//! \param type Type of the fill
|
||||
//! \param indices Array of material indices for cell
|
||||
//! \param n Length of indices array
|
||||
//! \return Status (negative if an error occurred)
|
||||
int openmc_cell_get_fill(
|
||||
int32_t index, int* type, int32_t** indices, int32_t* n);
|
||||
|
||||
//! Get the ID of a cell
|
||||
//
|
||||
//! \param index Index in the cells array
|
||||
//! \param id ID of the cell
|
||||
//! \return Status (negative if an error occurred)
|
||||
int openmc_cell_get_id(int32_t index, int32_t* id);
|
||||
|
||||
//! Get the temperature of a cell
|
||||
//
|
||||
//! \param index Index in the cells array
|
||||
//! \param instance Which instance of the cell. If a null pointer is
|
||||
//! passed, the temperature of the first instance is returned.
|
||||
//! \param T temperature of the cell
|
||||
//!\return Status (negative if an error occurred)
|
||||
int openmc_cell_get_temperature(
|
||||
int32_t index, const int32_t* instance, double* T);
|
||||
|
||||
int openmc_cell_get_density(
|
||||
int32_t index, const int32_t* instance, double* rho);
|
||||
int openmc_cell_get_translation(int32_t index, double xyz[]);
|
||||
|
|
@ -280,7 +307,7 @@ int openmc_zernike_filter_set_params(
|
|||
int openmc_particle_filter_get_bins(int32_t idx, int32_t bins[]);
|
||||
|
||||
//! Sets the mesh and energy grid for CMFD reweight
|
||||
//! \param[in] meshtyally_id id of CMFD Mesh Tally
|
||||
//! \param[in] meshtally_id id of CMFD Mesh Tally
|
||||
//! \param[in] cmfd_indices indices storing spatial and energy dimensions of
|
||||
//! CMFD problem \param[in] norm CMFD normalization factor
|
||||
void openmc_initialize_mesh_egrid(
|
||||
|
|
|
|||
|
|
@ -64,14 +64,14 @@ void write_message(
|
|||
int level, const std::string& message, const Params&... fmt_args)
|
||||
{
|
||||
if (settings::verbosity >= level) {
|
||||
write_message(fmt::format(message, fmt_args...));
|
||||
write_message(fmt::format(fmt::runtime(message), fmt_args...));
|
||||
}
|
||||
}
|
||||
|
||||
template<typename... Params>
|
||||
void write_message(const std::string& message, const Params&... fmt_args)
|
||||
{
|
||||
write_message(fmt::format(message, fmt_args...));
|
||||
write_message(fmt::format(fmt::runtime(message), fmt_args...));
|
||||
}
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
|
|
|
|||
|
|
@ -3,9 +3,12 @@
|
|||
|
||||
#include <cmath>
|
||||
#include <cstdint>
|
||||
#include <unordered_map>
|
||||
#include <vector>
|
||||
|
||||
#include "openmc/array.h"
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/random_ray/source_region.h" // For hash_combine
|
||||
#include "openmc/vector.h"
|
||||
|
||||
namespace openmc {
|
||||
|
|
@ -13,6 +16,34 @@ namespace openmc {
|
|||
class BoundaryInfo;
|
||||
class GeometryState;
|
||||
|
||||
//==============================================================================
|
||||
//! OverlapKey to store cell and universe data of a single overlap, along with
|
||||
//! a functor for hashing an OverlapKey into an unordered_map.
|
||||
//==============================================================================
|
||||
|
||||
struct OverlapKey {
|
||||
int universe_id;
|
||||
int cell1_id;
|
||||
int cell2_id;
|
||||
|
||||
bool operator==(const OverlapKey& other) const
|
||||
{
|
||||
return universe_id == other.universe_id && cell1_id == other.cell1_id &&
|
||||
cell2_id == other.cell2_id;
|
||||
}
|
||||
};
|
||||
|
||||
struct OverlapKeyHash {
|
||||
std::size_t operator()(const OverlapKey& k) const
|
||||
{
|
||||
size_t seed = 0;
|
||||
hash_combine(seed, k.universe_id);
|
||||
hash_combine(seed, k.cell1_id);
|
||||
hash_combine(seed, k.cell2_id);
|
||||
return seed;
|
||||
}
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
// Global variables
|
||||
//==============================================================================
|
||||
|
|
@ -24,6 +55,10 @@ extern "C" int n_coord_levels; //!< Number of CSG coordinate levels
|
|||
|
||||
extern vector<int64_t> overlap_check_count;
|
||||
|
||||
// Overlap data structures get cleared every slice_data run
|
||||
extern vector<OverlapKey> overlap_keys;
|
||||
extern std::unordered_map<OverlapKey, int, OverlapKeyHash> overlap_key_index;
|
||||
|
||||
} // namespace model
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -38,8 +73,7 @@ inline bool coincident(double d1, double d2)
|
|||
//==============================================================================
|
||||
//! Check for overlapping cells at a particle's position.
|
||||
//==============================================================================
|
||||
|
||||
bool check_cell_overlap(GeometryState& p, bool error = true);
|
||||
int check_cell_overlap(GeometryState& p, bool error = true);
|
||||
|
||||
//==============================================================================
|
||||
//! Get the cell instance for a particle at the specified universe level
|
||||
|
|
|
|||
|
|
@ -213,6 +213,20 @@ double exprel(double x);
|
|||
//! \return log(1+x)/x without loss of precision near 0
|
||||
double log1prel(double x);
|
||||
|
||||
//! Evaluate the cylindrical Bessel function of the first kind J_n(x)
|
||||
//!
|
||||
//! Uses std::cyl_bessel_j where available (e.g., libstdc++). On standard
|
||||
//! library implementations lacking the C++17 special math functions (e.g.,
|
||||
//! libc++ on Apple Clang/LLVM), falls back to the ascending power series,
|
||||
//! which converges to machine precision for the small arguments (|x| <= 2)
|
||||
//! used in OpenMC. Unlike std::cyl_bessel_j, negative arguments are handled
|
||||
//! via the parity relation J_n(-x) = (-1)^n J_n(x).
|
||||
//!
|
||||
//! \param n Non-negative integer order of the Bessel function
|
||||
//! \param x Real argument
|
||||
//! \return J_n(x)
|
||||
double cyl_bessel_j(int n, double x);
|
||||
|
||||
//! Helper function to get index and interpolation function on an incident
|
||||
//! energy grid
|
||||
//!
|
||||
|
|
@ -233,5 +247,17 @@ void get_energy_index(
|
|||
|
||||
double standard_normal_cdf(double z);
|
||||
|
||||
//==============================================================================
|
||||
//! Return true if two floating-point values are approximately equal within a
|
||||
//! combined relative and absolute tolerance.
|
||||
//!
|
||||
//! \param a first floating point value
|
||||
//! \param b second floating point value
|
||||
//! \param rel_tol relative tolerance
|
||||
//! \param abs_tol absolute tolerance
|
||||
//! \return true if a and b are approximately equal, false otherwise
|
||||
//==============================================================================
|
||||
bool isclose(double a, double b, double rel_tol, double abs_tol);
|
||||
|
||||
} // namespace openmc
|
||||
#endif // OPENMC_MATH_FUNCTIONS_H
|
||||
|
|
|
|||
|
|
@ -582,10 +582,8 @@ public:
|
|||
// Methods and accessors
|
||||
|
||||
// Cross section caches
|
||||
NuclideMicroXS& neutron_xs(int i)
|
||||
{
|
||||
return neutron_xs_[i];
|
||||
} // Microscopic neutron cross sections
|
||||
// Microscopic neutron cross sections
|
||||
NuclideMicroXS& neutron_xs(int i) { return neutron_xs_[i]; }
|
||||
const NuclideMicroXS& neutron_xs(int i) const { return neutron_xs_[i]; }
|
||||
|
||||
// Microscopic photon cross sections
|
||||
|
|
|
|||
|
|
@ -5,6 +5,7 @@
|
|||
#include <sstream>
|
||||
#include <unordered_map>
|
||||
#include <unordered_set>
|
||||
#include <vector>
|
||||
|
||||
#include "openmc/tensor.h"
|
||||
#include "pugixml.hpp"
|
||||
|
|
@ -155,7 +156,7 @@ struct IdData {
|
|||
// Methods
|
||||
void set_value(size_t y, size_t x, const Particle& p, int level,
|
||||
Filter* filter = nullptr, FilterMatch* match = nullptr);
|
||||
void set_overlap(size_t y, size_t x);
|
||||
void set_overlap(size_t y, size_t x, int overlap_idx);
|
||||
|
||||
// Members
|
||||
tensor::Tensor<int32_t> data_; //!< 2D array of cell & material ids
|
||||
|
|
@ -168,7 +169,7 @@ struct PropertyData {
|
|||
// Methods
|
||||
void set_value(size_t y, size_t x, const Particle& p, int level,
|
||||
Filter* filter = nullptr, FilterMatch* match = nullptr);
|
||||
void set_overlap(size_t y, size_t x);
|
||||
void set_overlap(size_t y, size_t x, int overlap_idx);
|
||||
|
||||
// Members
|
||||
tensor::Tensor<double> data_; //!< 2D array of temperature & density data
|
||||
|
|
@ -181,7 +182,7 @@ struct RasterData {
|
|||
// Methods
|
||||
void set_value(size_t y, size_t x, const Particle& p, int level,
|
||||
Filter* filter = nullptr, FilterMatch* match = nullptr);
|
||||
void set_overlap(size_t y, size_t x);
|
||||
void set_overlap(size_t y, size_t x, int overlap_idx);
|
||||
|
||||
// Members
|
||||
tensor::Tensor<int32_t>
|
||||
|
|
@ -278,8 +279,11 @@ T SlicePlotBase::get_map(int32_t filter_index) const
|
|||
if (found_cell) {
|
||||
data.set_value(y, x, p, j, filter, &match);
|
||||
}
|
||||
if (show_overlaps_ && check_cell_overlap(p, false)) {
|
||||
data.set_overlap(y, x);
|
||||
if (show_overlaps_) {
|
||||
int overlap_idx = check_cell_overlap(p, false);
|
||||
if (overlap_idx >= 0) {
|
||||
data.set_overlap(y, x, overlap_idx);
|
||||
}
|
||||
}
|
||||
} // inner for
|
||||
}
|
||||
|
|
|
|||
|
|
@ -114,6 +114,32 @@ public:
|
|||
data_[size_++] = value;
|
||||
}
|
||||
|
||||
//! Increase the size of the container by count elements without assigning
|
||||
//! values to the new elements. Existing elements are preserved if the
|
||||
//! container needs to grow. This does not perform any thread safety checks.
|
||||
//
|
||||
//! \param count The number of elements to append
|
||||
//! \return The starting index of the appended range
|
||||
int64_t extend_uninitialized(int64_t count)
|
||||
{
|
||||
int64_t offset = size_;
|
||||
int64_t new_size = size_ + count;
|
||||
if (new_size > capacity_) {
|
||||
int64_t new_capacity = capacity_ == 0 ? 8 : capacity_;
|
||||
while (new_capacity < new_size) {
|
||||
new_capacity *= 2;
|
||||
}
|
||||
unique_ptr<T[]> new_data = make_unique<T[]>(new_capacity);
|
||||
if (size_ > 0) {
|
||||
std::copy_n(data_.get(), size_, new_data.get());
|
||||
}
|
||||
data_ = std::move(new_data);
|
||||
capacity_ = new_capacity;
|
||||
}
|
||||
size_ = new_size;
|
||||
return offset;
|
||||
}
|
||||
|
||||
//! Return the number of elements in the container
|
||||
int64_t size() { return size_; }
|
||||
int64_t size() const { return size_; }
|
||||
|
|
|
|||
|
|
@ -7,9 +7,12 @@
|
|||
#include <atomic>
|
||||
#include <limits>
|
||||
#include <unordered_set>
|
||||
#include <utility> // for pair
|
||||
|
||||
#include "pugixml.hpp"
|
||||
|
||||
#include "openmc/array.h"
|
||||
#include "openmc/distribution.h"
|
||||
#include "openmc/distribution_multi.h"
|
||||
#include "openmc/distribution_spatial.h"
|
||||
#include "openmc/memory.h"
|
||||
|
|
@ -260,6 +263,139 @@ private:
|
|||
vector<unique_ptr<IndependentSource>> sources_; //!< Source distributions
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
//! Parametric tokamak plasma neutron source
|
||||
//!
|
||||
//! This source samples neutron positions from a tokamak plasma geometry using
|
||||
//! Miller-style flux surface parameterization with user-specified emission
|
||||
//! profiles and energy distributions.
|
||||
//!
|
||||
//! Flux surface parameterization:
|
||||
//! R = R0 + r*cos(alpha + delta*sin(alpha)) + Delta*(1 - (r/a)^2)
|
||||
//! Z = kappa * r * sin(alpha)
|
||||
//!
|
||||
//! The sampling algorithm:
|
||||
//! 1. Sample minor radius r from precomputed CDF of S(r) * Jacobian
|
||||
//! 2. Sample poloidal angle alpha from conditional P(alpha|r) using mixture
|
||||
//! of precomputed CDFs weighted by functions of r
|
||||
//! 3. Sample energy and time from user-provided distribution(s)
|
||||
//! 4. Sample isotropic direction
|
||||
//! 5. Sample toroidal angle phi uniformly in [phi_start, phi_start +
|
||||
//! phi_extent]
|
||||
//! 6. Transform (r, alpha, phi) to Cartesian (x, y, z), applying the optional
|
||||
//! vertical shift of the plasma center
|
||||
//!
|
||||
//! The user provides the emission density S(r) directly (e.g., from transport
|
||||
//! codes like TRANSP, ASTRA, etc.), allowing full flexibility in reaction
|
||||
//! physics calculations. S(r) is a profile in arbitrary units sampled on the
|
||||
//! r_over_a grid; only its shape matters, since it is normalized internally.
|
||||
//! Energy distributions can be specified as either a single distribution for
|
||||
//! all r, or one distribution per radial point.
|
||||
//==============================================================================
|
||||
|
||||
class TokamakSource : public Source {
|
||||
public:
|
||||
// Constructors
|
||||
explicit TokamakSource(pugi::xml_node node);
|
||||
|
||||
//! Sample from the tokamak source distribution
|
||||
//! \param[inout] seed Pseudorandom seed pointer
|
||||
//! \return Sampled site
|
||||
SourceSite sample(uint64_t* seed) const override;
|
||||
|
||||
private:
|
||||
//==========================================================================
|
||||
// Private methods
|
||||
|
||||
//! Precompute data structures for efficient sampling
|
||||
void precompute_sampling_distributions();
|
||||
|
||||
//! Sample minor radius from marginal CDF
|
||||
//! \param seed Pseudorandom seed pointer
|
||||
//! \return Sampled r/a value
|
||||
double sample_r_over_a(uint64_t* seed) const;
|
||||
|
||||
//! Sample poloidal angle given r using mixture of precomputed CDFs
|
||||
//! \param r_norm Normalized minor radius r/a
|
||||
//! \param seed Pseudorandom seed pointer
|
||||
//! \return Sampled poloidal angle alpha [rad]
|
||||
double sample_poloidal_angle(double r_norm, uint64_t* seed) const;
|
||||
|
||||
//! Compute the k-th mixture weight w_k(r) * I_hat_k for poloidal sampling
|
||||
//! \param k Basis function index (0-5)
|
||||
//! \param r Normalized minor radius r/a
|
||||
//! \return Mixture weight for component k
|
||||
double mixture_weight(int k, double r) const;
|
||||
|
||||
//! Sample energy from the distribution(s)
|
||||
//! \param r_norm Normalized minor radius r/a (for distribution selection)
|
||||
//! \param seed Pseudorandom seed pointer
|
||||
//! \return (Sampled energy [eV], importance weight)
|
||||
std::pair<double, double> sample_energy(double r_norm, uint64_t* seed) const;
|
||||
|
||||
//! Transform from flux coordinates (r, alpha, phi) to Cartesian (x, y, z)
|
||||
//! \param r Minor radius [cm]
|
||||
//! \param alpha Poloidal angle [rad]
|
||||
//! \param phi Toroidal angle [rad]
|
||||
//! \return Position in Cartesian coordinates [cm]
|
||||
Position flux_to_cartesian(double r, double alpha, double phi) const;
|
||||
|
||||
//==========================================================================
|
||||
// Data members
|
||||
|
||||
// Emission profile (input)
|
||||
vector<double> r_over_a_; //!< Normalized minor radius grid points
|
||||
vector<double> emission_density_; //!< Emission density S(r) at grid points
|
||||
|
||||
// Energy distribution(s): either 1 for all r, or one per r point
|
||||
vector<unique_ptr<Distribution>> energy_dists_;
|
||||
|
||||
// Time distribution (defaults to a delta distribution at t=0)
|
||||
UPtrDist time_;
|
||||
|
||||
// Angular distribution (isotropic)
|
||||
UPtrAngle angle_;
|
||||
|
||||
// Tokamak geometry parameters
|
||||
double major_radius_; //!< Major radius R0 [cm]
|
||||
double minor_radius_; //!< Minor radius a [cm]
|
||||
double elongation_; //!< Elongation kappa
|
||||
double triangularity_; //!< Triangularity delta
|
||||
double shafranov_shift_; //!< Shafranov shift Delta [cm]
|
||||
double vertical_shift_; //!< Vertical shift of plasma center [cm]
|
||||
|
||||
// Normalized geometry parameters (precomputed for efficiency)
|
||||
double epsilon_; //!< Inverse aspect ratio a/R0
|
||||
double delta_tilde_; //!< Normalized Shafranov shift Delta/a
|
||||
|
||||
// Toroidal angle bounds
|
||||
double phi_start_; //!< Starting toroidal angle [rad]
|
||||
double phi_extent_; //!< Toroidal angle extent [rad]
|
||||
|
||||
// Precomputed distribution for radial sampling
|
||||
unique_ptr<Tabular> radial_dist_;
|
||||
|
||||
// Coefficients of the radial geometric polynomial: A*r - B*r^2 - C*r^3
|
||||
// Also used as the analytical normalization for poloidal mixture weights
|
||||
double radial_poly_a_; //!< 1 + ε·Δ̃
|
||||
double radial_poly_b_; //!< (3/8)·c₁·ε
|
||||
double radial_poly_c_; //!< 2·ε·Δ̃
|
||||
|
||||
// Precomputed Bernstein basis functions for poloidal angle sampling.
|
||||
// Using the factorization f(r_tilde, alpha) = R_tilde x J_tilde where:
|
||||
// R_tilde = b0*(1-r)^2 + 2*b1*r*(1-r) + b2*r^2 (Bernstein quadratic)
|
||||
// J_tilde = b3*(1-r) + b4*r (Bernstein linear)
|
||||
// The product gives 6 non-negative basis functions g_k(alpha) with
|
||||
// weights w_k(r_tilde) that are products of Bernstein polynomials.
|
||||
// Distributions are tabulated on [0, pi] exploiting up-down symmetry.
|
||||
static constexpr int N_POLOIDAL_BASIS = 6; //!< Number of basis functions
|
||||
int n_alpha_; //!< Number of poloidal angle grid points
|
||||
array<unique_ptr<Tabular>, N_POLOIDAL_BASIS>
|
||||
poloidal_dists_; //!< Distributions for each basis function g_k(alpha)
|
||||
array<double, N_POLOIDAL_BASIS>
|
||||
poloidal_integrals_; //!< Integrals of g_k(alpha) over [0, pi]
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
// Functions
|
||||
//==============================================================================
|
||||
|
|
|
|||
|
|
@ -52,8 +52,12 @@ struct WeightWindow {
|
|||
double weight_cutoff {DEFAULT_WEIGHT_CUTOFF};
|
||||
int max_split {10};
|
||||
|
||||
//! Whether the weight window is in a valid state
|
||||
bool is_valid() const { return lower_weight >= 0.0; }
|
||||
//! Whether the weight window is in a valid state. A non-positive lower
|
||||
//! bound indicates that no weight window information exists at this
|
||||
//! location (generators mark such cells with -1, and a lower bound of zero
|
||||
//! conventionally turns the weight window game off in a cell, as in MCNP
|
||||
//! wwinp files), in which case no weight window game is played.
|
||||
bool is_valid() const { return lower_weight > 0.0; }
|
||||
|
||||
//! Adjust the weight window by a constant factor
|
||||
void scale(double factor)
|
||||
|
|
|
|||
|
|
@ -1,14 +1,23 @@
|
|||
from __future__ import annotations
|
||||
|
||||
import itertools
|
||||
import json
|
||||
import os
|
||||
import re
|
||||
from pathlib import Path
|
||||
from math import sqrt, log
|
||||
from typing import TYPE_CHECKING, Literal
|
||||
from warnings import warn
|
||||
|
||||
from endf.data import (ATOMIC_NUMBER, ATOMIC_SYMBOL, ELEMENT_SYMBOL,
|
||||
EV_PER_MEV, K_BOLTZMANN, gnds_name, zam)
|
||||
|
||||
import openmc
|
||||
from openmc.checkvalue import PathLike
|
||||
|
||||
if TYPE_CHECKING:
|
||||
from openmc.deplete import Chain
|
||||
|
||||
gnds_name.__module__ = __name__
|
||||
zam.__module__ = __name__
|
||||
|
||||
|
|
@ -296,25 +305,47 @@ def atomic_weight(element):
|
|||
raise ValueError(f"No naturally-occurring isotopes for element '{element}'.")
|
||||
|
||||
|
||||
def half_life(isotope):
|
||||
def half_life(
|
||||
isotope: str,
|
||||
chain_file: Literal[False] | None | PathLike | Chain = False
|
||||
) -> float | None:
|
||||
"""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>`_.
|
||||
By default, half-life values are from the `ENDF/B-VIII.0 decay sublibrary
|
||||
<https://www.nndc.bnl.gov/endf-b8.0/download.html>`_. A depletion chain can
|
||||
also be used as the source of half-life values.
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
.. versionchanged:: 0.15.4
|
||||
Added the ``chain_file`` argument.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
isotope : str
|
||||
Name of isotope, e.g., 'Pu239'
|
||||
chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
|
||||
Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
|
||||
used. If ``None``, the chain specified by
|
||||
``openmc.config['chain_file']`` is used when available. If a path or
|
||||
:class:`openmc.deplete.Chain` is given, that chain is used. For ``None``
|
||||
or an explicit chain, nuclides absent from the chain fall back to
|
||||
ENDF/B-VIII.0 data.
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Half-life of isotope in [s]
|
||||
float or None
|
||||
Half-life of isotope in [s], or None if the isotope is stable
|
||||
|
||||
"""
|
||||
if chain_file is not False:
|
||||
if chain_file is not None or openmc.config.get('chain_file') is not None:
|
||||
# Local import avoids a circular dependency
|
||||
from openmc.deplete.chain import _get_chain
|
||||
chain = _get_chain(chain_file)
|
||||
if isotope in chain:
|
||||
return chain[isotope].half_life
|
||||
|
||||
global _HALF_LIFE
|
||||
if not _HALF_LIFE:
|
||||
# Load ENDF/B-VIII.0 data from JSON file
|
||||
|
|
@ -324,7 +355,10 @@ def half_life(isotope):
|
|||
return _HALF_LIFE.get(isotope.lower())
|
||||
|
||||
|
||||
def decay_constant(isotope):
|
||||
def decay_constant(
|
||||
isotope: str,
|
||||
chain_file: Literal[False] | None | PathLike | Chain = False
|
||||
) -> float:
|
||||
"""Return decay constant of isotope in [s^-1]
|
||||
|
||||
Decay constants are based on half-life values from the
|
||||
|
|
@ -333,10 +367,20 @@ def decay_constant(isotope):
|
|||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
.. versionchanged:: 0.15.4
|
||||
Added the ``chain_file`` argument.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
isotope : str
|
||||
Name of isotope, e.g., 'Pu239'
|
||||
chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
|
||||
Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
|
||||
used. If ``None``, the chain specified by
|
||||
``openmc.config['chain_file']`` is used when available. If a path or
|
||||
:class:`openmc.deplete.Chain` is given, that chain is used. For ``None``
|
||||
or an explicit chain, nuclides absent from the chain fall back to
|
||||
ENDF/B-VIII.0 data.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -348,7 +392,7 @@ def decay_constant(isotope):
|
|||
openmc.data.half_life
|
||||
|
||||
"""
|
||||
t = half_life(isotope)
|
||||
t = half_life(isotope, chain_file)
|
||||
return _LOG_TWO / t if t else 0.0
|
||||
|
||||
|
||||
|
|
@ -496,5 +540,3 @@ def isotopes(element: str) -> list[tuple[str, float]]:
|
|||
result.append(kv)
|
||||
|
||||
return result
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -4,50 +4,76 @@ import numpy as np
|
|||
|
||||
import openmc.checkvalue as cv
|
||||
|
||||
_FILES = {
|
||||
('icrp74', 'neutron'): Path('icrp74') / 'neutrons.txt',
|
||||
('icrp74', 'photon'): Path('icrp74') / 'photons.txt',
|
||||
('icrp116', 'electron'): Path('icrp116') / 'electrons.txt',
|
||||
('icrp116', 'helium'): Path('icrp116') / 'helium_ions.txt',
|
||||
('icrp116', 'mu-'): Path('icrp116') / 'negative_muons.txt',
|
||||
('icrp116', 'pi-'): Path('icrp116') / 'negative_pions.txt',
|
||||
('icrp116', 'neutron'): Path('icrp116') / 'neutrons.txt',
|
||||
('icrp116', 'photon'): Path('icrp116') / 'photons.txt',
|
||||
('icrp116', 'photon kerma'): Path('icrp116') / 'photons_kerma.txt',
|
||||
('icrp116', 'mu+'): Path('icrp116') / 'positive_muons.txt',
|
||||
('icrp116', 'pi+'): Path('icrp116') / 'positive_pions.txt',
|
||||
('icrp116', 'positron'): Path('icrp116') / 'positrons.txt',
|
||||
('icrp116', 'proton'): Path('icrp116') / 'protons.txt',
|
||||
_FULL_GEOMETRIES = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO')
|
||||
_LIMITED_GEOMETRIES = ('AP', 'PA', 'ISO')
|
||||
|
||||
_TABLES = {
|
||||
('icrp74', 'effective', 'neutron'): (
|
||||
Path('icrp74') / 'neutrons.txt', _FULL_GEOMETRIES),
|
||||
('icrp74', 'effective', 'photon'): (
|
||||
Path('icrp74') / 'photons.txt', _FULL_GEOMETRIES),
|
||||
('icrp116', 'effective', 'electron'): (
|
||||
Path('icrp116') / 'electrons.txt', _LIMITED_GEOMETRIES),
|
||||
('icrp116', 'effective', 'helium'): (
|
||||
Path('icrp116') / 'helium_ions.txt', _LIMITED_GEOMETRIES),
|
||||
('icrp116', 'effective', 'mu-'): (
|
||||
Path('icrp116') / 'negative_muons.txt', _LIMITED_GEOMETRIES),
|
||||
('icrp116', 'effective', 'pi-'): (
|
||||
Path('icrp116') / 'negative_pions.txt', _LIMITED_GEOMETRIES),
|
||||
('icrp116', 'effective', 'neutron'): (
|
||||
Path('icrp116') / 'neutrons.txt', _FULL_GEOMETRIES),
|
||||
('icrp116', 'effective', 'photon'): (
|
||||
Path('icrp116') / 'photons.txt', _FULL_GEOMETRIES),
|
||||
('icrp116', 'effective', 'photon kerma'): (
|
||||
Path('icrp116') / 'photons_kerma.txt', _FULL_GEOMETRIES),
|
||||
('icrp116', 'effective', 'mu+'): (
|
||||
Path('icrp116') / 'positive_muons.txt', _LIMITED_GEOMETRIES),
|
||||
('icrp116', 'effective', 'pi+'): (
|
||||
Path('icrp116') / 'positive_pions.txt', _LIMITED_GEOMETRIES),
|
||||
('icrp116', 'effective', 'positron'): (
|
||||
Path('icrp116') / 'positrons.txt', _LIMITED_GEOMETRIES),
|
||||
('icrp116', 'effective', 'proton'): (
|
||||
Path('icrp116') / 'protons.txt', _FULL_GEOMETRIES),
|
||||
('icrp74', 'ambient', 'neutron'): (
|
||||
Path('icrp74') / 'neutrons_H10.txt', None),
|
||||
('icrp74', 'ambient', 'photon'): (
|
||||
Path('icrp74') / 'photons_H10.txt', None),
|
||||
}
|
||||
|
||||
_DOSE_TABLES = {}
|
||||
|
||||
|
||||
def _load_dose_icrp(data_source: str, particle: str):
|
||||
"""Load effective dose tables from text files.
|
||||
def _load_dose_table(data_source: str, dose_quantity: str, particle: str):
|
||||
"""Load dose tables from text files.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
data_source : {'icrp74', 'icrp116'}
|
||||
The dose conversion data source to use
|
||||
dose_quantity : {'effective', 'ambient'}
|
||||
Dose quantity to load. 'ambient' corresponds to ambient dose
|
||||
equivalent H*(10).
|
||||
particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'}
|
||||
Incident particle
|
||||
|
||||
"""
|
||||
path = Path(__file__).parent / _FILES[data_source, particle]
|
||||
key = (data_source, dose_quantity, particle)
|
||||
path = Path(__file__).parent / _TABLES[key][0]
|
||||
data = np.loadtxt(path, skiprows=3, encoding='utf-8')
|
||||
data[:, 0] *= 1e6 # Change energies to eV
|
||||
_DOSE_TABLES[data_source, particle] = data
|
||||
_DOSE_TABLES[key] = data
|
||||
|
||||
|
||||
def dose_coefficients(particle, geometry='AP', data_source='icrp116'):
|
||||
"""Return effective dose conversion coefficients.
|
||||
def dose_coefficients(
|
||||
particle, geometry='AP', data_source='icrp116', dose_quantity='effective'
|
||||
):
|
||||
"""Return dose conversion coefficients.
|
||||
|
||||
This function provides fluence (and air kerma) to effective or ambient dose
|
||||
(H*(10)) conversion coefficients for various types of external exposures
|
||||
based on values in ICRP publications. Corrected values found in a
|
||||
corrigendum are used rather than the values in the original report.
|
||||
Available libraries include `ICRP Publication 74
|
||||
This function provides fluence (and air kerma) to effective dose or ambient
|
||||
dose equivalent (H*(10)) conversion coefficients for various types of
|
||||
external exposures based on values in ICRP publications. Corrected values
|
||||
found in a corrigendum are used rather than the values in the original
|
||||
report. Available libraries include `ICRP Publication 74
|
||||
<https://doi.org/10.1016/S0146-6453(96)90010-X>` and `ICRP Publication 116
|
||||
<https://doi.org/10.1016/j.icrp.2011.10.001>`.
|
||||
|
||||
|
|
@ -63,45 +89,58 @@ def dose_coefficients(particle, geometry='AP', data_source='icrp116'):
|
|||
particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'}
|
||||
Incident particle
|
||||
geometry : {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'}
|
||||
Irradiation geometry assumed. Refer to ICRP-116 (Section 3.2) for the
|
||||
meaning of the options here.
|
||||
Irradiation geometry assumed for effective dose coefficients. Refer to
|
||||
ICRP-116 (Section 3.2) for the meaning of the options here. This
|
||||
argument does not apply when ``dose_quantity`` is 'ambient'.
|
||||
data_source : {'icrp74', 'icrp116'}
|
||||
The data source for the effective dose conversion coefficients.
|
||||
The data source for the dose conversion coefficients.
|
||||
dose_quantity : {'effective', 'ambient'}
|
||||
Dose quantity to return. 'effective' returns effective dose
|
||||
coefficients; 'ambient' returns ambient dose equivalent (H*(10))
|
||||
coefficients.
|
||||
|
||||
Returns
|
||||
-------
|
||||
energy : numpy.ndarray
|
||||
Energies at which dose conversion coefficients are given
|
||||
dose_coeffs : numpy.ndarray
|
||||
Effective dose coefficients in [pSv cm^2] at provided energies. For
|
||||
'photon kerma', the coefficients are given in [Sv/Gy].
|
||||
Dose coefficients in [pSv cm^2] at provided energies. For 'photon
|
||||
kerma', the coefficients are given in [Sv/Gy].
|
||||
|
||||
"""
|
||||
|
||||
cv.check_value('geometry', geometry, {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'})
|
||||
cv.check_value('geometry', geometry, _FULL_GEOMETRIES)
|
||||
cv.check_value('data_source', data_source, {'icrp74', 'icrp116'})
|
||||
cv.check_value('dose_quantity', dose_quantity, {'effective', 'ambient'})
|
||||
|
||||
if (data_source, particle) not in _FILES:
|
||||
available_particles = sorted({p for (ds, p) in _FILES if ds == data_source})
|
||||
key = (data_source, dose_quantity, particle)
|
||||
if key not in _TABLES:
|
||||
available_particles = sorted(
|
||||
p for ds, dq, p in _TABLES
|
||||
if ds == data_source and dq == dose_quantity
|
||||
)
|
||||
msg = (
|
||||
f"'{particle}' has no dose data in data source {data_source}. "
|
||||
f"Available particles for {data_source} are: {available_particles}"
|
||||
f"'{particle}' has no {dose_quantity} dose data in data source "
|
||||
f"{data_source}. Available particles for {data_source} "
|
||||
f"with dose quantity {dose_quantity} are: {available_particles}"
|
||||
)
|
||||
raise ValueError(msg)
|
||||
elif (data_source, particle) not in _DOSE_TABLES:
|
||||
_load_dose_icrp(data_source, particle)
|
||||
elif key not in _DOSE_TABLES:
|
||||
_load_dose_table(data_source, dose_quantity, particle)
|
||||
|
||||
# Get all data for selected particle
|
||||
data = _DOSE_TABLES[data_source, particle]
|
||||
data = _DOSE_TABLES[key]
|
||||
columns = _TABLES[key][1]
|
||||
|
||||
# Determine index for selected geometry
|
||||
if particle in ('neutron', 'photon', 'proton', 'photon kerma'):
|
||||
columns = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO')
|
||||
if columns is None:
|
||||
if geometry != 'AP':
|
||||
raise ValueError(
|
||||
"Irradiation geometry is not defined for ambient dose "
|
||||
"equivalent coefficients. Use the default geometry='AP'."
|
||||
)
|
||||
index = 0
|
||||
else:
|
||||
columns = ('AP', 'PA', 'ISO')
|
||||
index = columns.index(geometry)
|
||||
index = columns.index(geometry)
|
||||
|
||||
# Pull out energy and dose from table
|
||||
energy = data[:, 0].copy()
|
||||
dose_coeffs = data[:, index + 1].copy()
|
||||
return energy, dose_coeffs
|
||||
|
|
|
|||
50
openmc/data/dose/icrp74/neutrons_H10.txt
Normal file
50
openmc/data/dose/icrp74/neutrons_H10.txt
Normal file
|
|
@ -0,0 +1,50 @@
|
|||
Neutrons: Ambient per fluence, in units of pSv cm², for monoenergetic particles incident.
|
||||
|
||||
Energy (MeV) Dose
|
||||
1.00E-09 6.60
|
||||
1.00E-08 9.00
|
||||
2.53E-08 10.6
|
||||
1.00E-07 12.9
|
||||
2.00E-07 13.5
|
||||
5.00E-07 13.6
|
||||
1.00E-06 13.3
|
||||
2.00E-06 12.9
|
||||
5.00E-06 12.0
|
||||
1.00E-05 11.3
|
||||
2.00E-05 10.6
|
||||
5.00E-05 9.90
|
||||
1.00E-04 9.40
|
||||
2.00E-04 8.90
|
||||
5.00E-04 8.30
|
||||
1.00E-03 7.90
|
||||
2.00E-03 7.70
|
||||
5.00E-03 8.00
|
||||
1.00E-02 10.5
|
||||
2.00E-02 16.6
|
||||
3.00E-02 23.7
|
||||
5.00E-02 41.1
|
||||
7.00E-02 60.0
|
||||
1.00E-01 88.0
|
||||
1.50E-01 132
|
||||
2.00E-01 170
|
||||
3.00E-01 233
|
||||
5.00E-01 322
|
||||
7.00E-01 375
|
||||
9.00E-01 400
|
||||
1 416
|
||||
1.2 425
|
||||
2 420
|
||||
3 412
|
||||
4 408
|
||||
5 405
|
||||
6 400
|
||||
7 405
|
||||
8 409
|
||||
9 420
|
||||
10 440
|
||||
12 480
|
||||
14 520
|
||||
15 540
|
||||
16 555
|
||||
18 570
|
||||
20 600
|
||||
28
openmc/data/dose/icrp74/photons_H10.txt
Normal file
28
openmc/data/dose/icrp74/photons_H10.txt
Normal file
|
|
@ -0,0 +1,28 @@
|
|||
Photons: Ambient dose (H*10) per fluence, in units of pSv cm²
|
||||
|
||||
Energy (MeV) Dose
|
||||
0.010 0.061
|
||||
0.015 0.83
|
||||
0.020 1.05
|
||||
0.030 0.81
|
||||
0.040 0.64
|
||||
0.050 0.55
|
||||
0.060 0.51
|
||||
0.080 0.53
|
||||
0.100 0.61
|
||||
0.150 0.89
|
||||
0.200 1.20
|
||||
0.300 1.80
|
||||
0.400 2.38
|
||||
0.500 2.93
|
||||
0.600 3.44
|
||||
0.800 4.38
|
||||
1 5.20
|
||||
1.5 6.90
|
||||
2 8.60
|
||||
3 11.1
|
||||
4 13.4
|
||||
5 15.5
|
||||
6 17.6
|
||||
8 21.6
|
||||
10 25.6
|
||||
|
|
@ -233,7 +233,8 @@ def get_microxs_and_flux(
|
|||
model.export_to_model_xml()
|
||||
comm.barrier()
|
||||
# Reinitialize with tallies
|
||||
openmc.lib.init(intracomm=comm)
|
||||
output = run_kwargs.get('output', True) if run_kwargs else True
|
||||
openmc.lib.init(intracomm=comm, output=output)
|
||||
|
||||
with TemporaryDirectory() as temp_dir:
|
||||
# Indicate to run in temporary directory unless being executed through
|
||||
|
|
|
|||
|
|
@ -4,6 +4,7 @@ from contextlib import nullcontext
|
|||
import copy
|
||||
from datetime import datetime
|
||||
import json
|
||||
from numbers import Integral
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
|
|
@ -156,6 +157,7 @@ class R2SManager:
|
|||
mat_vol_kwargs: dict | None = None,
|
||||
run_kwargs: dict | None = None,
|
||||
operator_kwargs: dict | None = None,
|
||||
by_parent_nuclide: bool = False,
|
||||
):
|
||||
"""Run the R2S calculation.
|
||||
|
||||
|
|
@ -179,7 +181,7 @@ class R2SManager:
|
|||
timesteps. For example, if two timesteps are specified, the array of
|
||||
times would contain three entries, and [2] would indicate computing
|
||||
photon results at the last time. A value of None indicates to run
|
||||
photon transport for each time.
|
||||
photon transport at each time that has a decay photon source.
|
||||
output_dir : PathLike, optional
|
||||
Path to directory where R2S calculation outputs will be saved. If
|
||||
not provided, a timestamped directory 'r2s_YYYY-MM-DDTHH-MM-SS' is
|
||||
|
|
@ -207,6 +209,11 @@ class R2SManager:
|
|||
operator_kwargs : dict, optional
|
||||
Additional keyword arguments passed to
|
||||
:class:`openmc.deplete.IndependentOperator`.
|
||||
by_parent_nuclide : bool, optional
|
||||
Whether to score photon tallies separately for each parent
|
||||
radionuclide. A :class:`~openmc.ParentNuclideFilter` is added to
|
||||
tallies that do not already contain one, with bins determined from
|
||||
the prepared decay photon sources.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -256,9 +263,13 @@ class R2SManager:
|
|||
timesteps, source_rates, timestep_units,
|
||||
output_dir / 'activation', operator_kwargs=operator_kwargs
|
||||
)
|
||||
self.step3_photon_transport(
|
||||
self.step3_photon_source(
|
||||
photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
|
||||
mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
|
||||
mat_vol_kwargs=mat_vol_kwargs,
|
||||
)
|
||||
self.step4_photon_transport(
|
||||
output_dir / 'photon_transport', run_kwargs=run_kwargs,
|
||||
by_parent_nuclide=by_parent_nuclide,
|
||||
)
|
||||
|
||||
return output_dir
|
||||
|
|
@ -347,7 +358,7 @@ class R2SManager:
|
|||
|
||||
# Run neutron transport and get fluxes and micros. Run via openmc.lib to
|
||||
# maintain a consistent parallelism strategy with the activation step.
|
||||
with TemporarySession():
|
||||
with TemporarySession(output=False):
|
||||
self.results['fluxes'], self.results['micros'] = get_microxs_and_flux(
|
||||
self.neutron_model, domains, **micro_kwargs)
|
||||
|
||||
|
|
@ -438,23 +449,20 @@ class R2SManager:
|
|||
# Get depletion results
|
||||
self.results['depletion_results'] = Results(output_path)
|
||||
|
||||
def step3_photon_transport(
|
||||
def step3_photon_source(
|
||||
self,
|
||||
time_indices: Sequence[int] | None = None,
|
||||
bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
|
||||
output_dir: PathLike = 'photon_transport',
|
||||
mat_vol_kwargs: dict | None = None,
|
||||
run_kwargs: dict | None = None,
|
||||
):
|
||||
"""Run the photon transport step.
|
||||
"""Create decay photon sources.
|
||||
|
||||
This step performs photon transport calculations using decay photon
|
||||
sources created from the activated materials. For each specified time,
|
||||
it creates appropriate photon sources and runs a transport calculation.
|
||||
In mesh-based mode, the sources are created using the mesh material
|
||||
volumes, while in cell-based mode, they are created using bounding boxes
|
||||
for each cell. This step will populate the 'photon_tallies' key in the
|
||||
results dictionary.
|
||||
This step creates decay photon sources from the activated materials for
|
||||
each specified time. In mesh-based mode, the sources are created using
|
||||
mesh material volumes, while in cell-based mode, they are created using
|
||||
bounding boxes for each cell. This step will populate the
|
||||
'photon_sources' key in the results dictionary.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
|
|
@ -464,40 +472,54 @@ class R2SManager:
|
|||
timesteps. For example, if two timesteps are specified, the array of
|
||||
times would contain three entries, and [2] would indicate computing
|
||||
photon results at the last time. A value of None indicates to run
|
||||
photon transport for each time.
|
||||
photon transport at each time that has a decay photon source.
|
||||
bounding_boxes : dict[int, openmc.BoundingBox], optional
|
||||
Dictionary mapping cell IDs to bounding boxes used for spatial
|
||||
source sampling in cell-based R2S calculations. Required if method
|
||||
is 'cell-based'.
|
||||
output_dir : PathLike, optional
|
||||
Path to directory where photon transport outputs will be saved.
|
||||
Path to directory where photon source outputs will be saved.
|
||||
mat_vol_kwargs : dict, optional
|
||||
Additional keyword arguments passed to
|
||||
:meth:`openmc.MeshBase.material_volumes`.
|
||||
run_kwargs : dict, optional
|
||||
Additional keyword arguments passed to :meth:`openmc.Model.run`
|
||||
during the photon transport step. By default, output is disabled.
|
||||
"""
|
||||
|
||||
# Do not retain sources from an earlier successful call if this source
|
||||
# preparation attempt fails.
|
||||
self.results.pop('photon_sources', None)
|
||||
|
||||
# TODO: Automatically determine bounding box for each cell
|
||||
if bounding_boxes is None and self.method == 'cell-based':
|
||||
raise ValueError("bounding_boxes must be provided for cell-based "
|
||||
"R2S calculations.")
|
||||
|
||||
# Set default run arguments if not provided
|
||||
if run_kwargs is None:
|
||||
run_kwargs = {}
|
||||
run_kwargs.setdefault('output', False)
|
||||
|
||||
# Write out JSON file with tally IDs that can be used for loading
|
||||
# results
|
||||
output_dir = Path(output_dir)
|
||||
output_dir.mkdir(parents=True, exist_ok=True)
|
||||
|
||||
# Get default time indices if not provided
|
||||
# Determine and validate time indices before preparing source data.
|
||||
n_steps = len(self.results['depletion_results'])
|
||||
implicit_time_indices = time_indices is None
|
||||
if time_indices is None:
|
||||
n_steps = len(self.results['depletion_results'])
|
||||
time_indices = list(range(n_steps))
|
||||
else:
|
||||
time_indices = list(time_indices)
|
||||
if not time_indices:
|
||||
raise ValueError('time_indices must contain at least one index')
|
||||
|
||||
normalized_indices = []
|
||||
for index in time_indices:
|
||||
if isinstance(index, bool) or not isinstance(index, Integral):
|
||||
raise TypeError('time_indices must contain only integers')
|
||||
index = int(index)
|
||||
if index < -n_steps or index >= n_steps:
|
||||
raise IndexError(
|
||||
f'Photon time index {index} is out of range for '
|
||||
f'{n_steps} depletion results')
|
||||
normalized_index = index % n_steps
|
||||
normalized_indices.append(normalized_index)
|
||||
|
||||
# Remove duplicates while preserving order
|
||||
time_indices = list(dict.fromkeys(normalized_indices))
|
||||
|
||||
# Check whether the photon model is different
|
||||
neutron_univ = self.neutron_model.geometry.root_universe
|
||||
|
|
@ -523,13 +545,6 @@ class R2SManager:
|
|||
|
||||
self.results['mesh_material_volumes_photon'] = photon_mmv_list
|
||||
|
||||
if comm.rank == 0:
|
||||
tally_ids = [tally.id for tally in self.photon_model.tallies]
|
||||
with open(output_dir / 'tally_ids.json', 'w') as f:
|
||||
json.dump(tally_ids, f)
|
||||
|
||||
self.results['photon_tallies'] = {}
|
||||
|
||||
# Get dictionary of cells in the photon model
|
||||
if different_photon_model:
|
||||
photon_cells = self.photon_model.geometry.get_all_cells()
|
||||
|
|
@ -547,16 +562,100 @@ class R2SManager:
|
|||
continue
|
||||
work_items.append((cell, original_mat, bounding_boxes[cell.id]))
|
||||
|
||||
# Ensure photon transport is enabled in settings
|
||||
# Create decay photon sources for each time index
|
||||
photon_sources = {
|
||||
time_index: self._create_photon_sources(time_index, work_items)
|
||||
for time_index in time_indices
|
||||
}
|
||||
|
||||
# Determine if any times have no decay photon sources. If the user
|
||||
# didn't specify any specific time indices, remove those times from the
|
||||
# photon_sources dictionary. If the user did specify time indices, raise
|
||||
# an error if any of those times have no decay photon sources.
|
||||
empty_indices = [
|
||||
time_index for time_index, sources in photon_sources.items()
|
||||
if not sources
|
||||
]
|
||||
if implicit_time_indices:
|
||||
for time_index in empty_indices:
|
||||
del photon_sources[time_index]
|
||||
if not photon_sources:
|
||||
raise RuntimeError(
|
||||
'No decay photon sources were found at any depletion time')
|
||||
elif empty_indices:
|
||||
indices = ', '.join(str(index) for index in empty_indices)
|
||||
raise RuntimeError(
|
||||
f'No decay photon source was found for requested time '
|
||||
f'indices: {indices}')
|
||||
|
||||
self.results['photon_sources'] = photon_sources
|
||||
|
||||
def step4_photon_transport(
|
||||
self,
|
||||
output_dir: PathLike = 'photon_transport',
|
||||
run_kwargs: dict | None = None,
|
||||
by_parent_nuclide: bool = False,
|
||||
):
|
||||
"""Run photon transport using prepared decay photon sources.
|
||||
|
||||
This step runs a photon transport calculation for each source list
|
||||
created by :meth:`step3_photon_source`. It will populate the
|
||||
'photon_tallies' key in the results dictionary.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
output_dir : PathLike, optional
|
||||
Path to directory where photon transport outputs will be saved.
|
||||
run_kwargs : dict, optional
|
||||
Additional keyword arguments passed to :meth:`openmc.Model.run`.
|
||||
By default, output is disabled.
|
||||
by_parent_nuclide : bool, optional
|
||||
Whether to score photon tallies separately for each parent
|
||||
radionuclide. A :class:`~openmc.ParentNuclideFilter` is added to
|
||||
tallies that do not already contain one, with bins determined from
|
||||
the prepared decay photon sources.
|
||||
"""
|
||||
if 'photon_sources' not in self.results:
|
||||
raise RuntimeError(
|
||||
'Photon sources must be created with step3_photon_source '
|
||||
'before running photon transport.')
|
||||
photon_sources = self.results['photon_sources']
|
||||
if not photon_sources:
|
||||
raise RuntimeError(
|
||||
'No decay photon sources are available for transport')
|
||||
|
||||
if by_parent_nuclide:
|
||||
radionuclides = sorted({
|
||||
nuclide
|
||||
for sources in photon_sources.values()
|
||||
for source in sources
|
||||
for nuclide in source.energy.nuclides
|
||||
})
|
||||
|
||||
if radionuclides:
|
||||
parent_filter = openmc.ParentNuclideFilter(radionuclides)
|
||||
for tally in self.photon_model.tallies:
|
||||
if not tally.contains_filter(openmc.ParentNuclideFilter):
|
||||
tally.filters.append(parent_filter)
|
||||
|
||||
if run_kwargs is None:
|
||||
run_kwargs = {}
|
||||
run_kwargs.setdefault('output', False)
|
||||
|
||||
output_dir = Path(output_dir)
|
||||
output_dir.mkdir(parents=True, exist_ok=True)
|
||||
|
||||
if comm.rank == 0:
|
||||
tally_ids = [tally.id for tally in self.photon_model.tallies]
|
||||
with open(output_dir / 'tally_ids.json', 'w') as f:
|
||||
json.dump(tally_ids, f)
|
||||
|
||||
self.results['photon_tallies'] = {}
|
||||
|
||||
# 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'])
|
||||
|
||||
# Build decay photon sources and assign to the photon model
|
||||
sources = self._create_photon_sources(time_index, work_items)
|
||||
for time_index, sources in photon_sources.items():
|
||||
self.photon_model.settings.source = sources
|
||||
|
||||
# Run photon transport calculation
|
||||
|
|
|
|||
|
|
@ -1,12 +1,17 @@
|
|||
from __future__ import annotations
|
||||
|
||||
import numbers
|
||||
import bisect
|
||||
import math
|
||||
from collections.abc import Iterable
|
||||
from typing import Literal
|
||||
from warnings import warn
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
|
||||
import openmc
|
||||
from .chain import Chain, _get_chain
|
||||
from .stepresult import StepResult, VERSION_RESULTS
|
||||
import openmc.checkvalue as cv
|
||||
from openmc.data import atomic_mass, AVOGADRO
|
||||
|
|
@ -103,7 +108,8 @@ class Results(list):
|
|||
mat: Material | str,
|
||||
units: str = "Bq/cm3",
|
||||
by_nuclide: bool = False,
|
||||
volume: float | None = None
|
||||
volume: float | None = None,
|
||||
chain_file: Literal[False] | None | PathLike | Chain = None
|
||||
) -> tuple[np.ndarray, np.ndarray | list[dict]]:
|
||||
"""Get activity of material over time.
|
||||
|
||||
|
|
@ -115,22 +121,31 @@ class Results(list):
|
|||
Material object or material id to evaluate
|
||||
units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3'}
|
||||
Specifies the type of activity to return, options include total
|
||||
activity [Bq], specific [Bq/g, Bq/kg] or volumetric activity [Bq/cm3].
|
||||
activity [Bq], specific [Bq/g, Bq/kg] or volumetric activity
|
||||
[Bq/cm3].
|
||||
by_nuclide : bool
|
||||
Specifies if the activity should be returned for the material as a
|
||||
whole or per nuclide. Default is False.
|
||||
volume : float, optional
|
||||
Volume of the material. If not passed, defaults to using the
|
||||
:attr:`Material.volume` attribute.
|
||||
chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
|
||||
Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
|
||||
used. If ``None``, the chain specified by
|
||||
``openmc.config['chain_file']`` is used when available. If a path or
|
||||
:class:`openmc.deplete.Chain` is given, that chain is used. For
|
||||
``None`` or an explicit chain, nuclides absent from the chain fall
|
||||
back to ENDF/B-VIII.0 data.
|
||||
|
||||
.. versionadded:: 0.15.4
|
||||
|
||||
Returns
|
||||
-------
|
||||
times : numpy.ndarray
|
||||
Array of times in [s]
|
||||
activities : numpy.ndarray or List[dict]
|
||||
Array of total activities if by_nuclide = False (default)
|
||||
or list of dictionaries of activities by nuclide if
|
||||
by_nuclide = True.
|
||||
Array of total activities if by_nuclide = False (default) or list of
|
||||
dictionaries of activities by nuclide if by_nuclide = True.
|
||||
|
||||
"""
|
||||
if isinstance(mat, Material):
|
||||
|
|
@ -140,6 +155,13 @@ class Results(list):
|
|||
else:
|
||||
raise TypeError('mat should be of type openmc.Material or str')
|
||||
|
||||
if chain_file is not False:
|
||||
if chain_file is None:
|
||||
if openmc.config.get('chain_file') is not None:
|
||||
chain_file = _get_chain(None)
|
||||
else:
|
||||
chain_file = _get_chain(chain_file)
|
||||
|
||||
times = np.empty_like(self, dtype=float)
|
||||
if by_nuclide:
|
||||
activities = [None] * len(self)
|
||||
|
|
@ -149,7 +171,8 @@ class Results(list):
|
|||
# Evaluate activity for each depletion time
|
||||
for i, result in enumerate(self):
|
||||
times[i] = result.time[0]
|
||||
activities[i] = result.get_material(mat_id).get_activity(units, by_nuclide, volume)
|
||||
activities[i] = result.get_material(mat_id).get_activity(
|
||||
units, by_nuclide, volume, chain_file=chain_file)
|
||||
|
||||
return times, activities
|
||||
|
||||
|
|
|
|||
|
|
@ -87,7 +87,7 @@ _dll.openmc_slice_data.errcheck = _error_handler
|
|||
|
||||
|
||||
def slice_data(origin, width=None, basis='xy', u_span=None, v_span=None,
|
||||
pixels=None, show_overlaps=False, level=-1, filter=None,
|
||||
pixels=None, show_overlaps=False, level=None, filter=None,
|
||||
include_properties=True):
|
||||
"""Generate a 2D raster of geometry and property data for plotting.
|
||||
|
||||
|
|
@ -111,7 +111,7 @@ def slice_data(origin, width=None, basis='xy', u_span=None, v_span=None,
|
|||
show_overlaps : bool, optional
|
||||
Whether to detect overlapping cells
|
||||
level : int, optional
|
||||
Universe level (-1 for deepest)
|
||||
Universe level (None for deepest)
|
||||
filter : openmc.lib.Filter, optional
|
||||
Filter for bin index lookup
|
||||
include_properties : bool, optional
|
||||
|
|
@ -127,6 +127,12 @@ def slice_data(origin, width=None, basis='xy', u_span=None, v_span=None,
|
|||
Array of shape (v_res, h_res, 2) with float64 dtype containing
|
||||
[temperature, density], or None if include_properties=False
|
||||
"""
|
||||
# Set deepest level as default
|
||||
if level is None:
|
||||
level = -1
|
||||
if not isinstance(level, int):
|
||||
raise TypeError("level must be an integer.")
|
||||
|
||||
if pixels is None:
|
||||
raise ValueError("pixels must be specified.")
|
||||
if len(pixels) != 2:
|
||||
|
|
@ -255,6 +261,52 @@ def property_map(plot):
|
|||
return prop_data
|
||||
|
||||
|
||||
_dll.openmc_slice_data_overlap_count.argtypes = [POINTER(c_size_t)]
|
||||
_dll.openmc_slice_data_overlap_count.restype = c_int
|
||||
_dll.openmc_slice_data_overlap_count.errcheck = _error_handler
|
||||
|
||||
_dll.openmc_slice_data_overlap_info.argtypes = [c_size_t, POINTER(c_int32)]
|
||||
_dll.openmc_slice_data_overlap_info.restype = c_int
|
||||
_dll.openmc_slice_data_overlap_info.errcheck = _error_handler
|
||||
|
||||
|
||||
# Python wrappings for overlap functions
|
||||
def slice_data_overlap_count() -> int:
|
||||
"""Return the number of unique overlaps from the last slice plot.
|
||||
|
||||
Returns
|
||||
-------
|
||||
int
|
||||
Number of unique overlapping cell pairs detected by the most recent
|
||||
:func:`slice_data` call with overlap checking enabled.
|
||||
"""
|
||||
count = c_size_t()
|
||||
_dll.openmc_slice_data_overlap_count(count)
|
||||
return count.value
|
||||
|
||||
|
||||
def slice_data_overlap_info() -> np.ndarray:
|
||||
"""Return identifying information for overlaps from the last slice plot.
|
||||
|
||||
Returns
|
||||
-------
|
||||
numpy.ndarray
|
||||
Array of shape ``(n, 3)`` with int32 dtype, where ``n`` is the number
|
||||
of unique overlaps detected by the most recent :func:`slice_data` call
|
||||
with overlap checking enabled. Each row contains ``[universe_id,
|
||||
cell1_id, cell2_id]``.
|
||||
"""
|
||||
n = slice_data_overlap_count()
|
||||
overlap_info = np.empty((n, 3), dtype=np.int32)
|
||||
|
||||
if n > 0:
|
||||
_dll.openmc_slice_data_overlap_info(
|
||||
n,
|
||||
overlap_info.ctypes.data_as(POINTER(c_int32)),
|
||||
)
|
||||
return overlap_info
|
||||
|
||||
|
||||
_dll.openmc_get_plot_index.argtypes = [c_int32, POINTER(c_int32)]
|
||||
_dll.openmc_get_plot_index.restype = c_int
|
||||
_dll.openmc_get_plot_index.errcheck = _error_handler
|
||||
|
|
|
|||
|
|
@ -8,7 +8,7 @@ from pathlib import Path
|
|||
import re
|
||||
import sys
|
||||
import tempfile
|
||||
from typing import Sequence, Dict
|
||||
from typing import TYPE_CHECKING, Literal, Sequence, Dict
|
||||
import warnings
|
||||
|
||||
import lxml.etree as ET
|
||||
|
|
@ -28,6 +28,9 @@ from openmc.data.data import _get_element_symbol, JOULE_PER_EV
|
|||
from openmc.data.function import Tabulated1D
|
||||
from openmc.data import mass_energy_absorption_coefficient, dose_coefficients
|
||||
|
||||
if TYPE_CHECKING:
|
||||
from openmc.deplete import Chain
|
||||
|
||||
|
||||
# Units for density supported by OpenMC
|
||||
DENSITY_UNITS = ('g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum',
|
||||
|
|
@ -1385,8 +1388,13 @@ class Material(IDManagerMixin):
|
|||
return densities
|
||||
|
||||
|
||||
def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False,
|
||||
volume: float | None = None) -> dict[str, float] | float:
|
||||
def get_activity(
|
||||
self,
|
||||
units: str = 'Bq/cm3',
|
||||
by_nuclide: bool = False,
|
||||
volume: float | None = None,
|
||||
chain_file: Literal[False] | None | PathLike | Chain = None
|
||||
) -> dict[str, float] | float:
|
||||
"""Return the activity of the material or each nuclide within.
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
|
@ -1405,13 +1413,22 @@ class Material(IDManagerMixin):
|
|||
:attr:`Material.volume` attribute.
|
||||
|
||||
.. versionadded:: 0.13.3
|
||||
chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
|
||||
Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
|
||||
used. If ``None``, the chain specified by
|
||||
``openmc.config['chain_file']`` is used when available. If a path or
|
||||
:class:`openmc.deplete.Chain` is given, that chain is used. For
|
||||
``None`` or an explicit chain, nuclides absent from the chain fall
|
||||
back to ENDF/B-VIII.0 data.
|
||||
|
||||
.. versionadded:: 0.15.4
|
||||
|
||||
Returns
|
||||
-------
|
||||
Union[dict, float]
|
||||
If by_nuclide is True then a dictionary whose keys are nuclide
|
||||
names and values are activity is returned. Otherwise the activity
|
||||
of the material is returned as a float.
|
||||
If by_nuclide is True then a dictionary whose keys are nuclide names
|
||||
and values are activity is returned. Otherwise the activity of the
|
||||
material is returned as a float.
|
||||
"""
|
||||
|
||||
cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3', 'Ci', 'Ci/m3'})
|
||||
|
|
@ -1440,7 +1457,8 @@ class Material(IDManagerMixin):
|
|||
|
||||
activity = {}
|
||||
for nuclide, atoms_per_bcm in self.get_nuclide_atom_densities().items():
|
||||
inv_seconds = openmc.data.decay_constant(nuclide)
|
||||
inv_seconds = openmc.data.decay_constant(
|
||||
nuclide, chain_file=chain_file)
|
||||
activity[nuclide] = inv_seconds * 1e24 * atoms_per_bcm * multiplier
|
||||
|
||||
return activity if by_nuclide else sum(activity.values())
|
||||
|
|
|
|||
|
|
@ -21,7 +21,8 @@ import openmc
|
|||
import openmc._xml as xml
|
||||
from openmc.dummy_comm import DummyCommunicator
|
||||
from openmc.executor import _process_CLI_arguments
|
||||
from openmc.checkvalue import check_type, check_value, PathLike
|
||||
from openmc.checkvalue import (check_type, check_value, check_greater_than,
|
||||
check_length, PathLike)
|
||||
from openmc.exceptions import InvalidIDError
|
||||
from openmc.plots import add_plot_params, _BASIS_INDICES, id_map_to_rgb
|
||||
from openmc.utility_funcs import change_directory
|
||||
|
|
@ -33,6 +34,15 @@ class ModelModifier(Protocol):
|
|||
...
|
||||
|
||||
|
||||
def _check_pixels(pixels: int | Sequence[int]) -> None:
|
||||
if isinstance(pixels, Integral):
|
||||
check_greater_than('pixels', pixels, 0)
|
||||
else:
|
||||
check_length('pixels', pixels, 2)
|
||||
for p in pixels:
|
||||
check_greater_than('pixels', p, 0)
|
||||
|
||||
|
||||
class Model:
|
||||
"""Model container.
|
||||
|
||||
|
|
@ -1050,6 +1060,8 @@ class Model:
|
|||
pixels: int | Sequence[int],
|
||||
basis: str
|
||||
):
|
||||
_check_pixels(pixels)
|
||||
|
||||
x, y, _ = _BASIS_INDICES[basis]
|
||||
|
||||
bb = self.bounding_box
|
||||
|
|
@ -1212,6 +1224,8 @@ class Model:
|
|||
"""
|
||||
import openmc.lib
|
||||
|
||||
_check_pixels(pixels)
|
||||
|
||||
if width is not None and (u_span is not None or v_span is not None):
|
||||
raise ValueError("width is mutually exclusive with u_span/v_span.")
|
||||
|
||||
|
|
@ -1301,7 +1315,11 @@ class Model:
|
|||
import matplotlib.pyplot as plt
|
||||
|
||||
check_type('n_samples', n_samples, int | None)
|
||||
if n_samples is not None:
|
||||
check_greater_than('n_samples', n_samples, 0, equality=True)
|
||||
check_type('plane_tolerance', plane_tolerance, Real)
|
||||
check_greater_than('plane_tolerance', plane_tolerance, 0.0)
|
||||
|
||||
if legend_kwargs is None:
|
||||
legend_kwargs = {}
|
||||
legend_kwargs.setdefault('bbox_to_anchor', (1.05, 1))
|
||||
|
|
|
|||
|
|
@ -383,7 +383,7 @@ def id_map_to_rgb(
|
|||
"""
|
||||
# Initialize RGB array with white background (values between 0 and 1 for matplotlib)
|
||||
img = np.ones(id_map.shape, dtype=float)
|
||||
|
||||
|
||||
# Get the appropriate index based on color_by
|
||||
if color_by == 'cell':
|
||||
id_index = 0 # Cell IDs are in the first channel
|
||||
|
|
@ -391,14 +391,14 @@ def id_map_to_rgb(
|
|||
id_index = 2 # Material IDs are in the third channel
|
||||
else:
|
||||
raise ValueError("color_by must be either 'cell' or 'material'")
|
||||
|
||||
|
||||
# Get all unique IDs in the plot
|
||||
unique_ids = np.unique(id_map[:, :, id_index])
|
||||
|
||||
|
||||
# Generate default colors if not provided
|
||||
if colors is None:
|
||||
colors = {}
|
||||
|
||||
|
||||
# Convert colors dict to use IDs as keys
|
||||
color_map = {}
|
||||
for key, color in colors.items():
|
||||
|
|
@ -406,13 +406,13 @@ def id_map_to_rgb(
|
|||
color_map[key.id] = color
|
||||
else:
|
||||
color_map[key] = color
|
||||
|
||||
|
||||
# Generate random colors for IDs not in color_map
|
||||
rng = np.random.RandomState(1)
|
||||
for uid in unique_ids:
|
||||
if uid > 0 and uid not in color_map:
|
||||
color_map[uid] = rng.randint(0, 256, (3,))
|
||||
|
||||
|
||||
# Apply colors to each pixel
|
||||
for uid in unique_ids:
|
||||
if uid == -1: # Background/void
|
||||
|
|
@ -432,7 +432,7 @@ def id_map_to_rgb(
|
|||
rgb = color
|
||||
mask = id_map[:, :, id_index] == uid
|
||||
img[mask] = np.array(rgb) / 255.0
|
||||
|
||||
|
||||
return img
|
||||
|
||||
class PlotBase(IDManagerMixin):
|
||||
|
|
|
|||
430
openmc/source.py
430
openmc/source.py
|
|
@ -1,7 +1,7 @@
|
|||
from __future__ import annotations
|
||||
from abc import ABC, abstractmethod
|
||||
from collections.abc import Iterable, Sequence
|
||||
from numbers import Real
|
||||
from numbers import Integral, Real
|
||||
from pathlib import Path
|
||||
import warnings
|
||||
from typing import Any
|
||||
|
|
@ -46,7 +46,7 @@ class SourceBase(ABC):
|
|||
|
||||
Attributes
|
||||
----------
|
||||
type : {'independent', 'file', 'compiled', 'mesh'}
|
||||
type : {'independent', 'file', 'compiled', 'mesh', 'tokamak'}
|
||||
Indicator of source type.
|
||||
strength : float
|
||||
Strength of the source
|
||||
|
|
@ -203,6 +203,8 @@ class SourceBase(ABC):
|
|||
return FileSource.from_xml_element(elem)
|
||||
elif source_type == 'mesh':
|
||||
return MeshSource.from_xml_element(elem, meshes)
|
||||
elif source_type == 'tokamak':
|
||||
return TokamakSource.from_xml_element(elem)
|
||||
else:
|
||||
raise ValueError(
|
||||
f'Source type {source_type} is not recognized')
|
||||
|
|
@ -896,6 +898,430 @@ class FileSource(SourceBase):
|
|||
return cls(**kwargs)
|
||||
|
||||
|
||||
class TokamakSource(SourceBase):
|
||||
r"""A source representing neutron emission from a tokamak plasma.
|
||||
|
||||
This source samples neutron positions from a tokamak plasma geometry using
|
||||
Miller-style flux surface parameterization. The user provides an emission
|
||||
profile S(r/a) as a function of normalized minor radius, along with one or
|
||||
more energy distributions.
|
||||
|
||||
The flux surface parameterization is
|
||||
|
||||
.. math::
|
||||
|
||||
\begin{aligned}
|
||||
R &= R_0 + r \cos\left(\alpha + \delta \sin\alpha\right)
|
||||
+ \Delta \left[1 - \left(\frac{r}{a}\right)^2\right] \\
|
||||
Z &= Z_\mathrm{shift} + \kappa r \sin\alpha
|
||||
\end{aligned}
|
||||
|
||||
where :math:`R_0` is major radius, :math:`a` is minor radius,
|
||||
:math:`\kappa` is elongation, :math:`\delta` is triangularity,
|
||||
:math:`\Delta` is the Shafranov shift, and :math:`Z_\mathrm{shift}` is
|
||||
the vertical shift.
|
||||
|
||||
.. versionadded:: 0.15.4
|
||||
|
||||
Parameters
|
||||
----------
|
||||
major_radius : float
|
||||
Major radius R0 in [cm]
|
||||
minor_radius : float
|
||||
Minor radius a in [cm]
|
||||
elongation : float
|
||||
Plasma elongation κ (must be > 0)
|
||||
triangularity : float
|
||||
Plasma triangularity δ (must be in [-1, 1])
|
||||
shafranov_shift : float
|
||||
Shafranov shift Δ in [cm] (must be >= 0 and < a/2)
|
||||
r_over_a : numpy.ndarray
|
||||
Normalized minor radius grid points, must start at 0 and end at 1
|
||||
emission_density : numpy.ndarray
|
||||
Emission density S(r) at each r/a point (arbitrary units, must be >= 0).
|
||||
Values are linearly interpolated between grid points and refined on an
|
||||
internal grid for radial sampling. Must have the same length as
|
||||
``r_over_a`` and contain at least one positive value.
|
||||
energy : openmc.stats.Univariate or Sequence[openmc.stats.Univariate]
|
||||
Energy distribution(s). Either a single distribution used at all radii,
|
||||
or one distribution per ``r_over_a`` grid point. When one distribution
|
||||
per grid point is given, the energy of a sampled particle is drawn from
|
||||
one of the two distributions bracketing its sampled radius, selected
|
||||
stochastically with probability proportional to the proximity of the
|
||||
radius to each grid point (stochastic interpolation).
|
||||
time : openmc.stats.Univariate, optional
|
||||
Time distribution of the source. If None, particles are born at
|
||||
:math:`t=0`, matching the default behavior of
|
||||
:class:`openmc.IndependentSource`.
|
||||
phi_start : float
|
||||
Starting toroidal angle in [rad] (default: 0)
|
||||
phi_extent : float
|
||||
Toroidal angle extent in [rad] (default: 2π)
|
||||
n_alpha : int
|
||||
Number of poloidal angle grid points for CDF sampling (default: 101)
|
||||
vertical_shift : float
|
||||
Vertical shift of the plasma center in [cm] (default: 0)
|
||||
strength : float
|
||||
Strength of the source (default: 1.0)
|
||||
constraints : dict
|
||||
Constraints on sampled source particles. See :class:`SourceBase` for
|
||||
valid keys and values.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
major_radius : float
|
||||
Major radius R0 in [cm]
|
||||
minor_radius : float
|
||||
Minor radius a in [cm]
|
||||
elongation : float
|
||||
Plasma elongation κ
|
||||
triangularity : float
|
||||
Plasma triangularity δ
|
||||
shafranov_shift : float
|
||||
Shafranov shift Δ in [cm]
|
||||
r_over_a : numpy.ndarray
|
||||
Normalized minor radius grid points
|
||||
emission_density : numpy.ndarray
|
||||
Emission density S(r) at each r/a point
|
||||
energy : list of openmc.stats.Univariate
|
||||
Energy distribution(s)
|
||||
time : openmc.stats.Univariate or None
|
||||
Time distribution of the source
|
||||
phi_start : float
|
||||
Starting toroidal angle in [rad]
|
||||
phi_extent : float
|
||||
Toroidal angle extent in [rad]
|
||||
n_alpha : int
|
||||
Number of poloidal angle grid points
|
||||
vertical_shift : float
|
||||
Vertical shift of the plasma center in [cm]
|
||||
strength : float
|
||||
Strength of the source
|
||||
type : str
|
||||
Indicator of source type: 'tokamak'
|
||||
constraints : dict
|
||||
Constraints on sampled source particles
|
||||
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
major_radius: float,
|
||||
minor_radius: float,
|
||||
elongation: float,
|
||||
triangularity: float,
|
||||
shafranov_shift: float,
|
||||
r_over_a: Sequence[float],
|
||||
emission_density: Sequence[float],
|
||||
energy: Univariate | Sequence[Univariate],
|
||||
time: Univariate | None = None,
|
||||
phi_start: float = 0.0,
|
||||
phi_extent: float = 2.0 * np.pi,
|
||||
n_alpha: int = 101,
|
||||
vertical_shift: float = 0.0,
|
||||
strength: float = 1.0,
|
||||
constraints: dict[str, Any] | None = None
|
||||
):
|
||||
super().__init__(strength=strength, constraints=constraints)
|
||||
self.major_radius = major_radius
|
||||
self.minor_radius = minor_radius
|
||||
self.elongation = elongation
|
||||
self.triangularity = triangularity
|
||||
self.shafranov_shift = shafranov_shift
|
||||
self.r_over_a = r_over_a
|
||||
self.emission_density = emission_density
|
||||
self.phi_start = phi_start
|
||||
self.phi_extent = phi_extent
|
||||
self.n_alpha = n_alpha
|
||||
self.vertical_shift = vertical_shift
|
||||
self.energy = energy
|
||||
self.time = time
|
||||
|
||||
self._validate()
|
||||
|
||||
def _validate(self):
|
||||
"""Validate relationships between tokamak source parameters."""
|
||||
if self.minor_radius >= self.major_radius:
|
||||
raise ValueError(
|
||||
f"minor_radius ({self.minor_radius}) must be smaller than "
|
||||
f"major_radius ({self.major_radius})")
|
||||
if self.shafranov_shift >= 0.5 * self.minor_radius:
|
||||
raise ValueError(
|
||||
f"shafranov_shift ({self.shafranov_shift}) must be smaller "
|
||||
f"than half the minor_radius ({0.5 * self.minor_radius})")
|
||||
if len(self.emission_density) != len(self.r_over_a):
|
||||
raise ValueError(
|
||||
f"emission_density (length {len(self.emission_density)}) must "
|
||||
f"have the same length as r_over_a (length {len(self.r_over_a)})")
|
||||
if not np.any(self.emission_density > 0.0):
|
||||
raise ValueError("emission_density must contain a positive value")
|
||||
if len(self.energy) not in (1, len(self.r_over_a)):
|
||||
raise ValueError(
|
||||
f"Number of energy distributions ({len(self.energy)}) must be "
|
||||
f"either 1 or equal to the number of r_over_a grid points "
|
||||
f"({len(self.r_over_a)})")
|
||||
|
||||
@property
|
||||
def type(self) -> str:
|
||||
return "tokamak"
|
||||
|
||||
@property
|
||||
def major_radius(self) -> float:
|
||||
return self._major_radius
|
||||
|
||||
@major_radius.setter
|
||||
def major_radius(self, value: float):
|
||||
cv.check_type('major radius', value, Real)
|
||||
cv.check_greater_than('major radius', value, 0.0)
|
||||
self._major_radius = value
|
||||
|
||||
@property
|
||||
def minor_radius(self) -> float:
|
||||
return self._minor_radius
|
||||
|
||||
@minor_radius.setter
|
||||
def minor_radius(self, value: float):
|
||||
cv.check_type('minor radius', value, Real)
|
||||
cv.check_greater_than('minor radius', value, 0.0)
|
||||
self._minor_radius = value
|
||||
|
||||
@property
|
||||
def elongation(self) -> float:
|
||||
return self._elongation
|
||||
|
||||
@elongation.setter
|
||||
def elongation(self, value: float):
|
||||
cv.check_type('elongation', value, Real)
|
||||
cv.check_greater_than('elongation', value, 0.0)
|
||||
self._elongation = value
|
||||
|
||||
@property
|
||||
def triangularity(self) -> float:
|
||||
return self._triangularity
|
||||
|
||||
@triangularity.setter
|
||||
def triangularity(self, value: float):
|
||||
cv.check_type('triangularity', value, Real)
|
||||
cv.check_greater_than('triangularity', value, -1.0, equality=True)
|
||||
cv.check_less_than('triangularity', value, 1.0, equality=True)
|
||||
self._triangularity = value
|
||||
|
||||
@property
|
||||
def shafranov_shift(self) -> float:
|
||||
return self._shafranov_shift
|
||||
|
||||
@shafranov_shift.setter
|
||||
def shafranov_shift(self, value: float):
|
||||
cv.check_type('Shafranov shift', value, Real)
|
||||
cv.check_greater_than('Shafranov shift', value, 0.0, equality=True)
|
||||
self._shafranov_shift = value
|
||||
|
||||
@property
|
||||
def r_over_a(self) -> np.ndarray:
|
||||
return self._r_over_a
|
||||
|
||||
@r_over_a.setter
|
||||
def r_over_a(self, value: Sequence[float]):
|
||||
value = np.asarray(value, dtype=float)
|
||||
if value.ndim != 1 or len(value) < 2:
|
||||
raise ValueError("r_over_a must be a 1-D array with at least 2 points")
|
||||
if value[0] != 0.0:
|
||||
raise ValueError("r_over_a must start at 0")
|
||||
if value[-1] != 1.0:
|
||||
raise ValueError("r_over_a must end at 1")
|
||||
if not np.all(np.diff(value) > 0):
|
||||
raise ValueError("r_over_a must be strictly increasing")
|
||||
self._r_over_a = value
|
||||
|
||||
@property
|
||||
def emission_density(self) -> np.ndarray:
|
||||
return self._emission_density
|
||||
|
||||
@emission_density.setter
|
||||
def emission_density(self, value: Sequence[float]):
|
||||
value = np.asarray(value, dtype=float)
|
||||
if value.ndim != 1:
|
||||
raise ValueError("emission_density must be a 1-D array")
|
||||
if np.any(value < 0):
|
||||
raise ValueError("emission_density values cannot be negative")
|
||||
self._emission_density = value
|
||||
|
||||
@property
|
||||
def energy(self) -> list[Univariate]:
|
||||
return self._energy
|
||||
|
||||
@energy.setter
|
||||
def energy(self, value: Univariate | Sequence[Univariate]):
|
||||
if isinstance(value, Univariate):
|
||||
self._energy = [value]
|
||||
else:
|
||||
cv.check_iterable_type('energy distributions', value, Univariate)
|
||||
self._energy = list(value)
|
||||
|
||||
@property
|
||||
def time(self) -> Univariate | None:
|
||||
return self._time
|
||||
|
||||
@time.setter
|
||||
def time(self, value: Univariate | None):
|
||||
if value is not None:
|
||||
cv.check_type('time distribution', value, Univariate)
|
||||
self._time = value
|
||||
|
||||
@property
|
||||
def phi_start(self) -> float:
|
||||
return self._phi_start
|
||||
|
||||
@phi_start.setter
|
||||
def phi_start(self, value: float):
|
||||
cv.check_type('phi_start', value, Real)
|
||||
self._phi_start = value
|
||||
|
||||
@property
|
||||
def phi_extent(self) -> float:
|
||||
return self._phi_extent
|
||||
|
||||
@phi_extent.setter
|
||||
def phi_extent(self, value: float):
|
||||
cv.check_type('phi_extent', value, Real)
|
||||
cv.check_greater_than('phi_extent', value, 0.0)
|
||||
cv.check_less_than('phi_extent', value, 2.0 * np.pi, equality=True)
|
||||
self._phi_extent = value
|
||||
|
||||
@property
|
||||
def n_alpha(self) -> int:
|
||||
return self._n_alpha
|
||||
|
||||
@n_alpha.setter
|
||||
def n_alpha(self, value: int):
|
||||
cv.check_type('n_alpha', value, Integral)
|
||||
cv.check_greater_than('n_alpha', value, 2)
|
||||
if value < 51:
|
||||
warnings.warn(
|
||||
"n_alpha values below 51 may introduce noticeable "
|
||||
"discretization bias in tokamak source sampling", stacklevel=2)
|
||||
self._n_alpha = value
|
||||
|
||||
@property
|
||||
def vertical_shift(self) -> float:
|
||||
return self._vertical_shift
|
||||
|
||||
@vertical_shift.setter
|
||||
def vertical_shift(self, value: float):
|
||||
cv.check_type('vertical shift', value, Real)
|
||||
self._vertical_shift = value
|
||||
|
||||
def populate_xml_element(self, element):
|
||||
"""Add necessary tokamak source information to an XML element
|
||||
|
||||
Returns
|
||||
-------
|
||||
element : lxml.etree._Element
|
||||
XML element containing source data
|
||||
|
||||
"""
|
||||
self._validate()
|
||||
|
||||
# Geometry parameters
|
||||
ET.SubElement(element, "major_radius").text = str(self.major_radius)
|
||||
ET.SubElement(element, "minor_radius").text = str(self.minor_radius)
|
||||
ET.SubElement(element, "elongation").text = str(self.elongation)
|
||||
ET.SubElement(element, "triangularity").text = str(self.triangularity)
|
||||
ET.SubElement(element, "shafranov_shift").text = str(self.shafranov_shift)
|
||||
|
||||
# Toroidal angle bounds
|
||||
ET.SubElement(element, "phi_start").text = str(self.phi_start)
|
||||
ET.SubElement(element, "phi_extent").text = str(self.phi_extent)
|
||||
|
||||
# Poloidal sampling resolution
|
||||
ET.SubElement(element, "n_alpha").text = str(self.n_alpha)
|
||||
|
||||
# Vertical shift
|
||||
if self.vertical_shift != 0.0:
|
||||
ET.SubElement(element, "vertical_shift").text = str(self.vertical_shift)
|
||||
|
||||
# Emission profile
|
||||
ET.SubElement(element, "r_over_a").text = ' '.join(str(r) for r in self.r_over_a)
|
||||
ET.SubElement(element, "emission_density").text = ' '.join(str(s) for s in self.emission_density)
|
||||
|
||||
# Energy distribution(s)
|
||||
for dist in self.energy:
|
||||
element.append(dist.to_xml_element('energy'))
|
||||
|
||||
# Time distribution
|
||||
if self.time is not None:
|
||||
element.append(self.time.to_xml_element('time'))
|
||||
|
||||
@classmethod
|
||||
def from_xml_element(cls, elem: ET.Element) -> TokamakSource:
|
||||
"""Generate tokamak source from an XML element
|
||||
|
||||
Parameters
|
||||
----------
|
||||
elem : lxml.etree._Element
|
||||
XML element
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.TokamakSource
|
||||
Source generated from XML element
|
||||
|
||||
"""
|
||||
# Read geometry parameters
|
||||
major_radius = float(get_text(elem, 'major_radius'))
|
||||
minor_radius = float(get_text(elem, 'minor_radius'))
|
||||
elongation = float(get_text(elem, 'elongation'))
|
||||
triangularity = float(get_text(elem, 'triangularity'))
|
||||
shafranov_shift = float(get_text(elem, 'shafranov_shift'))
|
||||
|
||||
# Read optional parameters
|
||||
phi_start_text = get_text(elem, 'phi_start')
|
||||
phi_start = float(phi_start_text) if phi_start_text else 0.0
|
||||
|
||||
phi_extent_text = get_text(elem, 'phi_extent')
|
||||
phi_extent = float(phi_extent_text) if phi_extent_text else 2.0 * np.pi
|
||||
|
||||
n_alpha_text = get_text(elem, 'n_alpha')
|
||||
n_alpha = int(n_alpha_text) if n_alpha_text else 101
|
||||
|
||||
vertical_shift_text = get_text(elem, 'vertical_shift')
|
||||
vertical_shift = float(vertical_shift_text) if vertical_shift_text else 0.0
|
||||
|
||||
# Read emission profile
|
||||
r_over_a = np.array([float(x) for x in get_text(elem, 'r_over_a').split()])
|
||||
emission_density = np.array([float(x) for x in get_text(elem, 'emission_density').split()])
|
||||
|
||||
# Read energy distributions
|
||||
energy = [Univariate.from_xml_element(e) for e in elem.findall('energy')]
|
||||
if len(energy) == 1:
|
||||
energy = energy[0]
|
||||
|
||||
# Read time distribution
|
||||
time_elem = elem.find('time')
|
||||
time = Univariate.from_xml_element(time_elem) if time_elem is not None else None
|
||||
|
||||
# Read constraints and strength
|
||||
constraints = cls._get_constraints(elem)
|
||||
strength_text = get_text(elem, 'strength')
|
||||
strength = float(strength_text) if strength_text else 1.0
|
||||
|
||||
return cls(
|
||||
major_radius=major_radius,
|
||||
minor_radius=minor_radius,
|
||||
elongation=elongation,
|
||||
triangularity=triangularity,
|
||||
shafranov_shift=shafranov_shift,
|
||||
r_over_a=r_over_a,
|
||||
emission_density=emission_density,
|
||||
energy=energy,
|
||||
time=time,
|
||||
phi_start=phi_start,
|
||||
phi_extent=phi_extent,
|
||||
n_alpha=n_alpha,
|
||||
vertical_shift=vertical_shift,
|
||||
strength=strength,
|
||||
constraints=constraints
|
||||
)
|
||||
|
||||
|
||||
class SourceParticle:
|
||||
|
|
|
|||
|
|
@ -40,6 +40,7 @@ dependencies = [
|
|||
[project.optional-dependencies]
|
||||
depletion-mpi = ["mpi4py"]
|
||||
docs = [
|
||||
"breathe",
|
||||
"sphinx",
|
||||
"sphinxcontrib-katex",
|
||||
"sphinx-numfig",
|
||||
|
|
|
|||
|
|
@ -38,6 +38,9 @@ void VacuumBC::handle_particle(Particle& p, const Surface& surf) const
|
|||
void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const
|
||||
{
|
||||
Direction u = surf.reflect(p.r(), p.u(), &p);
|
||||
|
||||
// normalize reflected u to ensure no floating point error leads to
|
||||
// unnormalized directions
|
||||
u /= u.norm();
|
||||
|
||||
// Handle the effects of the surface albedo on the particle's weight.
|
||||
|
|
@ -53,6 +56,9 @@ void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const
|
|||
void WhiteBC::handle_particle(Particle& p, const Surface& surf) const
|
||||
{
|
||||
Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed());
|
||||
|
||||
// normalize outgoing u to ensure no floating point error leads to
|
||||
// unnormalized directions
|
||||
u /= u.norm();
|
||||
|
||||
// Handle the effects of the surface albedo on the particle's weight.
|
||||
|
|
@ -241,6 +247,10 @@ void RotationalPeriodicBC::handle_particle(
|
|||
new_u[axis_1_idx_] = cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_];
|
||||
new_u[axis_2_idx_] = sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_];
|
||||
|
||||
// normalize new_u to ensure no floating point error leads to unnormalized
|
||||
// directions
|
||||
new_u /= new_u.norm();
|
||||
|
||||
// Handle the effects of the surface albedo on the particle's weight.
|
||||
BoundaryCondition::handle_albedo(p, surf);
|
||||
|
||||
|
|
|
|||
|
|
@ -552,14 +552,9 @@ double Tabular::sample_unbiased(uint64_t* seed) const
|
|||
double c = prn(seed);
|
||||
|
||||
// Find first CDF bin which is above the sampled value
|
||||
double c_i = c_[0];
|
||||
int i;
|
||||
std::size_t n = c_.size();
|
||||
for (i = 0; i < n - 1; ++i) {
|
||||
if (c <= c_[i + 1])
|
||||
break;
|
||||
c_i = c_[i + 1];
|
||||
}
|
||||
auto c_iter = std::lower_bound(c_.begin() + 1, c_.end(), c);
|
||||
int i = std::distance(c_.begin(), c_iter) - 1;
|
||||
double c_i = c_[i];
|
||||
|
||||
// Determine bounding PDF values
|
||||
double x_i = x_[i];
|
||||
|
|
|
|||
|
|
@ -262,9 +262,11 @@ std::pair<Position, double> SphericalIndependent::sample(uint64_t* seed) const
|
|||
MeshSpatial::MeshSpatial(pugi::xml_node node)
|
||||
{
|
||||
|
||||
if (get_node_value(node, "type", true, true) != "mesh") {
|
||||
fatal_error(fmt::format(
|
||||
"Incorrect spatial type '{}' for a MeshSpatial distribution"));
|
||||
auto spatial_type = get_node_value(node, "type", true, true);
|
||||
if (spatial_type != "mesh") {
|
||||
fatal_error(
|
||||
fmt::format("Incorrect spatial type '{}' for a MeshSpatial distribution",
|
||||
spatial_type));
|
||||
}
|
||||
|
||||
// No in-tet distributions implemented, could include distributions for the
|
||||
|
|
|
|||
|
|
@ -26,16 +26,22 @@ int n_coord_levels;
|
|||
|
||||
vector<int64_t> overlap_check_count;
|
||||
|
||||
vector<OverlapKey> overlap_keys;
|
||||
std::unordered_map<OverlapKey, int, OverlapKeyHash> overlap_key_index;
|
||||
|
||||
} // namespace model
|
||||
|
||||
//==============================================================================
|
||||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
||||
bool check_cell_overlap(GeometryState& p, bool error)
|
||||
int check_cell_overlap(GeometryState& p, bool error)
|
||||
{
|
||||
int n_coord = p.n_coord();
|
||||
|
||||
// If no overlap found, return a nonphysical index
|
||||
int overlap_index = -1;
|
||||
|
||||
// Loop through each coordinate level
|
||||
for (int j = 0; j < n_coord; j++) {
|
||||
Universe& univ = *model::universes[p.coord(j).universe()];
|
||||
|
|
@ -44,21 +50,40 @@ bool check_cell_overlap(GeometryState& p, bool error)
|
|||
for (auto index_cell : univ.cells_) {
|
||||
Cell& c = *model::cells[index_cell];
|
||||
if (c.contains(p.coord(j).r(), p.coord(j).u(), p.surface())) {
|
||||
#pragma omp atomic
|
||||
++model::overlap_check_count[index_cell];
|
||||
if (index_cell != p.coord(j).cell()) {
|
||||
if (error) {
|
||||
fatal_error(
|
||||
fmt::format("Overlapping cells detected: {}, {} on universe {}",
|
||||
c.id_, model::cells[p.coord(j).cell()]->id_, univ.id_));
|
||||
}
|
||||
return true;
|
||||
|
||||
// With no fatal error (plotter is calling), now adds overlaps and
|
||||
// ensures order does not matter when making overlap key
|
||||
int cell_a = model::cells[index_cell]->id_;
|
||||
int cell_b = model::cells[p.coord(j).cell()]->id_;
|
||||
int a = std::min(cell_a, cell_b);
|
||||
int b = std::max(cell_a, cell_b);
|
||||
OverlapKey key {univ.id_, a, b};
|
||||
#pragma omp critical(overlap_key_update)
|
||||
{
|
||||
auto it = model::overlap_key_index.find(key);
|
||||
if (it != model::overlap_key_index.end()) {
|
||||
overlap_index = it->second; // already exists, reuse index
|
||||
} else {
|
||||
int idx = int(model::overlap_keys.size());
|
||||
model::overlap_keys.push_back(key);
|
||||
model::overlap_key_index[key] = idx;
|
||||
overlap_index = idx;
|
||||
}
|
||||
}
|
||||
break;
|
||||
}
|
||||
#pragma omp atomic
|
||||
++model::overlap_check_count[index_cell];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return false;
|
||||
return overlap_index;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
|
|||
|
|
@ -10,6 +10,7 @@
|
|||
#include "openmc/geometry.h"
|
||||
#include "openmc/geometry_aux.h"
|
||||
#include "openmc/hdf5_interface.h"
|
||||
#include "openmc/math_functions.h"
|
||||
#include "openmc/string_utils.h"
|
||||
#include "openmc/vector.h"
|
||||
#include "openmc/xml_interface.h"
|
||||
|
|
@ -260,25 +261,33 @@ std::pair<double, array<int, 3>> RectLattice::distance(
|
|||
// Determine the oncoming edge.
|
||||
double x0 {copysign(0.5 * pitch_[0], u.x)};
|
||||
double y0 {copysign(0.5 * pitch_[1], u.y)};
|
||||
double z0;
|
||||
|
||||
double d = std::min(
|
||||
u.x != 0.0 ? (x0 - x) / u.x : INFTY, u.y != 0.0 ? (y0 - y) / u.y : INFTY);
|
||||
// Evaluate the distance to each oncoming edge independently. Comparing these
|
||||
// distances directly (rather than reconstructing the crossing position)
|
||||
// avoids the floating-point cancellation that occurs for large pitches.
|
||||
double dx = u.x != 0.0 ? (x0 - x) / u.x : INFTY;
|
||||
double dy = u.y != 0.0 ? (y0 - y) / u.y : INFTY;
|
||||
double dz = INFTY;
|
||||
if (is_3d_) {
|
||||
z0 = copysign(0.5 * pitch_[2], u.z);
|
||||
d = std::min(d, u.z != 0.0 ? (z0 - z) / u.z : INFTY);
|
||||
double z0 {copysign(0.5 * pitch_[2], u.z)};
|
||||
dz = u.z != 0.0 ? (z0 - z) / u.z : INFTY;
|
||||
}
|
||||
|
||||
// Determine which lattice boundaries are being crossed
|
||||
// The distance to the nearest lattice boundary is the smallest axial
|
||||
// distance.
|
||||
double d = std::min({dx, dy, dz});
|
||||
|
||||
// Determine which lattice boundaries are being crossed. The axis attaining
|
||||
// the minimum is exactly equal to d, so at least one translation is always
|
||||
// set for a finite crossing; a near-equal second axis indicates a corner
|
||||
// crossing.
|
||||
array<int, 3> lattice_trans = {0, 0, 0};
|
||||
if (u.x != 0.0 && std::abs(x + u.x * d - x0) < FP_PRECISION)
|
||||
if (isclose(d, dx, FP_COINCIDENT, FP_PRECISION))
|
||||
lattice_trans[0] = copysign(1, u.x);
|
||||
if (u.y != 0.0 && std::abs(y + u.y * d - y0) < FP_PRECISION)
|
||||
if (isclose(d, dy, FP_COINCIDENT, FP_PRECISION))
|
||||
lattice_trans[1] = copysign(1, u.y);
|
||||
if (is_3d_) {
|
||||
if (u.z != 0.0 && std::abs(z + u.z * d - z0) < FP_PRECISION)
|
||||
lattice_trans[2] = copysign(1, u.z);
|
||||
}
|
||||
if (is_3d_ && isclose(d, dz, FP_COINCIDENT, FP_PRECISION))
|
||||
lattice_trans[2] = copysign(1, u.z);
|
||||
|
||||
return {d, lattice_trans};
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,5 +1,7 @@
|
|||
#include "openmc/math_functions.h"
|
||||
|
||||
#include <limits> // for numeric_limits
|
||||
|
||||
#include "openmc/external/Faddeeva.hh"
|
||||
|
||||
#include "openmc/constants.h"
|
||||
|
|
@ -944,6 +946,46 @@ double log1prel(double x)
|
|||
}
|
||||
}
|
||||
|
||||
double cyl_bessel_j(int n, double x)
|
||||
{
|
||||
// Handle negative arguments via the parity relation
|
||||
// J_n(-x) = (-1)^n J_n(x); std::cyl_bessel_j has a domain error for x < 0.
|
||||
double sign = 1.0;
|
||||
if (x < 0.0) {
|
||||
x = -x;
|
||||
if (n % 2 == 1)
|
||||
sign = -1.0;
|
||||
}
|
||||
|
||||
#if defined(__cpp_lib_math_special_functions) && \
|
||||
__cpp_lib_math_special_functions >= 201603L
|
||||
return sign * std::cyl_bessel_j(static_cast<double>(n), x);
|
||||
#else
|
||||
// Ascending power series (e.g., Abramowitz & Stegun eq. 9.1.10):
|
||||
// J_n(x) = sum_{m=0}^inf (-1)^m / (m! (m+n)!) * (x/2)^(2m+n)
|
||||
// The term ratio is -(x/2)^2 / (m*(m+n)), so for |x| <= 2 the series
|
||||
// converges to machine precision within ~20 terms.
|
||||
double half_x = 0.5 * x;
|
||||
|
||||
// First term: (x/2)^n / n!
|
||||
double term = 1.0;
|
||||
for (int k = 1; k <= n; ++k) {
|
||||
term *= half_x / k;
|
||||
}
|
||||
|
||||
double sum = term;
|
||||
double neg_half_x_sq = -half_x * half_x;
|
||||
for (int m = 1; m <= 50; ++m) {
|
||||
term *= neg_half_x_sq / (m * (m + n));
|
||||
sum += term;
|
||||
if (std::abs(term) <=
|
||||
std::numeric_limits<double>::epsilon() * std::abs(sum))
|
||||
break;
|
||||
}
|
||||
return sign * sum;
|
||||
#endif
|
||||
}
|
||||
|
||||
// Helper function to get index and interpolation function on an incident energy
|
||||
// grid
|
||||
void get_energy_index(
|
||||
|
|
@ -959,4 +1001,12 @@ void get_energy_index(
|
|||
}
|
||||
}
|
||||
|
||||
// Return true if two floating-point values are approximately equal within a
|
||||
// combined relative and absolute tolerance.
|
||||
bool isclose(double a, double b, double rel_tol, double abs_tol)
|
||||
{
|
||||
return std::abs(a - b) <=
|
||||
std::max(rel_tol * std::max(std::abs(a), std::abs(b)), abs_tol);
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -338,13 +338,14 @@ void Particle::event_cross_surface()
|
|||
boundary().lattice_translation()[2] != 0) {
|
||||
// Particle crosses lattice boundary
|
||||
|
||||
int i_lattice = coord(boundary().coord_level() - 1).lattice();
|
||||
bool verbose = settings::verbosity >= 10 || trace();
|
||||
cross_lattice(*this, boundary(), verbose);
|
||||
event() = TallyEvent::LATTICE;
|
||||
|
||||
// Score cell to cell partial currents
|
||||
if (!model::active_surface_tallies.empty()) {
|
||||
auto& lat {*model::lattices[lowest_coord().lattice()]};
|
||||
auto& lat {*model::lattices[i_lattice]};
|
||||
bool is_valid;
|
||||
Direction normal =
|
||||
lat.get_normal(boundary().lattice_translation(), is_valid);
|
||||
|
|
@ -553,15 +554,30 @@ void Particle::event_death()
|
|||
finalize_particle_track(*this);
|
||||
}
|
||||
|
||||
// Contribute tally reduction variables to global accumulator
|
||||
// Contribute tally reduction variables to global accumulator
|
||||
const auto k_absorption = keff_tally_absorption();
|
||||
const auto k_collision = keff_tally_collision();
|
||||
const auto k_tracklength = keff_tally_tracklength();
|
||||
const auto leakage = keff_tally_leakage();
|
||||
|
||||
if (settings::run_mode == RunMode::EIGENVALUE) {
|
||||
if (k_absorption != 0.0) {
|
||||
#pragma omp atomic
|
||||
global_tally_absorption += keff_tally_absorption();
|
||||
global_tally_absorption += k_absorption;
|
||||
}
|
||||
if (k_collision != 0.0) {
|
||||
#pragma omp atomic
|
||||
global_tally_collision += keff_tally_collision();
|
||||
global_tally_collision += k_collision;
|
||||
}
|
||||
if (k_tracklength != 0.0) {
|
||||
#pragma omp atomic
|
||||
global_tally_tracklength += keff_tally_tracklength();
|
||||
global_tally_tracklength += k_tracklength;
|
||||
}
|
||||
}
|
||||
if (leakage != 0.0) {
|
||||
#pragma omp atomic
|
||||
global_tally_leakage += keff_tally_leakage();
|
||||
global_tally_leakage += leakage;
|
||||
}
|
||||
|
||||
// Reset particle tallies once accumulated
|
||||
keff_tally_absorption() = 0.0;
|
||||
|
|
|
|||
|
|
@ -359,8 +359,8 @@ void sample_photon_reaction(Particle& p)
|
|||
// Allow electrons to fill orbital and produce Auger electrons and
|
||||
// fluorescent photons. Since Compton subshell data does not match atomic
|
||||
// relaxation data, use the mapping between the data to find the subshell
|
||||
if (settings::atomic_relaxation && i_shell >= 0 &&
|
||||
element.subshell_map_[i_shell] >= 0) {
|
||||
if (settings::atomic_relaxation && element.has_atomic_relaxation_ &&
|
||||
i_shell >= 0 && element.subshell_map_[i_shell] >= 0) {
|
||||
element.atomic_relaxation(element.subshell_map_[i_shell], p);
|
||||
}
|
||||
|
||||
|
|
|
|||
53
src/plot.cpp
53
src/plot.cpp
|
|
@ -73,7 +73,7 @@ void IdData::set_value(size_t y, size_t x, const Particle& p, int level,
|
|||
}
|
||||
}
|
||||
|
||||
void IdData::set_overlap(size_t y, size_t x)
|
||||
void IdData::set_overlap(size_t y, size_t x, int /*overlap_idx*/)
|
||||
{
|
||||
for (size_t k = 0; k < data_.shape(2); ++k)
|
||||
data_(y, x, k) = OVERLAP;
|
||||
|
|
@ -88,13 +88,10 @@ void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level,
|
|||
{
|
||||
Cell* c = model::cells.at(p.lowest_coord().cell()).get();
|
||||
data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN;
|
||||
if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) {
|
||||
Material* m = model::materials.at(p.material()).get();
|
||||
data_(y, x, 1) = m->density_gpcc_;
|
||||
}
|
||||
data_(y, x, 1) = c->density(p.cell_instance());
|
||||
}
|
||||
|
||||
void PropertyData::set_overlap(size_t y, size_t x)
|
||||
void PropertyData::set_overlap(size_t y, size_t x, int /*overlap_idx*/)
|
||||
{
|
||||
data_(y, x) = OVERLAP;
|
||||
}
|
||||
|
|
@ -150,18 +147,18 @@ void RasterData::set_value(size_t y, size_t x, const Particle& p, int level,
|
|||
// set density (g/cm³)
|
||||
if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) {
|
||||
Material* m = model::materials.at(p.material()).get();
|
||||
property_data_(y, x, 1) = m->density_gpcc_;
|
||||
property_data_(y, x, 1) = c->density(p.cell_instance());
|
||||
}
|
||||
}
|
||||
|
||||
void RasterData::set_overlap(size_t y, size_t x)
|
||||
void RasterData::set_overlap(size_t y, size_t x, int overlap_idx)
|
||||
{
|
||||
// Set cell, instance, and material to OVERLAP, but preserve filter bin
|
||||
id_data_(y, x, 0) = OVERLAP;
|
||||
// Set cell, instance, and material to OVERLAP, but preserve filter bin for
|
||||
// tally plotting. Cell encodes the overlap index as a negative number so that
|
||||
// it can be used to look up overlap information in the plotter.
|
||||
id_data_(y, x, 0) = OVERLAP - overlap_idx - 1;
|
||||
id_data_(y, x, 1) = OVERLAP;
|
||||
id_data_(y, x, 2) = OVERLAP;
|
||||
// Note: id_data_(y, x, 3) is NOT overwritten - preserves filter bin for tally
|
||||
// plotting
|
||||
|
||||
property_data_(y, x, 0) = OVERLAP;
|
||||
property_data_(y, x, 1) = OVERLAP;
|
||||
|
|
@ -1991,10 +1988,12 @@ extern "C" int openmc_slice_data(const double origin[3], const double u_span[3],
|
|||
plot_params.show_overlaps_ = color_overlaps;
|
||||
plot_params.slice_level_ = level;
|
||||
|
||||
// Clear overlap data structures on new slice call
|
||||
model::overlap_keys.clear();
|
||||
model::overlap_key_index.clear();
|
||||
|
||||
// Use get_map<RasterData> to generate data
|
||||
auto data = plot_params.get_map<RasterData>(filter_index);
|
||||
|
||||
// Copy geometry data
|
||||
std::copy(data.id_data_.begin(), data.id_data_.end(), geom_data);
|
||||
|
||||
// Copy property data if requested
|
||||
|
|
@ -2010,6 +2009,32 @@ extern "C" int openmc_slice_data(const double origin[3], const double u_span[3],
|
|||
return 0;
|
||||
}
|
||||
|
||||
// Gets the number of overlaps that we need data for
|
||||
extern "C" int openmc_slice_data_overlap_count(size_t* count)
|
||||
{
|
||||
if (!count) {
|
||||
set_errmsg("Null pointer passed for overlap count.");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
*count = model::overlap_keys.size();
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
// Plotter pre-allocates array size based on what is returned with
|
||||
// overlap_count; populates an array of size 3*count
|
||||
extern "C" int openmc_slice_data_overlap_info(
|
||||
size_t count, int32_t* overlap_info)
|
||||
{
|
||||
for (size_t i = 0; i < count; ++i) {
|
||||
overlap_info[i * 3] = model::overlap_keys[i].universe_id;
|
||||
overlap_info[i * 3 + 1] = model::overlap_keys[i].cell1_id;
|
||||
overlap_info[i * 3 + 2] = model::overlap_keys[i].cell2_id;
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
extern "C" int openmc_get_plot_index(int32_t id, int32_t* index)
|
||||
{
|
||||
auto it = model::plot_map.find(id);
|
||||
|
|
|
|||
|
|
@ -12,6 +12,45 @@ constexpr uint64_t prn_mult {6364136223846793005ULL}; // multiplication
|
|||
constexpr uint64_t prn_add {1442695040888963407ULL}; // additive factor, c
|
||||
uint64_t prn_stride {DEFAULT_STRIDE}; // stride between particles
|
||||
|
||||
namespace {
|
||||
|
||||
struct SkipAheadCoefficients {
|
||||
uint64_t multiplier;
|
||||
uint64_t increment;
|
||||
};
|
||||
|
||||
SkipAheadCoefficients future_seed_coefficients(uint64_t n)
|
||||
{
|
||||
// The algorithm here to determine the parameters used to skip ahead is
|
||||
// described in F. Brown, "Random Number Generation with Arbitrary Stride,"
|
||||
// Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in
|
||||
// O(log2(N)) operations instead of O(N). Basically, it computes parameters G
|
||||
// and C which can then be used to find x_N = G*x_0 + C mod 2^M.
|
||||
|
||||
// Initialize constants
|
||||
uint64_t g {prn_mult};
|
||||
uint64_t c {prn_add};
|
||||
uint64_t g_new {1};
|
||||
uint64_t c_new {0};
|
||||
|
||||
while (n > 0) {
|
||||
// Check if the least significant bit is 1.
|
||||
if (n & 1) {
|
||||
g_new *= g;
|
||||
c_new = c_new * g + c;
|
||||
}
|
||||
c *= (g + 1);
|
||||
g *= g;
|
||||
|
||||
// Move bits right, dropping least significant bit.
|
||||
n >>= 1;
|
||||
}
|
||||
|
||||
return {g_new, c_new};
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
//==============================================================================
|
||||
// PRN
|
||||
//==============================================================================
|
||||
|
|
@ -69,9 +108,10 @@ uint64_t init_seed(int64_t id, int offset)
|
|||
|
||||
void init_particle_seeds(int64_t id, uint64_t* seeds)
|
||||
{
|
||||
auto [multiplier, increment] =
|
||||
future_seed_coefficients(static_cast<uint64_t>(id) * prn_stride);
|
||||
for (int i = 0; i < N_STREAMS; i++) {
|
||||
seeds[i] =
|
||||
future_seed(static_cast<uint64_t>(id) * prn_stride, master_seed + i);
|
||||
seeds[i] = multiplier * (master_seed + i) + increment;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -90,33 +130,9 @@ void advance_prn_seed(int64_t n, uint64_t* seed)
|
|||
|
||||
uint64_t future_seed(uint64_t n, uint64_t seed)
|
||||
{
|
||||
// The algorithm here to determine the parameters used to skip ahead is
|
||||
// described in F. Brown, "Random Number Generation with Arbitrary Stride,"
|
||||
// Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in
|
||||
// O(log2(N)) operations instead of O(N). Basically, it computes parameters G
|
||||
// and C which can then be used to find x_N = G*x_0 + C mod 2^M.
|
||||
|
||||
// Initialize constants
|
||||
uint64_t g {prn_mult};
|
||||
uint64_t c {prn_add};
|
||||
uint64_t g_new {1};
|
||||
uint64_t c_new {0};
|
||||
|
||||
while (n > 0) {
|
||||
// Check if the least significant bit is 1.
|
||||
if (n & 1) {
|
||||
g_new *= g;
|
||||
c_new = c_new * g + c;
|
||||
}
|
||||
c *= (g + 1);
|
||||
g *= g;
|
||||
|
||||
// Move bits right, dropping least significant bit.
|
||||
n >>= 1;
|
||||
}
|
||||
|
||||
// With G and C, we can now find the new seed.
|
||||
return g_new * seed + c_new;
|
||||
auto [multiplier, increment] = future_seed_coefficients(n);
|
||||
return multiplier * seed + increment;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
|
|||
|
|
@ -12,6 +12,7 @@
|
|||
#include "openmc/material.h"
|
||||
#include "openmc/message_passing.h"
|
||||
#include "openmc/nuclide.h"
|
||||
#include "openmc/openmp_interface.h"
|
||||
#include "openmc/output.h"
|
||||
#include "openmc/particle.h"
|
||||
#include "openmc/photon.h"
|
||||
|
|
@ -41,6 +42,7 @@
|
|||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <numeric>
|
||||
#include <string>
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -354,6 +356,94 @@ int64_t simulation_tracks_completed {0};
|
|||
|
||||
} // namespace simulation
|
||||
|
||||
namespace {
|
||||
|
||||
//! Collect thread-local secondary banks into the shared secondary bank in
|
||||
//! sorted order.
|
||||
//!
|
||||
//! \param thread_banks Secondary banks produced by each OpenMP thread
|
||||
void collect_sorted_history_secondary_banks(
|
||||
vector<vector<SourceSite>>& thread_banks)
|
||||
{
|
||||
// Count the total number of all secondary sites produced
|
||||
int64_t n_collected = 0;
|
||||
for (const auto& bank : thread_banks) {
|
||||
n_collected += bank.size();
|
||||
}
|
||||
|
||||
// Count the expected number of progeny from per-parent progeny counts
|
||||
int64_t n_progeny = 0;
|
||||
for (int64_t count : simulation::progeny_per_particle) {
|
||||
n_progeny += count;
|
||||
}
|
||||
|
||||
if (n_collected != n_progeny) {
|
||||
fatal_error("Mismatch detected between sum of all particle progeny and "
|
||||
"secondary bank size during collection.");
|
||||
}
|
||||
|
||||
// Convert per-parent progeny counts to offsets into the sorted bank
|
||||
std::exclusive_scan(simulation::progeny_per_particle.begin(),
|
||||
simulation::progeny_per_particle.end(),
|
||||
simulation::progeny_per_particle.begin(), 0);
|
||||
|
||||
// Allocate the shared bank once for the complete generation
|
||||
simulation::shared_secondary_bank_write.resize(0);
|
||||
simulation::shared_secondary_bank_write.extend_uninitialized(n_progeny);
|
||||
|
||||
// Place each secondary according to its parent and progeny identifiers
|
||||
for (const auto& bank : thread_banks) {
|
||||
for (const auto& site : bank) {
|
||||
if (site.parent_id < 0 ||
|
||||
site.parent_id >=
|
||||
static_cast<int64_t>(simulation::progeny_per_particle.size())) {
|
||||
fatal_error(fmt::format("Invalid parent_id {} for banked site "
|
||||
"(expected range [0, {})).",
|
||||
site.parent_id, simulation::progeny_per_particle.size()));
|
||||
}
|
||||
int64_t idx =
|
||||
simulation::progeny_per_particle[site.parent_id] + site.progeny_id;
|
||||
if (idx < 0 || idx >= n_progeny) {
|
||||
fatal_error("Mismatch detected between sum of all particle progeny and "
|
||||
"secondary bank size during collection.");
|
||||
}
|
||||
simulation::shared_secondary_bank_write[idx] = site;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//! Collect particle-local secondary banks into the shared secondary bank.
|
||||
//!
|
||||
//! \param n_particles Number of particles in the active event-based buffer
|
||||
void collect_event_secondary_banks(int64_t n_particles)
|
||||
{
|
||||
// Compute offsets for each particle's local secondary bank.
|
||||
vector<int64_t> offsets(n_particles);
|
||||
int64_t total = 0;
|
||||
for (int64_t i = 0; i < n_particles; ++i) {
|
||||
offsets[i] = total;
|
||||
total += simulation::particles[i].local_secondary_bank().size();
|
||||
}
|
||||
|
||||
// Extend the shared bank once for all collected secondaries
|
||||
int64_t bank_offset =
|
||||
simulation::shared_secondary_bank_write.extend_uninitialized(total);
|
||||
|
||||
// Copy each local bank into its assigned range and clear the local storage
|
||||
#pragma omp parallel for schedule(static)
|
||||
for (int64_t i = 0; i < n_particles; ++i) {
|
||||
auto& local_bank = simulation::particles[i].local_secondary_bank();
|
||||
if (!local_bank.empty()) {
|
||||
std::copy(local_bank.cbegin(), local_bank.cend(),
|
||||
simulation::shared_secondary_bank_write.data() + bank_offset +
|
||||
offsets[i]);
|
||||
local_bank.clear();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
//==============================================================================
|
||||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
|
@ -886,11 +976,14 @@ void transport_history_based_single_particle(Particle& p)
|
|||
|
||||
void transport_history_based()
|
||||
{
|
||||
#pragma omp parallel for schedule(runtime)
|
||||
for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
|
||||
#pragma omp parallel
|
||||
{
|
||||
Particle p;
|
||||
initialize_particle_track(p, i_work, false);
|
||||
transport_history_based_single_particle(p);
|
||||
#pragma omp for schedule(runtime)
|
||||
for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
|
||||
initialize_particle_track(p, i_work, false);
|
||||
transport_history_based_single_particle(p);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -918,30 +1011,27 @@ void transport_history_based_shared_secondary()
|
|||
std::fill(simulation::progeny_per_particle.begin(),
|
||||
simulation::progeny_per_particle.end(), 0);
|
||||
|
||||
vector<vector<SourceSite>> thread_banks(num_threads());
|
||||
|
||||
// Phase 1: Transport primary particles and deposit first generation of
|
||||
// secondaries in the shared secondary bank
|
||||
#pragma omp parallel
|
||||
{
|
||||
vector<SourceSite> thread_bank;
|
||||
auto& thread_bank = thread_banks[thread_num()];
|
||||
Particle p;
|
||||
|
||||
#pragma omp for schedule(runtime)
|
||||
for (int64_t i = 1; i <= simulation::work_per_rank; i++) {
|
||||
Particle p;
|
||||
initialize_particle_track(p, i, false);
|
||||
transport_history_based_single_particle(p);
|
||||
for (auto& site : p.local_secondary_bank()) {
|
||||
thread_bank.push_back(site);
|
||||
}
|
||||
}
|
||||
|
||||
// Drain thread-local bank into the shared secondary bank (once per thread)
|
||||
#pragma omp critical(SharedSecondaryBank)
|
||||
{
|
||||
for (auto& site : thread_bank) {
|
||||
simulation::shared_secondary_bank_write.thread_unsafe_append(site);
|
||||
}
|
||||
p.local_secondary_bank().clear();
|
||||
}
|
||||
}
|
||||
collect_sorted_history_secondary_banks(thread_banks);
|
||||
thread_banks.clear();
|
||||
|
||||
simulation::simulation_tracks_completed += settings::n_particles;
|
||||
|
||||
|
|
@ -951,10 +1041,6 @@ void transport_history_based_shared_secondary()
|
|||
int64_t alive_secondary = 1;
|
||||
while (alive_secondary) {
|
||||
|
||||
// Sort the shared secondary bank by parent ID then progeny ID to
|
||||
// ensure reproducibility.
|
||||
sort_bank(simulation::shared_secondary_bank_write, false);
|
||||
|
||||
// Synchronize the shared secondary bank amongst all MPI ranks, such
|
||||
// that each MPI rank has an approximately equal number of secondary
|
||||
// tracks. Also reports the total number of secondaries alive across
|
||||
|
|
@ -983,16 +1069,17 @@ void transport_history_based_shared_secondary()
|
|||
simulation::shared_secondary_bank_read.size());
|
||||
std::fill(simulation::progeny_per_particle.begin(),
|
||||
simulation::progeny_per_particle.end(), 0);
|
||||
thread_banks.resize(num_threads());
|
||||
|
||||
// Transport all secondary tracks from the shared secondary bank
|
||||
#pragma omp parallel
|
||||
{
|
||||
vector<SourceSite> thread_bank;
|
||||
auto& thread_bank = thread_banks[thread_num()];
|
||||
Particle p;
|
||||
|
||||
#pragma omp for schedule(runtime)
|
||||
for (int64_t i = 1; i <= simulation::shared_secondary_bank_read.size();
|
||||
i++) {
|
||||
Particle p;
|
||||
initialize_particle_track(p, i, true);
|
||||
SourceSite& site = simulation::shared_secondary_bank_read[i - 1];
|
||||
p.event_revive_from_secondary(site);
|
||||
|
|
@ -1000,18 +1087,14 @@ void transport_history_based_shared_secondary()
|
|||
for (auto& secondary_site : p.local_secondary_bank()) {
|
||||
thread_bank.push_back(secondary_site);
|
||||
}
|
||||
}
|
||||
|
||||
// Drain thread-local bank into the shared secondary bank (once per
|
||||
// thread)
|
||||
#pragma omp critical(SharedSecondaryBank)
|
||||
{
|
||||
for (auto& secondary_site : thread_bank) {
|
||||
simulation::shared_secondary_bank_write.thread_unsafe_append(
|
||||
secondary_site);
|
||||
}
|
||||
p.local_secondary_bank().clear();
|
||||
}
|
||||
} // End of transport loop over tracks in shared secondary bank
|
||||
simulation::shared_secondary_bank_write =
|
||||
std::move(simulation::shared_secondary_bank_read);
|
||||
simulation::shared_secondary_bank_read = SharedArray<SourceSite>();
|
||||
collect_sorted_history_secondary_banks(thread_banks);
|
||||
thread_banks.clear();
|
||||
n_generation_depth++;
|
||||
simulation::simulation_tracks_completed += alive_secondary;
|
||||
} // End of loop over secondary generations
|
||||
|
|
@ -1075,13 +1158,7 @@ void transport_event_based_shared_secondary()
|
|||
process_transport_events();
|
||||
process_death_events(n_particles);
|
||||
|
||||
// Collect secondaries from all particle buffers into shared bank
|
||||
for (int64_t i = 0; i < n_particles; i++) {
|
||||
for (auto& site : simulation::particles[i].local_secondary_bank()) {
|
||||
simulation::shared_secondary_bank_write.thread_unsafe_append(site);
|
||||
}
|
||||
simulation::particles[i].local_secondary_bank().clear();
|
||||
}
|
||||
collect_event_secondary_banks(n_particles);
|
||||
|
||||
remaining_work -= n_particles;
|
||||
source_offset += n_particles;
|
||||
|
|
@ -1145,13 +1222,7 @@ void transport_event_based_shared_secondary()
|
|||
process_transport_events();
|
||||
process_death_events(n_particles);
|
||||
|
||||
// Collect secondaries from all particle buffers into shared bank
|
||||
for (int64_t i = 0; i < n_particles; i++) {
|
||||
for (auto& site : simulation::particles[i].local_secondary_bank()) {
|
||||
simulation::shared_secondary_bank_write.thread_unsafe_append(site);
|
||||
}
|
||||
simulation::particles[i].local_secondary_bank().clear();
|
||||
}
|
||||
collect_event_secondary_banks(n_particles);
|
||||
|
||||
sec_remaining -= n_particles;
|
||||
sec_offset += n_particles;
|
||||
|
|
|
|||
474
src/source.cpp
474
src/source.cpp
|
|
@ -4,7 +4,9 @@
|
|||
#define HAS_DYNAMIC_LINKING
|
||||
#endif
|
||||
|
||||
#include <utility> // for move
|
||||
#include <algorithm> // for max
|
||||
#include <cmath> // for sin, cos, abs
|
||||
#include <utility> // for move
|
||||
|
||||
#ifdef HAS_DYNAMIC_LINKING
|
||||
#include <dlfcn.h> // for dlopen, dlsym, dlclose, dlerror
|
||||
|
|
@ -16,17 +18,20 @@
|
|||
#include "openmc/bank.h"
|
||||
#include "openmc/capi.h"
|
||||
#include "openmc/cell.h"
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/container_util.h"
|
||||
#include "openmc/error.h"
|
||||
#include "openmc/file_utils.h"
|
||||
#include "openmc/geometry.h"
|
||||
#include "openmc/hdf5_interface.h"
|
||||
#include "openmc/material.h"
|
||||
#include "openmc/math_functions.h"
|
||||
#include "openmc/mcpl_interface.h"
|
||||
#include "openmc/memory.h"
|
||||
#include "openmc/message_passing.h"
|
||||
#include "openmc/mgxs_interface.h"
|
||||
#include "openmc/nuclide.h"
|
||||
#include "openmc/random_dist.h"
|
||||
#include "openmc/random_lcg.h"
|
||||
#include "openmc/search.h"
|
||||
#include "openmc/settings.h"
|
||||
|
|
@ -100,6 +105,8 @@ unique_ptr<Source> Source::create(pugi::xml_node node)
|
|||
return make_unique<CompiledSourceWrapper>(node);
|
||||
} else if (source_type == "mesh") {
|
||||
return make_unique<MeshSource>(node);
|
||||
} else if (source_type == "tokamak") {
|
||||
return make_unique<TokamakSource>(node);
|
||||
} else {
|
||||
fatal_error(fmt::format("Invalid source type '{}' found.", source_type));
|
||||
}
|
||||
|
|
@ -673,6 +680,471 @@ SourceSite MeshSource::sample(uint64_t* seed) const
|
|||
return source(element)->sample_with_constraints(seed);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// TokamakSource implementation
|
||||
//==============================================================================
|
||||
|
||||
TokamakSource::TokamakSource(pugi::xml_node node) : Source(node)
|
||||
{
|
||||
// Read geometry parameters
|
||||
major_radius_ = std::stod(get_node_value(node, "major_radius"));
|
||||
minor_radius_ = std::stod(get_node_value(node, "minor_radius"));
|
||||
elongation_ = std::stod(get_node_value(node, "elongation"));
|
||||
triangularity_ = std::stod(get_node_value(node, "triangularity"));
|
||||
shafranov_shift_ = std::stod(get_node_value(node, "shafranov_shift"));
|
||||
|
||||
// Read optional vertical shift
|
||||
if (check_for_node(node, "vertical_shift")) {
|
||||
vertical_shift_ = std::stod(get_node_value(node, "vertical_shift"));
|
||||
} else {
|
||||
vertical_shift_ = 0.0;
|
||||
}
|
||||
|
||||
// Read optional toroidal angle bounds
|
||||
if (check_for_node(node, "phi_start")) {
|
||||
phi_start_ = std::stod(get_node_value(node, "phi_start"));
|
||||
} else {
|
||||
phi_start_ = 0.0;
|
||||
}
|
||||
if (check_for_node(node, "phi_extent")) {
|
||||
phi_extent_ = std::stod(get_node_value(node, "phi_extent"));
|
||||
} else {
|
||||
phi_extent_ = 2.0 * PI;
|
||||
}
|
||||
if (check_for_node(node, "n_alpha")) {
|
||||
n_alpha_ = std::stoi(get_node_value(node, "n_alpha"));
|
||||
} else {
|
||||
n_alpha_ = 101; // Default
|
||||
}
|
||||
|
||||
// Read emission profile
|
||||
r_over_a_ = get_node_array<double>(node, "r_over_a");
|
||||
emission_density_ = get_node_array<double>(node, "emission_density");
|
||||
|
||||
// Read energy distribution(s)
|
||||
for (auto energy_node : node.children("energy")) {
|
||||
energy_dists_.push_back(distribution_from_xml(energy_node));
|
||||
}
|
||||
|
||||
// Read optional time distribution; default to a delta distribution at t=0
|
||||
// for the same behavior as IndependentSource
|
||||
if (check_for_node(node, "time")) {
|
||||
time_ = distribution_from_xml(node.child("time"));
|
||||
} else {
|
||||
double T[] {0.0};
|
||||
double p[] {1.0};
|
||||
time_ = UPtrDist {new Discrete {T, p, 1}};
|
||||
}
|
||||
|
||||
// Validate inputs
|
||||
if (emission_density_.size() != r_over_a_.size()) {
|
||||
fatal_error("TokamakSource: emission_density and r_over_a must have the "
|
||||
"same length.");
|
||||
}
|
||||
if (r_over_a_.size() < 2) {
|
||||
fatal_error(
|
||||
"TokamakSource: At least 2 radial points are required for profiles.");
|
||||
}
|
||||
if (r_over_a_.front() != 0.0) {
|
||||
fatal_error("TokamakSource: r_over_a must start at 0.");
|
||||
}
|
||||
if (r_over_a_.back() != 1.0) {
|
||||
fatal_error("TokamakSource: r_over_a must end at 1.");
|
||||
}
|
||||
for (size_t i = 1; i < r_over_a_.size(); ++i) {
|
||||
if (r_over_a_[i] <= r_over_a_[i - 1]) {
|
||||
fatal_error("TokamakSource: r_over_a must be strictly increasing.");
|
||||
}
|
||||
}
|
||||
for (size_t i = 0; i < emission_density_.size(); ++i) {
|
||||
if (emission_density_[i] < 0.0) {
|
||||
fatal_error("TokamakSource: emission_density values cannot be negative.");
|
||||
}
|
||||
}
|
||||
if (major_radius_ <= 0.0) {
|
||||
fatal_error("TokamakSource: major_radius must be > 0.");
|
||||
}
|
||||
if (minor_radius_ <= 0.0) {
|
||||
fatal_error("TokamakSource: minor_radius must be > 0.");
|
||||
}
|
||||
if (minor_radius_ >= major_radius_) {
|
||||
fatal_error("TokamakSource: minor_radius must be less than major_radius.");
|
||||
}
|
||||
if (elongation_ <= 0.0) {
|
||||
fatal_error("TokamakSource: elongation must be > 0.");
|
||||
}
|
||||
if (triangularity_ < -1.0 || triangularity_ > 1.0) {
|
||||
fatal_error("TokamakSource: triangularity must be in the range [-1, 1].");
|
||||
}
|
||||
if (shafranov_shift_ < 0.0) {
|
||||
fatal_error("TokamakSource: shafranov_shift must be >= 0.");
|
||||
}
|
||||
if (shafranov_shift_ >= 0.5 * minor_radius_) {
|
||||
fatal_error("TokamakSource: shafranov_shift must be less than half the "
|
||||
"minor radius.");
|
||||
}
|
||||
if (phi_extent_ <= 0.0 || phi_extent_ > 2.0 * PI) {
|
||||
fatal_error("TokamakSource: phi_extent must be > 0 and <= 2*pi.");
|
||||
}
|
||||
if (n_alpha_ <= 2) {
|
||||
fatal_error("TokamakSource: n_alpha must be > 2.");
|
||||
}
|
||||
if (n_alpha_ < 51) {
|
||||
warning("TokamakSource: n_alpha values below 51 may introduce noticeable "
|
||||
"discretization bias in source sampling.");
|
||||
}
|
||||
if (energy_dists_.empty()) {
|
||||
fatal_error("TokamakSource: At least one energy distribution is required.");
|
||||
}
|
||||
if (energy_dists_.size() != 1 && energy_dists_.size() != r_over_a_.size()) {
|
||||
fatal_error("TokamakSource: energy distributions must be either 1 (for all "
|
||||
"r) or match the number of r_over_a points.");
|
||||
}
|
||||
|
||||
// Compute normalized geometry parameters
|
||||
epsilon_ = minor_radius_ / major_radius_;
|
||||
delta_tilde_ = shafranov_shift_ / minor_radius_;
|
||||
|
||||
// Initialize isotropic angular distribution
|
||||
angle_ = UPtrAngle {new Isotropic()};
|
||||
|
||||
precompute_sampling_distributions();
|
||||
}
|
||||
|
||||
void TokamakSource::precompute_sampling_distributions()
|
||||
{
|
||||
// Use precomputed normalized geometry parameters
|
||||
double eps = epsilon_; // Inverse aspect ratio (a/R0)
|
||||
double Dt = delta_tilde_; // Normalized Shafranov shift (Delta/a)
|
||||
double delta = triangularity_;
|
||||
|
||||
//==========================================================================
|
||||
// RADIAL CDF (computed first since it's simpler and sampled first)
|
||||
//==========================================================================
|
||||
// The marginal radial PDF is obtained by analytically integrating the joint
|
||||
// distribution f(r_tilde, alpha) over alpha. The result is:
|
||||
//
|
||||
// p(r_tilde) ~ S(r_tilde) * [(1 + eps*Dt)*r_tilde
|
||||
// - (3/8)*c1*eps*r_tilde^2
|
||||
// - 2*eps*Dt*r_tilde^3]
|
||||
//
|
||||
// where the Bessel function coefficients are:
|
||||
// c0 = J_0(delta) + J_2(delta)
|
||||
// c1 = (J_1(2*delta) + J_3(2*delta)) / c0
|
||||
//
|
||||
// For delta -> 0, c0 -> 1 and c1 -> 0, giving the circular cross-section
|
||||
// limit.
|
||||
|
||||
// Compute Bessel function coefficients. openmc::cyl_bessel_j handles
|
||||
// negative arguments (negative triangularity) via the parity relation
|
||||
// J_n(-x) = (-1)^n * J_n(x).
|
||||
double J0_d = cyl_bessel_j(0, delta);
|
||||
double J2_d = cyl_bessel_j(2, delta);
|
||||
double J1_2d = cyl_bessel_j(1, 2.0 * delta);
|
||||
double J3_2d = cyl_bessel_j(3, 2.0 * delta);
|
||||
double c0 = J0_d + J2_d;
|
||||
double c1 = (J1_2d + J3_2d) / c0;
|
||||
|
||||
// Coefficients for the radial polynomial: A*r - B*r^2 - C*r^3
|
||||
radial_poly_a_ = 1.0 + eps * Dt;
|
||||
radial_poly_b_ = 0.375 * c1 * eps; // 3/8 * c1 * eps
|
||||
radial_poly_c_ = 2.0 * eps * Dt;
|
||||
|
||||
// Build a refined radial grid that retains the user-specified grid points.
|
||||
// The emission density is interpreted as linear-linear between those points.
|
||||
constexpr int MIN_SUBINTERVALS = 8;
|
||||
constexpr double MAX_GRID_SPACING = 1.0e-3;
|
||||
vector<double> radial_grid {r_over_a_.front()};
|
||||
vector<double> radial_emission {emission_density_.front()};
|
||||
for (size_t i = 1; i < r_over_a_.size(); ++i) {
|
||||
double r_lo = r_over_a_[i - 1];
|
||||
double r_hi = r_over_a_[i];
|
||||
double s_lo = emission_density_[i - 1];
|
||||
double s_hi = emission_density_[i];
|
||||
int n_subintervals = std::max(MIN_SUBINTERVALS,
|
||||
static_cast<int>(std::ceil((r_hi - r_lo) / MAX_GRID_SPACING)));
|
||||
for (int j = 1; j <= n_subintervals; ++j) {
|
||||
double t = static_cast<double>(j) / n_subintervals;
|
||||
radial_grid.push_back(r_lo + t * (r_hi - r_lo));
|
||||
radial_emission.push_back(s_lo + t * (s_hi - s_lo));
|
||||
}
|
||||
}
|
||||
|
||||
vector<double> radial_pdf(radial_grid.size());
|
||||
for (size_t i = 0; i < radial_grid.size(); ++i) {
|
||||
double r = radial_grid[i];
|
||||
// p(r) ~ S(r) * [A*r - B*r^2 - C*r^3]
|
||||
double geometric_factor =
|
||||
radial_poly_a_ * r - radial_poly_b_ * r * r - radial_poly_c_ * r * r * r;
|
||||
radial_pdf[i] = radial_emission[i] * std::max(0.0, geometric_factor);
|
||||
}
|
||||
|
||||
// Check that the refined profile contains positive probability mass before
|
||||
// constructing the normalized tabular distribution.
|
||||
double total = 0.0;
|
||||
for (size_t i = 1; i < radial_grid.size(); ++i) {
|
||||
total += 0.5 * (radial_pdf[i - 1] + radial_pdf[i]) *
|
||||
(radial_grid[i] - radial_grid[i - 1]);
|
||||
}
|
||||
if (total <= 0.0) {
|
||||
fatal_error(
|
||||
"TokamakSource: Integrated emission density is zero or negative. "
|
||||
"Check emission_density profile.");
|
||||
}
|
||||
radial_dist_ = make_unique<Tabular>(radial_grid.data(), radial_pdf.data(),
|
||||
radial_grid.size(), Interpolation::lin_lin);
|
||||
|
||||
//==========================================================================
|
||||
// POLOIDAL CDFs (for conditional sampling of alpha given r)
|
||||
//==========================================================================
|
||||
// The conditional distribution P(alpha | r) is a mixture:
|
||||
// P(alpha | r) ~ sum_k w_k(r) * I_hat_k * p_k(alpha)
|
||||
// where:
|
||||
// - w_k(r) are the "dynamic" Bernstein weight functions (depend on r)
|
||||
// - I_hat_k are the "static" normalized integrals (precomputed constants)
|
||||
// - p_k(alpha) are the normalized basis distributions (precomputed CDFs)
|
||||
//
|
||||
// The static weights I_hat_k = I_k / (2*pi*c0) are:
|
||||
// I_hat_0 = 1 + eps*Dt
|
||||
// I_hat_1 = 1 + eps*Dt - (3/16)*c1*eps
|
||||
// I_hat_2 = 1 - (3/8)*c1*eps
|
||||
// I_hat_3 = 1 + eps*Dt
|
||||
// I_hat_4 = 1 + (1/2)*eps*Dt - (3/16)*c1*eps
|
||||
// I_hat_5 = 1 - eps*Dt - (3/8)*c1*eps
|
||||
|
||||
// Compute static weights analytically
|
||||
poloidal_integrals_[0] = 1.0 + eps * Dt;
|
||||
poloidal_integrals_[1] = 1.0 + eps * Dt - 0.1875 * c1 * eps; // 3/16 = 0.1875
|
||||
poloidal_integrals_[2] = 1.0 - 0.375 * c1 * eps; // 3/8 = 0.375
|
||||
poloidal_integrals_[3] = 1.0 + eps * Dt;
|
||||
poloidal_integrals_[4] = 1.0 + 0.5 * eps * Dt - 0.1875 * c1 * eps;
|
||||
poloidal_integrals_[5] = 1.0 - eps * Dt - 0.375 * c1 * eps;
|
||||
|
||||
// Build the alpha grid on [0, pi] (half domain due to up-down symmetry)
|
||||
int n_alpha = n_alpha_;
|
||||
vector<double> alpha_grid(n_alpha);
|
||||
double dalpha = PI / (n_alpha - 1);
|
||||
for (int i = 0; i < n_alpha; ++i) {
|
||||
alpha_grid[i] = i * dalpha;
|
||||
}
|
||||
|
||||
// Compute basis function values g_k(alpha) for tabular distributions
|
||||
// Using Bernstein form:
|
||||
// R_tilde = b0*(1-r)^2 + 2*b1*r*(1-r) + b2*r^2
|
||||
// J_tilde = b3*(1-r) + b4*r
|
||||
// with:
|
||||
// b0(alpha) = 1 + eps*Dt
|
||||
// b1(alpha) = b0 + (eps/2)*cos(psi), psi = alpha + delta*sin(alpha)
|
||||
// b2(alpha) = 1 + eps*cos(psi)
|
||||
// b3(alpha) = cos(delta*sin(alpha))
|
||||
// + (delta/4)*(cos(alpha - delta*sin(alpha))
|
||||
// - cos(3*alpha + delta*sin(alpha)))
|
||||
// b4(alpha) = b3(alpha) - 2*Dt*cos(alpha)
|
||||
|
||||
array<vector<double>, N_POLOIDAL_BASIS> basis;
|
||||
for (int k = 0; k < N_POLOIDAL_BASIS; ++k) {
|
||||
basis[k].resize(n_alpha);
|
||||
}
|
||||
|
||||
for (int i = 0; i < n_alpha; ++i) {
|
||||
double alpha = alpha_grid[i];
|
||||
double sin_alpha = std::sin(alpha);
|
||||
double cos_alpha = std::cos(alpha);
|
||||
double delta_sin_alpha = delta * sin_alpha;
|
||||
double psi = alpha + delta_sin_alpha;
|
||||
double cos_psi = std::cos(psi);
|
||||
|
||||
// Bernstein coefficients b0-b4
|
||||
double b0 = 1.0 + eps * Dt;
|
||||
double b1 = b0 + 0.5 * eps * cos_psi;
|
||||
double b2 = 1.0 + eps * cos_psi;
|
||||
double b3 =
|
||||
std::cos(delta_sin_alpha) + 0.25 * delta *
|
||||
(std::cos(alpha - delta_sin_alpha) -
|
||||
std::cos(3.0 * alpha + delta_sin_alpha));
|
||||
double b4 = b3 - 2.0 * Dt * cos_alpha;
|
||||
|
||||
// 6 basis functions g_k(alpha) = b_i * b_j
|
||||
basis[0][i] = b0 * b3; // w0 = (1-r)^3
|
||||
basis[1][i] = b1 * b3; // w1 = 2*r*(1-r)^2
|
||||
basis[2][i] = b2 * b3; // w2 = r^2*(1-r)
|
||||
basis[3][i] = b0 * b4; // w3 = r*(1-r)^2
|
||||
basis[4][i] = b1 * b4; // w4 = 2*r^2*(1-r)
|
||||
basis[5][i] = b2 * b4; // w5 = r^3
|
||||
}
|
||||
|
||||
// Build a linear-linear distribution for each basis function p_k(alpha)
|
||||
for (int k = 0; k < N_POLOIDAL_BASIS; ++k) {
|
||||
poloidal_dists_[k] = make_unique<Tabular>(
|
||||
alpha_grid.data(), basis[k].data(), n_alpha, Interpolation::lin_lin);
|
||||
}
|
||||
}
|
||||
|
||||
double TokamakSource::sample_r_over_a(uint64_t* seed) const
|
||||
{
|
||||
return radial_dist_->sample(seed).first;
|
||||
}
|
||||
|
||||
double TokamakSource::mixture_weight(int k, double r) const
|
||||
{
|
||||
double s = 1.0 - r;
|
||||
switch (k) {
|
||||
case 0:
|
||||
return s * s * s * poloidal_integrals_[0];
|
||||
case 1:
|
||||
return 2.0 * r * s * s * poloidal_integrals_[1];
|
||||
case 2:
|
||||
return r * r * s * poloidal_integrals_[2];
|
||||
case 3:
|
||||
return r * s * s * poloidal_integrals_[3];
|
||||
case 4:
|
||||
return 2.0 * r * r * s * poloidal_integrals_[4];
|
||||
case 5:
|
||||
return r * r * r * poloidal_integrals_[5];
|
||||
default:
|
||||
UNREACHABLE();
|
||||
}
|
||||
}
|
||||
|
||||
double TokamakSource::sample_poloidal_angle(double r_norm, uint64_t* seed) const
|
||||
{
|
||||
// Sample from the conditional distribution P(alpha | r_tilde) using
|
||||
// mixture sampling with 6 precomputed basis distributions.
|
||||
//
|
||||
// The conditional is: P(alpha | r) ~ sum_k w_k(r) * I_hat_k * p_k(alpha)
|
||||
// where:
|
||||
// - w_k(r) are the "dynamic" Bernstein weight functions
|
||||
// - I_hat_k are the "static" normalized integrals (precomputed in
|
||||
// poloidal_integrals_)
|
||||
// - p_k(alpha) are the normalized, precomputed basis distributions
|
||||
//
|
||||
// The normalization sum_k w_k(r) * I_hat_k equals the radial geometric
|
||||
// polynomial evaluated at r, which is known analytically.
|
||||
//
|
||||
// Algorithm:
|
||||
// 1. Compute total from analytical normalization
|
||||
// 2. Lazily evaluate mixture weights with early exit to select component k
|
||||
// 3. Sample alpha from the selected basis distribution
|
||||
|
||||
// Analytical normalization: sum_k w_k(r) * I_hat_k
|
||||
double total =
|
||||
radial_poly_a_ - radial_poly_b_ * r_norm - radial_poly_c_ * r_norm * r_norm;
|
||||
double xi = prn(seed) * total;
|
||||
|
||||
// Sample component via lazy evaluation with early exit
|
||||
// Order optimized for peaked emission profiles: 0, 1, 4, 5, 3, 2
|
||||
constexpr int order[] = {0, 1, 4, 5, 3, 2};
|
||||
double cumsum = 0.0;
|
||||
int component = order[N_POLOIDAL_BASIS - 1];
|
||||
for (int i = 0; i < N_POLOIDAL_BASIS; ++i) {
|
||||
cumsum += mixture_weight(order[i], r_norm);
|
||||
if (xi < cumsum) {
|
||||
component = order[i];
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
// Sample alpha from [0, pi]
|
||||
double alpha = poloidal_dists_[component]->sample(seed).first;
|
||||
|
||||
// Exploit up-down symmetry: randomly flip to [pi, 2*pi] with 50% probability
|
||||
// This is equivalent to flipping the sign of Z in the final position
|
||||
if (prn(seed) >= 0.5) {
|
||||
alpha = 2.0 * PI - alpha;
|
||||
}
|
||||
return alpha;
|
||||
}
|
||||
|
||||
std::pair<double, double> TokamakSource::sample_energy(
|
||||
double r_norm, uint64_t* seed) const
|
||||
{
|
||||
if (energy_dists_.size() == 1) {
|
||||
// Single distribution for all r
|
||||
return energy_dists_[0]->sample(seed);
|
||||
}
|
||||
|
||||
// Multiple distributions: stochastic selection between bracketing r points
|
||||
// Find the interval containing r_norm
|
||||
size_t i = lower_bound_index(r_over_a_.begin(), r_over_a_.end(), r_norm);
|
||||
|
||||
// Handle boundary cases
|
||||
if (i >= energy_dists_.size() - 1) {
|
||||
return energy_dists_.back()->sample(seed);
|
||||
}
|
||||
|
||||
// Stochastic interpolation: randomly select one of the two bracketing
|
||||
// distributions based on distance to each
|
||||
double t = (r_norm - r_over_a_[i]) / (r_over_a_[i + 1] - r_over_a_[i]);
|
||||
size_t idx = (prn(seed) < t) ? i + 1 : i;
|
||||
return energy_dists_[idx]->sample(seed);
|
||||
}
|
||||
|
||||
Position TokamakSource::flux_to_cartesian(
|
||||
double r, double alpha, double phi) const
|
||||
{
|
||||
// Flux surface parameterization:
|
||||
// R = R0 + r*cos(alpha + delta*sin(alpha)) + Delta*(1 - (r/a)^2)
|
||||
// Z = kappa * r * sin(alpha)
|
||||
// x = R * cos(phi)
|
||||
// y = R * sin(phi)
|
||||
// z = Z
|
||||
|
||||
double psi = alpha + triangularity_ * std::sin(alpha);
|
||||
double r_over_a_sq = (r * r) / (minor_radius_ * minor_radius_);
|
||||
|
||||
double R =
|
||||
major_radius_ + r * std::cos(psi) + shafranov_shift_ * (1.0 - r_over_a_sq);
|
||||
double Z = elongation_ * r * std::sin(alpha);
|
||||
|
||||
double x = R * std::cos(phi);
|
||||
double y = R * std::sin(phi);
|
||||
double z = Z;
|
||||
|
||||
return {x, y, z};
|
||||
}
|
||||
|
||||
SourceSite TokamakSource::sample(uint64_t* seed) const
|
||||
{
|
||||
SourceSite site;
|
||||
site.particle = ParticleType::neutron();
|
||||
site.wgt = 1.0;
|
||||
site.delayed_group = 0;
|
||||
|
||||
// 1. Sample r/a from radial CDF
|
||||
double r_norm = sample_r_over_a(seed);
|
||||
double r = r_norm * minor_radius_;
|
||||
|
||||
// 2. Sample poloidal angle from conditional distribution P(alpha|r)
|
||||
double alpha = sample_poloidal_angle(r_norm, seed);
|
||||
|
||||
// 3. Sample toroidal angle uniformly in [phi_start, phi_start + phi_extent]
|
||||
double phi = phi_start_ + phi_extent_ * prn(seed);
|
||||
|
||||
// 4. Convert to Cartesian coordinates
|
||||
site.r = flux_to_cartesian(r, alpha, phi);
|
||||
|
||||
// 4a. Apply vertical shift if non-zero
|
||||
if (vertical_shift_ != 0.0) {
|
||||
site.r.z += vertical_shift_;
|
||||
}
|
||||
|
||||
// 5. Sample isotropic direction
|
||||
site.u = angle_->sample(seed).first;
|
||||
|
||||
// 6. Sample energy from distribution(s), applying the importance weight so
|
||||
// that biased distributions are handled correctly
|
||||
auto [E, E_wgt] = sample_energy(r_norm, seed);
|
||||
site.E = E;
|
||||
|
||||
// 7. Sample particle creation time
|
||||
auto [time, time_wgt] = time_->sample(seed);
|
||||
site.time = time;
|
||||
|
||||
site.wgt *= E_wgt * time_wgt;
|
||||
|
||||
return site;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
|
|
|||
|
|
@ -2730,6 +2730,9 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
|
|||
int orig_cell = p.coord(0).cell();
|
||||
double orig_E_last = p.E_last();
|
||||
|
||||
// Set particle in top level
|
||||
p.n_coord() = 1;
|
||||
|
||||
for (auto i_tally : tallies) {
|
||||
auto& tally {*model::tallies[i_tally]};
|
||||
|
||||
|
|
@ -2740,7 +2743,6 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
|
|||
|
||||
for (auto cell_id : cells) {
|
||||
// Temporarily change cell of particle
|
||||
p.n_coord() = 1;
|
||||
p.coord(0).cell() = cell_id;
|
||||
|
||||
// Determine index of cell in model::pulse_height_cells
|
||||
|
|
@ -2756,20 +2758,20 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
|
|||
// we skip the assume_separate break below.
|
||||
auto filter_iter = FilterBinIter(tally, p);
|
||||
auto end = FilterBinIter(tally, true, &p.filter_matches());
|
||||
if (filter_iter == end)
|
||||
continue;
|
||||
if (filter_iter != end) {
|
||||
|
||||
// Loop over filter bins.
|
||||
for (; filter_iter != end; ++filter_iter) {
|
||||
auto filter_index = filter_iter.index_;
|
||||
auto filter_weight = filter_iter.weight_;
|
||||
// Loop over filter bins.
|
||||
for (; filter_iter != end; ++filter_iter) {
|
||||
auto filter_index = filter_iter.index_;
|
||||
auto filter_weight = filter_iter.weight_;
|
||||
|
||||
// Loop over scores.
|
||||
for (auto score_index = 0; score_index < tally.scores_.size();
|
||||
++score_index) {
|
||||
// Loop over scores.
|
||||
for (auto score_index = 0; score_index < tally.scores_.size();
|
||||
++score_index) {
|
||||
#pragma omp atomic
|
||||
tally.results_(filter_index, score_index, TallyResult::VALUE) +=
|
||||
filter_weight;
|
||||
tally.results_(filter_index, score_index, TallyResult::VALUE) +=
|
||||
filter_weight;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -2777,10 +2779,10 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
|
|||
for (auto& match : p.filter_matches())
|
||||
match.bins_present_ = false;
|
||||
}
|
||||
// Restore cell/energy
|
||||
p.n_coord() = orig_n_coord;
|
||||
p.coord(0).cell() = orig_cell;
|
||||
p.E_last() = orig_E_last;
|
||||
}
|
||||
// Restore cell/energy
|
||||
p.n_coord() = orig_n_coord;
|
||||
p.coord(0).cell() = orig_cell;
|
||||
p.E_last() = orig_E_last;
|
||||
}
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -837,7 +837,8 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
|
|||
ratio_));
|
||||
if (ratio_ <= 1.0)
|
||||
fatal_error(fmt::format("Invalid weight window ratio '{}' (<= 1.0) "
|
||||
"specified for weight window generation"));
|
||||
"specified for weight window generation",
|
||||
ratio_));
|
||||
|
||||
// create a matching weight windows object
|
||||
auto wws = WeightWindows::create();
|
||||
|
|
|
|||
|
|
@ -82,6 +82,31 @@ TEST_CASE("Test alias sampling method for pugixml constructor")
|
|||
}
|
||||
}
|
||||
|
||||
TEST_CASE("Test sampling a large linear-linear tabular distribution")
|
||||
{
|
||||
constexpr int n_points = 10001;
|
||||
constexpr int n_samples = 200000;
|
||||
openmc::vector<double> x(n_points);
|
||||
openmc::vector<double> p(n_points);
|
||||
for (int i = 0; i < n_points; ++i) {
|
||||
x[i] = static_cast<double>(i) / (n_points - 1);
|
||||
p[i] = 2.0 * x[i];
|
||||
}
|
||||
|
||||
openmc::Tabular dist(
|
||||
x.data(), p.data(), n_points, openmc::Interpolation::lin_lin);
|
||||
uint64_t seed = openmc::init_seed(0, 0);
|
||||
|
||||
double mean = 0.0;
|
||||
for (int i = 0; i < n_samples; ++i) {
|
||||
mean += dist.sample(&seed).first;
|
||||
}
|
||||
mean /= n_samples;
|
||||
|
||||
// The normalized PDF is 2x on [0, 1], which has a mean of 2/3.
|
||||
REQUIRE_THAT(mean, Catch::Matchers::WithinAbs(2.0 / 3.0, 0.003));
|
||||
}
|
||||
|
||||
TEST_CASE("Test construction of SpatialBox with parameters")
|
||||
{
|
||||
openmc::Position ll {-1, -2, -3};
|
||||
|
|
|
|||
|
|
@ -46,6 +46,21 @@ TEST_CASE("Test t_percentile")
|
|||
}
|
||||
}
|
||||
|
||||
TEST_CASE("Test cylindrical Bessel functions")
|
||||
{
|
||||
constexpr double x = 2.0;
|
||||
constexpr double expected[] {0.22389077914123567, 0.57672480775687339,
|
||||
0.35283402861563773, 0.12894324947440205};
|
||||
|
||||
for (int n = 0; n < 4; ++n) {
|
||||
REQUIRE_THAT(openmc::cyl_bessel_j(n, x),
|
||||
Catch::Matchers::WithinRel(expected[n], 1.0e-14));
|
||||
REQUIRE_THAT(openmc::cyl_bessel_j(n, -x),
|
||||
Catch::Matchers::WithinRel(
|
||||
(n % 2 == 0 ? 1.0 : -1.0) * expected[n], 1.0e-14));
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("Test calc_pn")
|
||||
{
|
||||
// The reference solutions come from scipy.special.eval_legendre
|
||||
|
|
@ -355,3 +370,24 @@ TEST_CASE("Test broaden_wmp_polynomials")
|
|||
REQUIRE_THAT(ref_val, Catch::Matchers::Approx(test_val));
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("Test isclose")
|
||||
{
|
||||
using openmc::isclose;
|
||||
|
||||
// Identical values are always close, regardless of tolerances.
|
||||
REQUIRE(isclose(1.0, 1.0, 0.0, 0.0));
|
||||
REQUIRE(isclose(0.0, 0.0, 0.0, 0.0));
|
||||
|
||||
// Absolute tolerance governs comparisons near zero.
|
||||
REQUIRE(isclose(0.0, 1e-15, 0.0, 1e-14));
|
||||
REQUIRE_FALSE(isclose(0.0, 1e-13, 0.0, 1e-14));
|
||||
|
||||
// Relative tolerance scales with the magnitude of the operands.
|
||||
REQUIRE(isclose(1.0e12, 1.0e12 + 1.0, 1e-12, 0.0));
|
||||
REQUIRE_FALSE(isclose(1.0e12, 1.0e12 + 10.0, 1e-12, 0.0));
|
||||
|
||||
// The looser of the two tolerances wins.
|
||||
REQUIRE(isclose(1.0, 1.0 + 1e-13, 0.0, 1e-12));
|
||||
REQUIRE(isclose(1.0e6, 1.0e6 + 1e-4, 1e-9, 0.0));
|
||||
}
|
||||
|
|
|
|||
0
tests/regression_tests/lattice_large_pitch/__init__.py
Normal file
0
tests/regression_tests/lattice_large_pitch/__init__.py
Normal file
51
tests/regression_tests/lattice_large_pitch/inputs_true.dat
Normal file
51
tests/regression_tests/lattice_large_pitch/inputs_true.dat
Normal file
|
|
@ -0,0 +1,51 @@
|
|||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<model>
|
||||
<materials>
|
||||
<material id="1">
|
||||
<density value="0.001" units="g/cm3"/>
|
||||
<nuclide name="N14" ao="1.0"/>
|
||||
</material>
|
||||
<material id="2">
|
||||
<density value="7.0" units="g/cm3"/>
|
||||
<nuclide name="Fe56" ao="1.0"/>
|
||||
</material>
|
||||
</materials>
|
||||
<geometry>
|
||||
<cell id="1" material="2" universe="1"/>
|
||||
<cell id="2" material="1" universe="2"/>
|
||||
<cell id="3" fill="3" region="1 -2 3 -4" universe="4"/>
|
||||
<lattice id="3">
|
||||
<pitch>100.0 100.0</pitch>
|
||||
<outer>2</outer>
|
||||
<dimension>3 3</dimension>
|
||||
<lower_left>-150.0 -150.0</lower_left>
|
||||
<universes>
|
||||
1 2 1
|
||||
2 1 2
|
||||
1 2 1 </universes>
|
||||
</lattice>
|
||||
<surface id="1" name="minimum x" type="x-plane" boundary="vacuum" coeffs="-150.0"/>
|
||||
<surface id="2" name="maximum x" type="x-plane" boundary="vacuum" coeffs="150.0"/>
|
||||
<surface id="3" name="minimum y" type="y-plane" boundary="vacuum" coeffs="-150.0"/>
|
||||
<surface id="4" name="maximum y" type="y-plane" boundary="vacuum" coeffs="150.0"/>
|
||||
</geometry>
|
||||
<settings>
|
||||
<run_mode>fixed source</run_mode>
|
||||
<particles>1000</particles>
|
||||
<batches>10</batches>
|
||||
</settings>
|
||||
<tallies>
|
||||
<mesh id="1">
|
||||
<dimension>6 6</dimension>
|
||||
<lower_left>-150.0 -150.0</lower_left>
|
||||
<upper_right>150.0 150.0</upper_right>
|
||||
</mesh>
|
||||
<filter id="1" type="mesh">
|
||||
<bins>1</bins>
|
||||
</filter>
|
||||
<tally id="1">
|
||||
<filters>1</filters>
|
||||
<scores>flux</scores>
|
||||
</tally>
|
||||
</tallies>
|
||||
</model>
|
||||
73
tests/regression_tests/lattice_large_pitch/results_true.dat
Normal file
73
tests/regression_tests/lattice_large_pitch/results_true.dat
Normal file
|
|
@ -0,0 +1,73 @@
|
|||
tally 1:
|
||||
1.749265E+01
|
||||
3.275925E+01
|
||||
4.159430E+01
|
||||
1.799568E+02
|
||||
7.027454E+01
|
||||
4.967568E+02
|
||||
6.789716E+01
|
||||
4.652258E+02
|
||||
4.272935E+01
|
||||
1.861129E+02
|
||||
2.073342E+01
|
||||
4.374279E+01
|
||||
4.028297E+01
|
||||
1.641642E+02
|
||||
6.734263E+01
|
||||
4.627171E+02
|
||||
1.037657E+02
|
||||
1.078682E+03
|
||||
1.108394E+02
|
||||
1.240927E+03
|
||||
7.236770E+01
|
||||
5.302034E+02
|
||||
4.356361E+01
|
||||
1.921990E+02
|
||||
6.731536E+01
|
||||
4.615977E+02
|
||||
9.850684E+01
|
||||
9.810272E+02
|
||||
5.164863E+02
|
||||
2.670336E+04
|
||||
5.097770E+02
|
||||
2.604770E+04
|
||||
1.064847E+02
|
||||
1.137536E+03
|
||||
7.052986E+01
|
||||
5.003220E+02
|
||||
7.065046E+01
|
||||
5.036723E+02
|
||||
1.044890E+02
|
||||
1.103250E+03
|
||||
5.121544E+02
|
||||
2.632389E+04
|
||||
5.065331E+02
|
||||
2.570752E+04
|
||||
1.044794E+02
|
||||
1.101249E+03
|
||||
6.569142E+01
|
||||
4.361120E+02
|
||||
4.739812E+01
|
||||
2.313289E+02
|
||||
7.402284E+01
|
||||
5.591689E+02
|
||||
1.066905E+02
|
||||
1.149521E+03
|
||||
1.038665E+02
|
||||
1.086813E+03
|
||||
7.160008E+01
|
||||
5.217744E+02
|
||||
4.219705E+01
|
||||
1.816100E+02
|
||||
2.089838E+01
|
||||
4.464899E+01
|
||||
4.154271E+01
|
||||
1.745866E+02
|
||||
7.081253E+01
|
||||
5.049997E+02
|
||||
6.673273E+01
|
||||
4.509038E+02
|
||||
4.184624E+01
|
||||
1.771546E+02
|
||||
1.853898E+01
|
||||
3.538608E+01
|
||||
70
tests/regression_tests/lattice_large_pitch/test.py
Normal file
70
tests/regression_tests/lattice_large_pitch/test.py
Normal file
|
|
@ -0,0 +1,70 @@
|
|||
"""Regression test for rectangular lattices with large pitch values.
|
||||
|
||||
Large pitches used to trigger a segmentation fault in ``RectLattice::distance``
|
||||
because the boundary-crossing check compared a reconstructed crossing position
|
||||
against an absolute tolerance that did not scale with the geometry (see #3852).
|
||||
This test transports particles across a lattice with a large pitch to ensure the
|
||||
crossing logic remains robust.
|
||||
|
||||
"""
|
||||
|
||||
import openmc
|
||||
import pytest
|
||||
|
||||
from tests.testing_harness import PyAPITestHarness
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def model():
|
||||
model = openmc.Model()
|
||||
|
||||
# Large pitch that previously caused floating-point cancellation
|
||||
pitch = 100.0
|
||||
n = 3
|
||||
|
||||
air = openmc.Material()
|
||||
air.set_density('g/cm3', 0.001)
|
||||
air.add_nuclide('N14', 1.0)
|
||||
metal = openmc.Material()
|
||||
metal.set_density('g/cm3', 7.0)
|
||||
metal.add_nuclide('Fe56', 1.0)
|
||||
|
||||
metal_cell = openmc.Cell(fill=metal)
|
||||
metal_uni = openmc.Universe(cells=[metal_cell])
|
||||
air_cell = openmc.Cell(fill=air)
|
||||
air_uni = openmc.Universe(cells=[air_cell])
|
||||
|
||||
lattice = openmc.RectLattice()
|
||||
lattice.lower_left = (-pitch*n/2, -pitch*n/2)
|
||||
lattice.pitch = (pitch, pitch)
|
||||
lattice.outer = air_uni
|
||||
lattice.universes = [
|
||||
[metal_uni, air_uni, metal_uni],
|
||||
[air_uni, metal_uni, air_uni],
|
||||
[metal_uni, air_uni, metal_uni],
|
||||
]
|
||||
|
||||
box = openmc.model.RectangularPrism(pitch*n, pitch*n, boundary_type='vacuum')
|
||||
root_cell = openmc.Cell(region=-box, fill=lattice)
|
||||
model.geometry = openmc.Geometry([root_cell])
|
||||
|
||||
model.settings.run_mode = 'fixed source'
|
||||
model.settings.batches = 10
|
||||
model.settings.particles = 1000
|
||||
|
||||
mesh = openmc.RegularMesh()
|
||||
mesh.dimension = (6, 6)
|
||||
mesh.lower_left = (-pitch*n/2, -pitch*n/2)
|
||||
mesh.upper_right = (pitch*n/2, pitch*n/2)
|
||||
mesh_filter = openmc.MeshFilter(mesh)
|
||||
tally = openmc.Tally(tally_id=1)
|
||||
tally.filters = [mesh_filter]
|
||||
tally.scores = ['flux']
|
||||
model.tallies = [tally]
|
||||
|
||||
return model
|
||||
|
||||
|
||||
def test_lattice_large_pitch(model):
|
||||
harness = PyAPITestHarness('statepoint.10.h5', model)
|
||||
harness.main()
|
||||
|
|
@ -119,6 +119,61 @@ def test_weightwindows(shared_secondary, subdir):
|
|||
test.main()
|
||||
|
||||
|
||||
def test_zero_bound_windows_play_no_game(tmp_path):
|
||||
# A weight window lower bound of zero means no weight window information
|
||||
# exists there (MCNP wwinp files use zero to turn the game off in a cell),
|
||||
# so transport must proceed as if weight windows were disabled. Previously,
|
||||
# zero-bound windows demanded a split at every checkpoint (weight/0 ->
|
||||
# max_split), multiplying the particle population until terminated by the
|
||||
# split or weight cutoff limits.
|
||||
model = build_model(False)
|
||||
for ww in model.settings.weight_windows:
|
||||
ww.lower_ww_bounds = np.zeros_like(ww.lower_ww_bounds)
|
||||
ww.upper_ww_bounds = np.zeros_like(ww.upper_ww_bounds)
|
||||
sp_zero = model.run(cwd=tmp_path / 'zero_windows')
|
||||
|
||||
model.settings.weight_windows_on = False
|
||||
sp_off = model.run(cwd=tmp_path / 'windows_off')
|
||||
|
||||
with openmc.StatePoint(sp_zero) as sp:
|
||||
flux_zero = list(sp.tallies.values())[0].mean
|
||||
with openmc.StatePoint(sp_off) as sp:
|
||||
flux_off = list(sp.tallies.values())[0].mean
|
||||
|
||||
np.testing.assert_allclose(flux_zero, flux_off, rtol=1e-12)
|
||||
|
||||
|
||||
def test_zero_and_negative_bounds_equivalent(tmp_path):
|
||||
# Zero and negative lower bounds both mean that no weight window
|
||||
# information exists in a cell (generators mark such cells with -1, and
|
||||
# MCNP wwinp files use zero), so they must produce identical transport.
|
||||
# Unlike the all-zero case above, here particles are born under valid
|
||||
# windows and encounter the no-information region in flight; previously a
|
||||
# zero lower bound in that situation demanded a split at every checkpoint
|
||||
# in the cell (weight/0 -> max_split), multiplying the particle population,
|
||||
# while -1 played no game.
|
||||
def run_with(bound_value, subdir):
|
||||
model = build_model(False)
|
||||
for ww in model.settings.weight_windows:
|
||||
lb = np.array(ww.lower_ww_bounds, copy=True)
|
||||
ub = np.array(ww.upper_ww_bounds, copy=True)
|
||||
lb[3:, :, :, :] = bound_value
|
||||
ub[3:, :, :, :] = bound_value
|
||||
ww.lower_ww_bounds = lb
|
||||
ww.upper_ww_bounds = ub
|
||||
return model.run(cwd=tmp_path / subdir)
|
||||
|
||||
sp_zero = run_with(0.0, 'zero_region')
|
||||
sp_negative = run_with(-1.0, 'negative_region')
|
||||
|
||||
with openmc.StatePoint(sp_zero) as sp:
|
||||
flux_zero = list(sp.tallies.values())[0].mean
|
||||
with openmc.StatePoint(sp_negative) as sp:
|
||||
flux_negative = list(sp.tallies.values())[0].mean
|
||||
|
||||
np.testing.assert_allclose(flux_zero, flux_negative, rtol=1e-12)
|
||||
|
||||
|
||||
def test_wwinp_cylindrical():
|
||||
|
||||
ww = openmc.WeightWindowsList.from_wwinp('ww_n_cyl.txt')[0]
|
||||
|
|
|
|||
|
|
@ -155,16 +155,16 @@ def test_decay_photon_energy():
|
|||
with pytest.raises(DataError):
|
||||
openmc.data.decay_photon_energy('I135')
|
||||
|
||||
# Set chain file to simple chain
|
||||
openmc.config['chain_file'] = Path(__file__).parents[1] / "chain_simple.xml"
|
||||
# Temporarily Set chain file to simple chain
|
||||
with openmc.config.patch('chain_file', Path(__file__).parents[1] / 'chain_simple.xml'):
|
||||
|
||||
# Check strength of I135 source and presence of specific spectral line
|
||||
src = openmc.data.decay_photon_energy('I135')
|
||||
assert isinstance(src, openmc.stats.Discrete)
|
||||
assert src.integral() == pytest.approx(3.920996223799345e-05)
|
||||
assert 1260409. in src.x
|
||||
# Check strength of I135 source and presence of specific spectral line
|
||||
src = openmc.data.decay_photon_energy('I135')
|
||||
assert isinstance(src, openmc.stats.Discrete)
|
||||
assert src.integral() == pytest.approx(3.920996223799345e-05)
|
||||
assert 1260409. in src.x
|
||||
|
||||
# Check Xe135 source, which should be tabular
|
||||
src = openmc.data.decay_photon_energy('Xe135')
|
||||
assert isinstance(src, openmc.stats.Tabular)
|
||||
assert src.integral() == pytest.approx(2.076506258964966e-05)
|
||||
# Check Xe135 source, which should be tabular
|
||||
src = openmc.data.decay_photon_energy('Xe135')
|
||||
assert isinstance(src, openmc.stats.Tabular)
|
||||
assert src.integral() == pytest.approx(2.076506258964966e-05)
|
||||
|
|
|
|||
|
|
@ -34,6 +34,20 @@ def test_dose_coefficients():
|
|||
assert energy[-1] == approx(20.0e6)
|
||||
assert dose[-1] == approx(338.0)
|
||||
|
||||
energy, dose = dose_coefficients(
|
||||
'neutron', data_source='icrp74', dose_quantity='ambient')
|
||||
assert energy[0] == approx(1e-3)
|
||||
assert dose[0] == approx(6.60)
|
||||
assert energy[-1] == approx(20.0e6)
|
||||
assert dose[-1] == approx(600)
|
||||
|
||||
energy, dose = dose_coefficients(
|
||||
'photon', data_source='icrp74', dose_quantity='ambient')
|
||||
assert energy[0] == approx(0.01e6)
|
||||
assert dose[0] == approx(0.061)
|
||||
assert energy[-1] == approx(10e6)
|
||||
assert dose[-1] == approx(25.6)
|
||||
|
||||
# Invalid particle/geometry should raise an exception
|
||||
with raises(ValueError):
|
||||
dose_coefficients('slime', 'LAT')
|
||||
|
|
@ -41,6 +55,14 @@ def test_dose_coefficients():
|
|||
dose_coefficients('neutron', 'ZZ')
|
||||
with raises(ValueError):
|
||||
dose_coefficients('neutron', data_source='icrp7000')
|
||||
with raises(ValueError):
|
||||
dose_coefficients('neutron', dose_quantity='banana')
|
||||
with raises(ValueError):
|
||||
dose_coefficients(
|
||||
'neutron', data_source='icrp116', dose_quantity='ambient')
|
||||
with raises(ValueError):
|
||||
dose_coefficients(
|
||||
'neutron', 'ISO', data_source='icrp74', dose_quantity='ambient')
|
||||
with raises(ValueError) as excinfo:
|
||||
dose_coefficients("photons", data_source="icrp116")
|
||||
expected_particles = [
|
||||
|
|
@ -57,7 +79,8 @@ def test_dose_coefficients():
|
|||
"proton",
|
||||
]
|
||||
expected_msg = (
|
||||
"'photons' has no dose data in data source icrp116. "
|
||||
f"Available particles for icrp116 are: {expected_particles}"
|
||||
"'photons' has no effective dose data in data source icrp116. "
|
||||
"Available particles for icrp116 with dose quantity effective are: "
|
||||
f"{expected_particles}"
|
||||
)
|
||||
assert str(excinfo.value) == expected_msg
|
||||
|
|
|
|||
|
|
@ -7,6 +7,7 @@ from pathlib import Path
|
|||
import numpy as np
|
||||
import pytest
|
||||
import openmc.data
|
||||
from openmc.deplete import Chain, Nuclide
|
||||
|
||||
|
||||
def test_data_library(tmpdir):
|
||||
|
|
@ -134,7 +135,8 @@ def test_zam():
|
|||
with pytest.raises(ValueError):
|
||||
openmc.data.zam('Am242-m1')
|
||||
|
||||
def test_half_life():
|
||||
|
||||
def test_half_life(tmp_path):
|
||||
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)
|
||||
|
|
@ -143,3 +145,32 @@ def test_half_life():
|
|||
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)
|
||||
|
||||
# Create minimal chain with H3 and Am242 to test half-life and decay
|
||||
# constant retrieval from chain file
|
||||
chain = Chain()
|
||||
h3 = Nuclide("H3")
|
||||
h3.half_life = 1.0
|
||||
chain.add_nuclide(h3)
|
||||
am242 = Nuclide("Am242")
|
||||
chain.add_nuclide(am242)
|
||||
|
||||
assert openmc.data.half_life('H3', chain_file=chain) == 1.0
|
||||
assert openmc.data.decay_constant('H3', chain_file=chain) == pytest.approx(log(2.0))
|
||||
|
||||
# Nuclides that are present but stable in the chain should not fall back to
|
||||
# ENDF/B-VIII.0 data.
|
||||
assert openmc.data.half_life('Am242', chain_file=chain) is None
|
||||
assert openmc.data.decay_constant('Am242', chain_file=chain) == 0.0
|
||||
|
||||
# Nuclides missing from the chain fall back to ENDF/B-VIII.0 data.
|
||||
assert openmc.data.half_life('U235', chain_file=chain) == pytest.approx(2.22102e16)
|
||||
|
||||
chain_path = tmp_path / "chain.xml"
|
||||
chain.export_to_xml(chain_path)
|
||||
assert openmc.data.half_life('H3', chain_file=chain_path) == 1.0
|
||||
|
||||
endf_h3 = openmc.data.half_life('H3')
|
||||
with openmc.config.patch('chain_file', chain_path):
|
||||
assert openmc.data.half_life('H3', chain_file=None) == 1.0
|
||||
assert openmc.data.half_life('H3', chain_file=False) == endf_h3
|
||||
|
|
|
|||
|
|
@ -6,6 +6,7 @@ from pathlib import Path
|
|||
import numpy as np
|
||||
import pytest
|
||||
|
||||
import openmc
|
||||
import openmc.deplete
|
||||
|
||||
|
||||
|
|
@ -38,6 +39,36 @@ def test_get_activity(res):
|
|||
np.testing.assert_allclose(a_xe135, a_xe135_ref)
|
||||
|
||||
|
||||
def test_get_activity_chain_file(res, tmp_path):
|
||||
"""Tests evaluating activity with chain half-life data"""
|
||||
_, a_endf = res.get_activity("1", by_nuclide=True, chain_file=False)
|
||||
xe135_endf = np.array([a["Xe135"] for a in a_endf])
|
||||
|
||||
chain = openmc.deplete.Chain()
|
||||
xe135 = openmc.deplete.Nuclide("Xe135")
|
||||
xe135.half_life = openmc.data.half_life("Xe135") / 2.0
|
||||
chain.add_nuclide(xe135)
|
||||
|
||||
t_chain, a_chain = res.get_activity("1", by_nuclide=True, chain_file=chain)
|
||||
xe135_chain = np.array([a["Xe135"] for a in a_chain])
|
||||
|
||||
t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0])
|
||||
np.testing.assert_allclose(t_chain, t_ref)
|
||||
np.testing.assert_allclose(xe135_chain, 2.0 * xe135_endf)
|
||||
|
||||
chain_path = tmp_path / "chain.xml"
|
||||
chain.export_to_xml(chain_path)
|
||||
with openmc.config.patch('chain_file', chain_path):
|
||||
_, a_config = res.get_activity("1", by_nuclide=True)
|
||||
xe135_config = np.array([a["Xe135"] for a in a_config])
|
||||
np.testing.assert_allclose(xe135_config, xe135_chain)
|
||||
|
||||
stable_chain = openmc.deplete.Chain()
|
||||
stable_chain.add_nuclide(openmc.deplete.Nuclide("Xe135"))
|
||||
_, a_stable = res.get_activity("1", by_nuclide=True, chain_file=stable_chain)
|
||||
assert all(a["Xe135"] == 0.0 for a in a_stable)
|
||||
|
||||
|
||||
def test_get_atoms(res):
|
||||
"""Tests evaluating single nuclide concentration."""
|
||||
t, n = res.get_atoms("1", "Xe135")
|
||||
|
|
|
|||
|
|
@ -926,6 +926,22 @@ def test_property_map(lib_init):
|
|||
assert np.allclose(expected_properties, properties, atol=1e-04)
|
||||
|
||||
|
||||
def test_slice_data(lib_init):
|
||||
expected_properties = np.array(
|
||||
[[(293.6, 0.740582), (293.6, 6.55), (293.6, 0.740582)],
|
||||
[ (293.6, 6.55), (293.6, 10.29769), (293.6, 6.55)],
|
||||
[(293.6, 0.740582), (293.6, 6.55), (293.6, 0.740582)]], dtype='float')
|
||||
origin = (0.0, 0.0, 0.0)
|
||||
_, properties = openmc.lib.slice_data(
|
||||
origin,
|
||||
width=(1.26, 1.26),
|
||||
basis='xy',
|
||||
pixels=(3, 3),
|
||||
include_properties=True
|
||||
)
|
||||
assert np.allclose(expected_properties, properties, atol=1e-04)
|
||||
|
||||
|
||||
def test_solid_raytrace_plot(lib_init, pincell_model):
|
||||
# Ensure plot mapping can be accessed and grows after allocation
|
||||
n0 = len(openmc.lib.plots)
|
||||
|
|
|
|||
|
|
@ -7,7 +7,7 @@ import numpy as np
|
|||
|
||||
import openmc
|
||||
from openmc.data import decay_photon_energy
|
||||
from openmc.deplete import Chain
|
||||
from openmc.deplete import Chain, Nuclide
|
||||
import openmc.examples
|
||||
import openmc.model
|
||||
import openmc.stats
|
||||
|
|
@ -614,6 +614,35 @@ def test_get_activity():
|
|||
assert m4.get_activity(units='Ci/m3') == pytest.approx(ci/m3)
|
||||
|
||||
|
||||
def test_get_activity_chain_file(tmp_path):
|
||||
m = openmc.Material()
|
||||
m.add_nuclide("H3", 1.0)
|
||||
m.set_density('g/cm3', 1.0)
|
||||
|
||||
chain = Chain()
|
||||
h3 = Nuclide("H3")
|
||||
h3.half_life = 1.0
|
||||
chain.add_nuclide(h3)
|
||||
|
||||
atoms_per_bcm = m.get_nuclide_atom_densities()["H3"]
|
||||
expected = np.log(2.0) * 1e24 * atoms_per_bcm
|
||||
|
||||
assert m.get_activity(chain_file=chain) == pytest.approx(expected)
|
||||
|
||||
chain_path = tmp_path / "chain.xml"
|
||||
chain.export_to_xml(chain_path)
|
||||
assert m.get_activity(chain_file=chain_path) == pytest.approx(expected)
|
||||
|
||||
endf_activity = m.get_activity(chain_file=False)
|
||||
with openmc.config.patch('chain_file', chain_path):
|
||||
assert m.get_activity() == pytest.approx(expected)
|
||||
assert m.get_activity(chain_file=False) == pytest.approx(endf_activity)
|
||||
|
||||
stable_chain = Chain()
|
||||
stable_chain.add_nuclide(Nuclide("H3"))
|
||||
assert m.get_activity(chain_file=stable_chain) == 0.0
|
||||
|
||||
|
||||
def test_get_decay_heat():
|
||||
# Set chain file for testing
|
||||
openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml'
|
||||
|
|
|
|||
|
|
@ -660,6 +660,29 @@ def test_model_plot():
|
|||
plt.close('all')
|
||||
|
||||
|
||||
def test_model_plot_invalid_inputs():
|
||||
surface = openmc.Sphere(r=10.0, boundary_type="vacuum")
|
||||
cell = openmc.Cell(region=-surface)
|
||||
model = openmc.Model(openmc.Geometry([cell]))
|
||||
|
||||
with pytest.raises(ValueError):
|
||||
model.plot(n_samples=-1)
|
||||
with pytest.raises(TypeError):
|
||||
model.plot(n_samples=1.5)
|
||||
with pytest.raises(ValueError):
|
||||
model.plot(plane_tolerance=0.0)
|
||||
with pytest.raises(TypeError):
|
||||
model.plot(plane_tolerance='1')
|
||||
with pytest.raises(ValueError):
|
||||
model.plot(pixels=-1)
|
||||
with pytest.raises(ValueError):
|
||||
model.plot(pixels=(0, 100))
|
||||
with pytest.raises(ValueError):
|
||||
model.plot(pixels=(100,))
|
||||
with pytest.raises(ValueError):
|
||||
model.slice_data(u_span=(2, 0, 0), v_span=(0, 2, 0), pixels=-1)
|
||||
|
||||
|
||||
def test_model_id_map_initialization(run_in_tmpdir):
|
||||
model = openmc.examples.pwr_assembly()
|
||||
model.init_lib(output=False)
|
||||
|
|
|
|||
59
tests/unit_tests/test_pulse_height.py
Normal file
59
tests/unit_tests/test_pulse_height.py
Normal file
|
|
@ -0,0 +1,59 @@
|
|||
import numpy as np
|
||||
import pytest
|
||||
import openmc
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def model():
|
||||
openmc.reset_auto_ids()
|
||||
model = openmc.Model()
|
||||
|
||||
# Define materials
|
||||
NaI = openmc.Material()
|
||||
NaI.set_density('g/cm3', 3.7)
|
||||
NaI.add_element('Na', 1.0)
|
||||
NaI.add_element('I', 1.0)
|
||||
|
||||
# Define geometry: two spheres in each other
|
||||
s1 = openmc.Sphere(r=1)
|
||||
s2 = openmc.Sphere(r=2, boundary_type='vacuum')
|
||||
inner_sphere = openmc.Cell(name='inner sphere', fill=NaI, region=-s1)
|
||||
outer_sphere = openmc.Cell(name='outer sphere', fill=NaI, region=+s1 & -s2)
|
||||
model.geometry = openmc.Geometry([inner_sphere, outer_sphere])
|
||||
|
||||
# Define settings
|
||||
model.settings.run_mode = 'fixed source'
|
||||
model.settings.batches = 1
|
||||
model.settings.particles = 10000
|
||||
model.settings.photon_transport = True
|
||||
model.settings.source = openmc.IndependentSource(
|
||||
energy=openmc.stats.delta_function(1e6),
|
||||
particle='photon'
|
||||
)
|
||||
|
||||
# Define tallies
|
||||
energy_filter = openmc.EnergyFilter([1e3, 1e7])
|
||||
|
||||
tally1 = openmc.Tally()
|
||||
tally1.scores = ['pulse-height']
|
||||
cell_filter1 = openmc.CellFilter([inner_sphere, outer_sphere])
|
||||
tally1.filters = [cell_filter1, energy_filter]
|
||||
|
||||
tally2 = openmc.Tally()
|
||||
tally2.scores = ['pulse-height']
|
||||
cell_filter2 = openmc.CellFilter([outer_sphere, inner_sphere])
|
||||
tally2.filters = [cell_filter2, energy_filter]
|
||||
|
||||
model.tallies = [tally1, tally2]
|
||||
return model
|
||||
|
||||
|
||||
def test_pulse_height(model, run_in_tmpdir):
|
||||
sp_path = model.run()
|
||||
sp = openmc.StatePoint(sp_path)
|
||||
t1 = sp.tallies[1].mean.squeeze()
|
||||
t2 = sp.tallies[2].mean.squeeze()
|
||||
|
||||
np.testing.assert_array_equal(t1, t2[::-1])
|
||||
|
||||
|
||||
|
|
@ -33,8 +33,8 @@ def simple_model_and_mesh():
|
|||
|
||||
# Simple settings with a point source
|
||||
settings = openmc.Settings()
|
||||
settings.batches = 10
|
||||
settings.particles = 1000
|
||||
settings.batches = 2
|
||||
settings.particles = 250
|
||||
settings.run_mode = 'fixed source'
|
||||
settings.source = openmc.IndependentSource()
|
||||
model = openmc.Model(geometry, settings=settings)
|
||||
|
|
@ -46,6 +46,16 @@ def simple_model_and_mesh():
|
|||
return model, (c1, c2), mesh
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def source_stage_manager(simple_model_and_mesh):
|
||||
model, (c1, c2), _ = simple_model_and_mesh
|
||||
r2s = R2SManager(model, [c1, c2])
|
||||
r2s.results['depletion_results'] = [None, None]
|
||||
r2s.results['activation_materials'] = [c1.fill, c2.fill]
|
||||
bounding_boxes = {c1.id: c1.bounding_box, c2.id: c2.bounding_box}
|
||||
return r2s, bounding_boxes
|
||||
|
||||
|
||||
def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
|
||||
model, (c1, c2), mesh = simple_model_and_mesh
|
||||
|
||||
|
|
@ -59,9 +69,9 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
|
|||
outdir = r2s.run(
|
||||
timesteps=[(1.0, 'd')],
|
||||
source_rates=[1.0],
|
||||
photon_time_indices=[1],
|
||||
output_dir=tmp_path,
|
||||
chain_file=chain,
|
||||
micro_kwargs={'nuclides': ['Ni58'], 'reactions': ['(n,p)']},
|
||||
)
|
||||
|
||||
# Check directories and files exist
|
||||
|
|
@ -73,7 +83,8 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
|
|||
assert (act / 'depletion_results.h5').exists()
|
||||
pt = Path(outdir) / 'photon_transport'
|
||||
assert (pt / 'tally_ids.json').exists()
|
||||
assert (pt / 'time_1' / 'statepoint.10.h5').exists()
|
||||
assert not (pt / 'time_0').exists()
|
||||
assert (pt / 'time_1' / 'statepoint.2.h5').exists()
|
||||
|
||||
# Basic results structure checks
|
||||
assert len(r2s.results['fluxes']) == 2
|
||||
|
|
@ -82,6 +93,8 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
|
|||
assert len(r2s.results['mesh_material_volumes'][0]) == 2
|
||||
assert len(r2s.results['activation_materials']) == 2
|
||||
assert len(r2s.results['depletion_results']) == 2
|
||||
assert list(r2s.results['photon_sources']) == [1]
|
||||
assert r2s.results['photon_sources'][1]
|
||||
|
||||
# Check activation materials
|
||||
amats = r2s.results['activation_materials']
|
||||
|
|
@ -125,6 +138,7 @@ def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path):
|
|||
photon_time_indices=[1],
|
||||
output_dir=tmp_path,
|
||||
chain_file=chain,
|
||||
micro_kwargs={'nuclides': ['Ni58'], 'reactions': ['(n,p)']},
|
||||
)
|
||||
|
||||
# Check that per-mesh MMV files were written
|
||||
|
|
@ -137,7 +151,7 @@ def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path):
|
|||
assert (act / 'depletion_results.h5').exists()
|
||||
pt = Path(outdir) / 'photon_transport'
|
||||
assert (pt / 'tally_ids.json').exists()
|
||||
assert (pt / 'time_1' / 'statepoint.10.h5').exists()
|
||||
assert (pt / 'time_1' / 'statepoint.2.h5').exists()
|
||||
|
||||
# Two meshes, each with 1 element containing both materials →
|
||||
# 2 element-material combinations per mesh, 4 total
|
||||
|
|
@ -167,6 +181,9 @@ def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path):
|
|||
|
||||
def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
|
||||
model, (c1, c2), _ = simple_model_and_mesh
|
||||
tally = openmc.Tally()
|
||||
tally.scores = ['flux']
|
||||
model.tallies = [tally]
|
||||
|
||||
# Use cell-based domains
|
||||
r2s = R2SManager(model, [c1, c2])
|
||||
|
|
@ -180,9 +197,11 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
|
|||
timesteps=[(1.0, 'd')],
|
||||
source_rates=[1.0],
|
||||
photon_time_indices=[1],
|
||||
by_parent_nuclide=True,
|
||||
output_dir=tmp_path,
|
||||
bounding_boxes=bounding_boxes,
|
||||
chain_file=chain
|
||||
chain_file=chain,
|
||||
micro_kwargs={'nuclides': ['Ni58'], 'reactions': ['(n,p)']},
|
||||
)
|
||||
|
||||
# Check directories and files exist
|
||||
|
|
@ -193,13 +212,15 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
|
|||
assert (act / 'depletion_results.h5').exists()
|
||||
pt = Path(outdir) / 'photon_transport'
|
||||
assert (pt / 'tally_ids.json').exists()
|
||||
assert (pt / 'time_1' / 'statepoint.10.h5').exists()
|
||||
assert (pt / 'time_1' / 'statepoint.2.h5').exists()
|
||||
|
||||
# Basic results structure checks
|
||||
assert len(r2s.results['fluxes']) == 2
|
||||
assert len(r2s.results['micros']) == 2
|
||||
assert len(r2s.results['activation_materials']) == 2
|
||||
assert len(r2s.results['depletion_results']) == 2
|
||||
assert r2s.photon_model.tallies[0].contains_filter(
|
||||
openmc.ParentNuclideFilter)
|
||||
|
||||
# Check activation materials
|
||||
amats = r2s.results['activation_materials']
|
||||
|
|
@ -217,3 +238,88 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
|
|||
assert len(r2s_loaded.results['micros']) == 2
|
||||
assert len(r2s_loaded.results['activation_materials']) == 2
|
||||
assert len(r2s_loaded.results['depletion_results']) == 2
|
||||
|
||||
|
||||
def test_step4_requires_photon_sources(simple_model_and_mesh, tmp_path):
|
||||
model, (c1, c2), _ = simple_model_and_mesh
|
||||
r2s = R2SManager(model, [c1, c2])
|
||||
output_dir = tmp_path / 'photon'
|
||||
|
||||
with pytest.raises(RuntimeError, match='step3_photon_source'):
|
||||
r2s.step4_photon_transport(output_dir)
|
||||
|
||||
r2s.results['photon_sources'] = {}
|
||||
with pytest.raises(RuntimeError, match='No decay photon sources'):
|
||||
r2s.step4_photon_transport(output_dir)
|
||||
|
||||
assert not output_dir.exists()
|
||||
|
||||
|
||||
def test_default_photon_times_skip_empty_sources(
|
||||
source_stage_manager, tmp_path, monkeypatch
|
||||
):
|
||||
r2s, bounding_boxes = source_stage_manager
|
||||
source = object()
|
||||
sources_by_time = {0: [], 1: [source]}
|
||||
monkeypatch.setattr(
|
||||
r2s, '_create_photon_sources',
|
||||
lambda time_index, work_items: sources_by_time[time_index])
|
||||
|
||||
r2s.step3_photon_source(
|
||||
bounding_boxes=bounding_boxes, output_dir=tmp_path)
|
||||
|
||||
assert r2s.results['photon_sources'] == {1: [source]}
|
||||
|
||||
|
||||
def test_explicit_empty_photon_source_fails(
|
||||
source_stage_manager, tmp_path, monkeypatch
|
||||
):
|
||||
r2s, bounding_boxes = source_stage_manager
|
||||
source = object()
|
||||
sources_by_time = {0: [], 1: [source]}
|
||||
monkeypatch.setattr(
|
||||
r2s, '_create_photon_sources',
|
||||
lambda time_index, work_items: sources_by_time[time_index])
|
||||
r2s.results['photon_sources'] = {99: [source]}
|
||||
|
||||
with pytest.raises(RuntimeError, match='requested time indices: 0'):
|
||||
r2s.step3_photon_source(
|
||||
[0, 1], bounding_boxes, output_dir=tmp_path)
|
||||
|
||||
assert 'photon_sources' not in r2s.results
|
||||
|
||||
|
||||
def test_default_photon_times_require_a_source(
|
||||
source_stage_manager, tmp_path, monkeypatch
|
||||
):
|
||||
r2s, bounding_boxes = source_stage_manager
|
||||
monkeypatch.setattr(
|
||||
r2s, '_create_photon_sources',
|
||||
lambda time_index, work_items: [])
|
||||
|
||||
with pytest.raises(RuntimeError, match='at any depletion time'):
|
||||
r2s.step3_photon_source(
|
||||
bounding_boxes=bounding_boxes, output_dir=tmp_path)
|
||||
|
||||
assert 'photon_sources' not in r2s.results
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
('time_indices', 'exception'),
|
||||
[
|
||||
([], ValueError),
|
||||
([2], IndexError),
|
||||
([-3], IndexError),
|
||||
([1.0], TypeError),
|
||||
],
|
||||
)
|
||||
def test_photon_time_index_validation(
|
||||
source_stage_manager, tmp_path, time_indices, exception
|
||||
):
|
||||
r2s, bounding_boxes = source_stage_manager
|
||||
|
||||
with pytest.raises(exception):
|
||||
r2s.step3_photon_source(
|
||||
time_indices, bounding_boxes, output_dir=tmp_path)
|
||||
|
||||
assert 'photon_sources' not in r2s.results
|
||||
|
|
|
|||
89
tests/unit_tests/test_slice_data_overlap_info.py
Normal file
89
tests/unit_tests/test_slice_data_overlap_info.py
Normal file
|
|
@ -0,0 +1,89 @@
|
|||
import pytest
|
||||
import numpy as np
|
||||
import openmc
|
||||
import openmc.lib
|
||||
|
||||
# Sentinel value matching _OVERLAP in plotmodel.py and OVERLAP in plot.cpp
|
||||
_OVERLAP = -3
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def overlap_model():
|
||||
openmc.reset_auto_ids()
|
||||
|
||||
# Three cylinders: cyl1 and cyl2 overlap near x=0, cyl2 and cyl3 overlap
|
||||
# near x=4. This gives us two spatially distinct overlap regions in one model.
|
||||
mat1 = openmc.Material(components={'H1': 1.0})
|
||||
mat2 = openmc.Material(components={'H1': 1.0})
|
||||
mat3 = openmc.Material(components={'H1': 1.0})
|
||||
|
||||
# cyl1 and cyl2 overlap on the left, cyl2 and cyl3 overlap on the right
|
||||
cyl1 = openmc.ZCylinder(x0=-2.0, r=2.5)
|
||||
cyl2 = openmc.ZCylinder(x0=0.0, r=2.5)
|
||||
cyl3 = openmc.ZCylinder(x0=2.0, r=2.5)
|
||||
boundary = openmc.Sphere(r=20.0, boundary_type='vacuum')
|
||||
cell1 = openmc.Cell(region=-cyl1, fill=mat1)
|
||||
cell2 = openmc.Cell(region=-cyl2, fill=mat2)
|
||||
cell3 = openmc.Cell(region=-cyl3, fill=mat3)
|
||||
cell_outside = openmc.Cell(region=+cyl1 & +cyl2 & +cyl3 & -boundary)
|
||||
geometry = openmc.Geometry([cell1, cell2, cell3, cell_outside])
|
||||
|
||||
settings = openmc.Settings()
|
||||
settings.run_mode = 'fixed source'
|
||||
settings.particles = 100
|
||||
settings.batches = 1
|
||||
model = openmc.Model(geometry=geometry, settings=settings)
|
||||
|
||||
with openmc.lib.TemporarySession(model, args=['-s', '1']):
|
||||
yield
|
||||
|
||||
|
||||
def run_slice(origin=(0.0, 0.0, 0.0), width=(10.0, 6.0), show_overlaps=True):
|
||||
# Helper that runs a slice over a region covering both overlap zones
|
||||
geom_data, _ = openmc.lib.slice_data(
|
||||
origin=origin,
|
||||
width=width,
|
||||
basis='xy',
|
||||
pixels=(100, 60),
|
||||
show_overlaps=show_overlaps,
|
||||
include_properties=False,
|
||||
)
|
||||
return geom_data
|
||||
|
||||
|
||||
def test_overlaps_enabled(overlap_model):
|
||||
# Run a single slice with overlap detection enabled and check all
|
||||
# expected properties in one pass.
|
||||
geom_data = run_slice()
|
||||
overlap_info = openmc.lib.slice_data_overlap_info()
|
||||
n = overlap_info.shape[0]
|
||||
cell_ids = geom_data[:, :, 0]
|
||||
|
||||
# cell_ids should contain values more negative than _OVERLAP; RasterData
|
||||
# encodes each unique overlap as OVERLAP - overlap_idx - 1 into slot 2.
|
||||
assert np.any(cell_ids < _OVERLAP)
|
||||
|
||||
# overlap_keys should have 2 entries for the two distinct overlapping
|
||||
# cylinder pairs in this model.
|
||||
assert n == 2, f"Expected exactly 2 overlap entries, got {n}"
|
||||
|
||||
# Each entry is a (universe_id, cell1_id, cell2_id) triple; verify values.
|
||||
for i in range(n):
|
||||
universe_id = int(overlap_info[i, 0])
|
||||
cell1_id = int(overlap_info[i, 1])
|
||||
cell2_id = int(overlap_info[i, 2])
|
||||
assert universe_id == 1
|
||||
assert cell1_id in {1, 2, 3}
|
||||
assert cell2_id in {1, 2, 3}
|
||||
assert cell1_id != cell2_id
|
||||
|
||||
|
||||
def test_overlaps_disabled(overlap_model):
|
||||
# With show_overlaps=False, set_overlap is never called and overlap_keys
|
||||
# is never written to, so the image and map should both be clean.
|
||||
geom_data = run_slice(show_overlaps=False)
|
||||
overlap_info = openmc.lib.slice_data_overlap_info()
|
||||
n = overlap_info.shape[0]
|
||||
|
||||
assert not np.any(geom_data[:, :, 2] < _OVERLAP)
|
||||
assert n == 0
|
||||
263
tests/unit_tests/test_source_tokamak.py
Normal file
263
tests/unit_tests/test_source_tokamak.py
Normal file
|
|
@ -0,0 +1,263 @@
|
|||
import numpy as np
|
||||
import pytest
|
||||
|
||||
import openmc
|
||||
import openmc.stats
|
||||
|
||||
from tests.unit_tests import assert_sample_mean
|
||||
|
||||
|
||||
def make_source(**kwargs):
|
||||
"""Build a valid TokamakSource, overriding defaults via kwargs."""
|
||||
r_over_a = np.linspace(0.0, 1.0, 10)
|
||||
params = dict(
|
||||
major_radius=620.0,
|
||||
minor_radius=200.0,
|
||||
elongation=1.8,
|
||||
triangularity=0.45,
|
||||
shafranov_shift=10.0,
|
||||
r_over_a=r_over_a,
|
||||
emission_density=(1.0 - r_over_a**2),
|
||||
energy=openmc.stats.muir(e0=14.08e6, m_rat=5.0, kt=2.0e4),
|
||||
)
|
||||
params.update(kwargs)
|
||||
return openmc.TokamakSource(**params)
|
||||
|
||||
|
||||
def test_tokamak_source_roundtrip():
|
||||
src = make_source(
|
||||
phi_start=0.1, phi_extent=np.pi, n_alpha=51, vertical_shift=5.0,
|
||||
strength=2.0, time=openmc.stats.Uniform(0.0, 1e-6))
|
||||
|
||||
elem = src.to_xml_element()
|
||||
assert elem.get('type') == 'tokamak'
|
||||
|
||||
new = openmc.SourceBase.from_xml_element(elem)
|
||||
assert isinstance(new, openmc.TokamakSource)
|
||||
assert new.major_radius == src.major_radius
|
||||
assert new.minor_radius == src.minor_radius
|
||||
assert new.elongation == src.elongation
|
||||
assert new.triangularity == src.triangularity
|
||||
assert new.shafranov_shift == src.shafranov_shift
|
||||
assert new.phi_start == src.phi_start
|
||||
assert new.phi_extent == src.phi_extent
|
||||
assert new.n_alpha == src.n_alpha
|
||||
assert new.vertical_shift == src.vertical_shift
|
||||
assert new.strength == src.strength
|
||||
np.testing.assert_allclose(new.r_over_a, src.r_over_a)
|
||||
np.testing.assert_allclose(new.emission_density, src.emission_density)
|
||||
assert len(new.energy) == 1
|
||||
assert isinstance(new.time, openmc.stats.Uniform)
|
||||
assert new.time.a == src.time.a
|
||||
assert new.time.b == src.time.b
|
||||
|
||||
|
||||
def test_tokamak_source_default_time():
|
||||
src = make_source()
|
||||
assert src.time is None
|
||||
|
||||
new = openmc.SourceBase.from_xml_element(src.to_xml_element())
|
||||
assert new.time is None
|
||||
|
||||
with pytest.raises(TypeError):
|
||||
make_source(time=1.0)
|
||||
|
||||
|
||||
def test_tokamak_source_multiple_energies():
|
||||
r_over_a = np.linspace(0.0, 1.0, 5)
|
||||
energies = [openmc.stats.muir(e0=14.08e6, m_rat=5.0, kt=kt)
|
||||
for kt in (1.0e4, 1.5e4, 2.0e4, 2.5e4, 3.0e4)]
|
||||
src = make_source(r_over_a=r_over_a,
|
||||
emission_density=np.ones_like(r_over_a),
|
||||
energy=energies)
|
||||
assert len(src.energy) == len(r_over_a)
|
||||
|
||||
new = openmc.SourceBase.from_xml_element(src.to_xml_element())
|
||||
assert len(new.energy) == len(r_over_a)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("kwargs, match", [
|
||||
(dict(minor_radius=700.0), "smaller than major_radius"),
|
||||
(dict(shafranov_shift=150.0), "half the minor_radius"),
|
||||
(dict(emission_density=np.ones(5)), "same length as r_over_a"),
|
||||
(dict(energy=[openmc.stats.muir(14.08e6, 5.0, 2.0e4)] * 2),
|
||||
"Number of energy distributions"),
|
||||
(dict(r_over_a=np.linspace(0.1, 1.0, 10)), "must start at 0"),
|
||||
(dict(r_over_a=np.linspace(0.0, 0.9, 10)), "must end at 1"),
|
||||
(dict(emission_density=-np.linspace(0.0, 1.0, 10)), "cannot be negative"),
|
||||
(dict(emission_density=np.zeros(10)), "must contain a positive value"),
|
||||
])
|
||||
def test_tokamak_source_invalid(kwargs, match):
|
||||
with pytest.raises(ValueError, match=match):
|
||||
make_source(**kwargs)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("value", [-1.5, 1.5])
|
||||
def test_tokamak_source_invalid_triangularity(value):
|
||||
with pytest.raises(ValueError):
|
||||
make_source(triangularity=value)
|
||||
|
||||
|
||||
def test_tokamak_source_invalid_n_alpha():
|
||||
with pytest.raises(ValueError):
|
||||
make_source(n_alpha=2)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(("attribute", "value", "match"), [
|
||||
("minor_radius", 700.0, "smaller than major_radius"),
|
||||
("shafranov_shift", 150.0, "half the minor_radius"),
|
||||
("emission_density", np.ones(5), "same length as r_over_a"),
|
||||
("energy", [openmc.stats.delta_function(1.0)] * 2,
|
||||
"Number of energy distributions"),
|
||||
])
|
||||
def test_tokamak_source_mutation_validation(attribute, value, match):
|
||||
src = make_source()
|
||||
setattr(src, attribute, value)
|
||||
with pytest.raises(ValueError, match=match):
|
||||
src.to_xml_element()
|
||||
|
||||
|
||||
@pytest.mark.parametrize(("emission_density", "expected_mean"), [
|
||||
([1.0, 1.0], 2.0 / 3.0),
|
||||
([1.0, 0.0], 0.5),
|
||||
])
|
||||
def test_tokamak_source_radial_sampling(
|
||||
run_in_tmpdir, emission_density, expected_mean
|
||||
):
|
||||
"""Check radial sampling for profiles on the coarsest valid grid."""
|
||||
major_radius = 620.0
|
||||
minor_radius = 200.0
|
||||
src = make_source(
|
||||
major_radius=major_radius,
|
||||
minor_radius=minor_radius,
|
||||
elongation=1.0,
|
||||
triangularity=0.0,
|
||||
shafranov_shift=0.0,
|
||||
r_over_a=[0.0, 1.0],
|
||||
emission_density=emission_density,
|
||||
energy=openmc.stats.delta_function(14.07e6),
|
||||
)
|
||||
|
||||
sphere = openmc.Sphere(r=2000.0, boundary_type='vacuum')
|
||||
model = openmc.Model(
|
||||
geometry=openmc.Geometry([openmc.Cell(region=-sphere)]),
|
||||
settings=openmc.Settings(
|
||||
particles=100, batches=1, run_mode='fixed source', source=src),
|
||||
)
|
||||
|
||||
sites = model.sample_external_source(20_000)
|
||||
xyz = np.array([site.r for site in sites])
|
||||
major_r = np.hypot(xyz[:, 0], xyz[:, 1])
|
||||
r_over_a = np.hypot(major_r - major_radius, xyz[:, 2]) / minor_radius
|
||||
assert_sample_mean(r_over_a, expected_mean)
|
||||
|
||||
|
||||
def test_tokamak_source_poloidal_sampling(run_in_tmpdir):
|
||||
"""Check linear-linear poloidal sampling on a coarse internal grid."""
|
||||
major_radius = 620.0
|
||||
minor_radius = 200.0
|
||||
with pytest.warns(UserWarning, match="below 51"):
|
||||
src = make_source(
|
||||
major_radius=major_radius,
|
||||
minor_radius=minor_radius,
|
||||
elongation=1.0,
|
||||
triangularity=0.0,
|
||||
shafranov_shift=0.0,
|
||||
emission_density=np.ones(10),
|
||||
n_alpha=3,
|
||||
energy=openmc.stats.delta_function(14.07e6),
|
||||
)
|
||||
|
||||
sphere = openmc.Sphere(r=2000.0, boundary_type='vacuum')
|
||||
model = openmc.Model(
|
||||
geometry=openmc.Geometry([openmc.Cell(region=-sphere)]),
|
||||
settings=openmc.Settings(
|
||||
particles=100, batches=1, run_mode='fixed source', source=src),
|
||||
)
|
||||
|
||||
sites = model.sample_external_source(100_000)
|
||||
xyz = np.array([site.r for site in sites])
|
||||
major_r = np.hypot(xyz[:, 0], xyz[:, 1])
|
||||
|
||||
# With three alpha points, linear interpolation of cos(alpha) on [0, pi]
|
||||
# gives 1 - 2*alpha/pi. Integrating the resulting density gives this mean.
|
||||
expected_R = (
|
||||
major_radius
|
||||
+ 2.0 * minor_radius**2 / (np.pi**2 * major_radius)
|
||||
)
|
||||
assert_sample_mean(major_r, expected_R)
|
||||
|
||||
|
||||
@pytest.mark.flaky(reruns=1)
|
||||
def test_tokamak_source_sampling(run_in_tmpdir):
|
||||
"""Exercise the compiled C++ sampling path and check invariants.
|
||||
|
||||
Sampled moments are compared against direct numerical quadrature of the
|
||||
exact source density S(r)*R*|J|, where the Jacobian J of the flux-surface
|
||||
map is computed from analytic partial derivatives. This check is
|
||||
independent of the Bernstein-mixture factorization used by the
|
||||
implementation.
|
||||
"""
|
||||
R0, a, kappa, delta = 620.0, 200.0, 1.8, -0.5
|
||||
shafranov, zshift = 40.0, 25.0
|
||||
phi_start, phi_extent = 0.5, np.pi / 2
|
||||
|
||||
# Fine grids to make discretization error negligible relative to
|
||||
# statistical uncertainty
|
||||
r_over_a = np.linspace(0.0, 1.0, 200)
|
||||
src = make_source(
|
||||
major_radius=R0, minor_radius=a, elongation=kappa,
|
||||
triangularity=delta, shafranov_shift=shafranov, vertical_shift=zshift,
|
||||
r_over_a=r_over_a, emission_density=1.0 - r_over_a**2,
|
||||
phi_start=phi_start, phi_extent=phi_extent, n_alpha=201,
|
||||
energy=openmc.stats.delta_function(14.07e6))
|
||||
|
||||
sphere = openmc.Sphere(r=2000.0, boundary_type='vacuum')
|
||||
cell = openmc.Cell(region=-sphere)
|
||||
settings = openmc.Settings(
|
||||
particles=100, batches=1, run_mode='fixed source', source=src)
|
||||
model = openmc.Model(geometry=openmc.Geometry([cell]), settings=settings)
|
||||
|
||||
n_samples = 20_000
|
||||
sites = model.sample_external_source(n_samples)
|
||||
|
||||
xyz = np.array([s.r for s in sites])
|
||||
R = np.hypot(xyz[:, 0], xyz[:, 1])
|
||||
z = xyz[:, 2]
|
||||
|
||||
# Energy, weight, and time invariants
|
||||
assert np.all([s.E == 14.07e6 for s in sites])
|
||||
assert np.all([s.wgt == 1.0 for s in sites])
|
||||
assert np.all([s.time == 0.0 for s in sites])
|
||||
|
||||
# Toroidal angle within the requested sector
|
||||
phi = np.arctan2(xyz[:, 1], xyz[:, 0])
|
||||
assert phi.min() >= phi_start
|
||||
assert phi.max() <= phi_start + phi_extent
|
||||
|
||||
# Positions bounded by the last closed flux surface
|
||||
assert R.min() >= R0 - a
|
||||
assert R.max() <= R0 + a + shafranov
|
||||
assert np.abs(z - zshift).max() <= kappa * a
|
||||
|
||||
# Up-down symmetry about the vertical shift
|
||||
assert_sample_mean(z, zshift)
|
||||
|
||||
# Reference moments by 2D quadrature of the exact density
|
||||
r = np.linspace(0.0, 1.0, 1001)[:, np.newaxis] # r/a
|
||||
alpha = np.linspace(0.0, 2 * np.pi, 2001)[np.newaxis, :]
|
||||
psi = alpha + delta * np.sin(alpha)
|
||||
R_map = R0 + a * r * np.cos(psi) + shafranov * (1.0 - r**2)
|
||||
Z_map = kappa * a * r * np.sin(alpha)
|
||||
dR_dr = a * np.cos(psi) - 2.0 * shafranov * r
|
||||
dR_da = -a * r * np.sin(psi) * (1.0 + delta * np.cos(alpha))
|
||||
dZ_dr = kappa * a * np.sin(alpha) * np.ones_like(psi)
|
||||
dZ_da = kappa * a * r * np.cos(alpha)
|
||||
jac = np.abs(dR_dr * dZ_da - dR_da * dZ_dr)
|
||||
dens = (1.0 - r**2) * R_map * jac
|
||||
norm = dens.sum()
|
||||
expected_R = (R_map * dens).sum() / norm
|
||||
expected_z2 = (Z_map**2 * dens).sum() / norm
|
||||
|
||||
assert_sample_mean(R, expected_R)
|
||||
assert_sample_mean((z - zshift)**2, expected_z2)
|
||||
|
|
@ -98,6 +98,50 @@ def test_surface_filter_flux_angled(two_cell_model, run_in_tmpdir):
|
|||
assert flux_mean == pytest.approx(1.0 / mu)
|
||||
|
||||
|
||||
def test_surface_tally_during_lattice_crossing(run_in_tmpdir):
|
||||
openmc.reset_auto_ids()
|
||||
model = openmc.Model()
|
||||
|
||||
xmin = openmc.XPlane(-1.0, boundary_type="vacuum")
|
||||
xmax = openmc.XPlane(1.0, boundary_type="vacuum")
|
||||
ymin = openmc.YPlane(-1.0, boundary_type="vacuum")
|
||||
ymax = openmc.YPlane(1.0, boundary_type="vacuum")
|
||||
zmin = openmc.ZPlane(-1.0, boundary_type="vacuum")
|
||||
zmax = openmc.ZPlane(1.0, boundary_type="vacuum")
|
||||
|
||||
inner_cell = openmc.Cell()
|
||||
inner_univ = openmc.Universe(cells=[inner_cell])
|
||||
|
||||
tile_cell = openmc.Cell(fill=inner_univ)
|
||||
tile_univ = openmc.Universe(cells=[tile_cell])
|
||||
|
||||
lattice = openmc.RectLattice()
|
||||
lattice.lower_left = (-1.0, -1.0)
|
||||
lattice.pitch = (1.0, 2.0)
|
||||
lattice.universes = [[tile_univ, tile_univ]]
|
||||
|
||||
root_cell = openmc.Cell(
|
||||
fill=lattice, region=+xmin & -xmax & +ymin & -ymax & +zmin & -zmax)
|
||||
model.geometry = openmc.Geometry([root_cell])
|
||||
|
||||
src = openmc.IndependentSource()
|
||||
src.space = openmc.stats.Point((-0.5, 0.0, 0.0))
|
||||
src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0))
|
||||
|
||||
model.settings.run_mode = 'fixed source'
|
||||
model.settings.batches = 1
|
||||
model.settings.particles = 5
|
||||
model.settings.source = src
|
||||
|
||||
current_tally = openmc.Tally()
|
||||
current_tally.filters = [openmc.SurfaceFilter(xmax)]
|
||||
current_tally.scores = ['current']
|
||||
model.tallies = [current_tally]
|
||||
|
||||
model.run(apply_tally_results=True)
|
||||
assert current_tally.mean.flat[0] == pytest.approx(1.0)
|
||||
|
||||
|
||||
def test_cellfrom_filter_flux_directional(two_cell_model, run_in_tmpdir):
|
||||
"""SurfaceFilter + CellFromFilter + flux scores only the correct direction."""
|
||||
model, xmid, cell1, cell2 = two_cell_model
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue