crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 21:25:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'openmc-dev:develop' into convert-to-multigroup-material-cell-wise

Browse source
This commit is contained in:
Jonathan Shimwell 2026-07-21 11:42:16 +02:00 • committed by GitHub
commit 3bc19ef017
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
90 changed files with 3586 additions and 451 deletions
Download patch file Download diff file Expand all files Collapse all files

1
.gitignore vendored
View file

@ -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/

2
.readthedocs.yaml
View file

@ -7,6 +7,8 @@ build:
jobs:
post_checkout:
- git fetch --unshallow || true
- cd docs/doxygen && doxygen && cd -
sphinx:
configuration: docs/source/conf.py

2
CMakeLists.txt
View file

@ -532,7 +532,7 @@ else()
endif()
if(OPENMC_USE_DAGMC)
target_compile_definitions(libopenmc PRIVATE OPENMC_DAGMC_ENABLED)
target_compile_definitions(libopenmc PUBLIC OPENMC_DAGMC_ENABLED)
target_link_libraries(libopenmc dagmc-shared)
if(OPENMC_USE_UWUW)

1
docs/Makefile
View file

@ -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
View 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

43
docs/source/capi/index.rst
View file

@ -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)

13
docs/source/conf.py
View file

@ -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 ---------------------------------------------------

3
docs/source/devguide/contributing.rst
View file

@ -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.

5
docs/source/devguide/docbuild.rst
View file

@ -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
-----------------------------------

6
docs/source/io_formats/geometry.rst
View file

@ -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

80
docs/source/io_formats/settings.rst
View file

@ -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).

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

@ -26,6 +26,7 @@ Simulation Settings
openmc.FileSource
openmc.CompiledSource
openmc.MeshSource
openmc.TokamakSource
openmc.SourceParticle
openmc.VolumeCalculation
openmc.Settings

24
docs/source/usersguide/decay_sources.rst
View file

@ -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)

13
docs/source/usersguide/random_ray.rst
View file

@ -1116,10 +1116,15 @@ external source is present in the problem. Simulation settings (e.g., number of
rays, batches, etc.) will be identical for both calculations. At the
conclusion of the run, all results (e.g., tallies, plots, etc.) will be
derived from the adjoint flux rather than the forward flux but are not labeled
any differently. The initial forward flux solution will not be stored or
available in the final statepoint file. Those wishing to do analysis requiring
both the forward and adjoint solutions will need to run two separate
simulations and load both statepoint files.
any differently. When an initial forward solve is performed (i.e., when no
user-specified adjoint source is present), its output files are also written to
disk with a ``forward`` infix, so they are not overwritten by the subsequent
adjoint solve. This applies to the statepoint, ``tallies.out``, and any voxel
plots, e.g., ``statepoint.forward.N.h5`` and ``tallies.forward.out``; the
adjoint solve keeps the usual file names. This allows analyses requiring both
the forward and adjoint solutions to be performed from a single run. When
generating FW-CADIS weight windows, no weight window file is written for the
forward solve, as only the final adjoint-derived weight windows are meaningful.
.. note::
Use of the automated

47
docs/source/usersguide/settings.rst
View file

@ -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
------------------

29
include/openmc/capi.h
View file

@ -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(

1
include/openmc/constants.h
View file

@ -368,6 +368,7 @@ enum class SolverType { MONTE_CARLO, RANDOM_RAY };
enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID };
enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY };
enum class RandomRaySampleMethod { PRNG, HALTON, S2 };
enum class RandomRaySolve { FORWARD, FORWARD_FOR_ADJOINT, ADJOINT };
//==============================================================================
// Geometry Constants

4
include/openmc/error.h
View file

@ -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

38
include/openmc/geometry.h
View file

@ -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

26
include/openmc/math_functions.h
View file

@ -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

6
include/openmc/particle_data.h
View file

@ -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

14
include/openmc/plot.h
View file

@ -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
}

5
include/openmc/random_ray/flat_source_domain.h
View file

@ -76,7 +76,10 @@ public:
//----------------------------------------------------------------------------
// Static Data members
static bool volume_normalized_flux_tallies_;
static bool adjoint_; // If the user wants outputs based on the adjoint flux
// If the user wants outputs based on the adjoint flux
static bool adjoint_requested_;
// The solve currently being executed
static RandomRaySolve solve_;
static bool fw_cadis_local_;
static double
diagonal_stabilization_rho_; // Adjusts strength of diagonal stabilization

2
include/openmc/random_ray/random_ray_simulation.h
View file

@ -22,7 +22,7 @@ public:
void apply_fixed_sources_and_mesh_domains();
void prepare_fw_fixed_sources_adjoint();
void prepare_local_fixed_sources_adjoint();
void prepare_adjoint_simulation(bool fw_adjoint);
void prepare_adjoint_simulation(bool from_forward);
void simulate();
void output_simulation_results() const;
void instability_check(

6
include/openmc/random_ray/source_region.h
View file

@ -51,10 +51,10 @@ inline void hash_combine(size_t& seed, const size_t v)
// every iteration.
struct TallyTask {
int tally_idx;
int filter_idx;
int64_t filter_idx;
int score_idx;
int score_type;
TallyTask(int tally_idx, int filter_idx, int score_idx, int score_type)
TallyTask(int tally_idx, int64_t filter_idx, int score_idx, int score_type)
: tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx),
score_type(score_type)
{}
@ -690,7 +690,7 @@ private:
// Private Methods
// Helper function for indexing
inline int index(int64_t sr, int g) const { return sr * negroups_ + g; }
inline int64_t index(int64_t sr, int g) const { return sr * negroups_ + g; }
};
} // namespace openmc

26
include/openmc/shared_array.h
View file

@ -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_; }

136
include/openmc/source.h
View file

@ -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
//==============================================================================

8
include/openmc/tallies/tally.h
View file

@ -97,9 +97,9 @@ public:
//! Given already-set filters, set the stride lengths
void set_strides();
int32_t strides(int i) const { return strides_[i]; }
int64_t strides(int i) const { return strides_[i]; }
int32_t n_filter_bins() const { return n_filter_bins_; }
int64_t n_filter_bins() const { return n_filter_bins_; }
bool multiply_density() const { return multiply_density_; }
@ -184,9 +184,9 @@ private:
vector<int32_t> filters_; //!< Filter indices in global filters array
//! Index strides assigned to each filter to support 1D indexing.
vector<int32_t> strides_;
vector<int64_t> strides_;
int32_t n_filter_bins_ {0};
int64_t n_filter_bins_ {0};
//! Whether to multiply by atom density for reaction rates
bool multiply_density_ {true};

2
include/openmc/tallies/tally_scoring.h
View file

@ -41,7 +41,7 @@ public:
FilterBinIter& operator++();
int index_ {1};
int64_t index_ {1};
double weight_ {1.};
vector<FilterMatch>& filter_matches_;

8
include/openmc/weight_windows.h
View file

@ -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)

60
openmc/data/data.py
View file

@ -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

129
openmc/data/dose/dose.py
View file

@ -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
View 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
View 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

10
openmc/data/njoy.py
View file

@ -259,14 +259,14 @@ broadr / %%%%%%%%%%%%%%%%%%%%%%% Doppler broaden XS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
_TEMPLATE_HEATR = """
heatr / %%%%%%%%%%%%%%%%%%%%%%%%% Add heating kerma %%%%%%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {nheatr_in} {nheatr} /
{mat} 4 0 0 0 /
{mat} 4 0 0 0 0 {ed}/
302 318 402 444 /
"""
_TEMPLATE_HEATR_LOCAL = """
heatr / %%%%%%%%%%%%%%%%% Add heating kerma (local photons) %%%%%%%%%%%%%%%%%%%%
{nendf} {nheatr_in} {nheatr_local} /
{mat} 4 0 0 1 /
{mat} 4 0 0 1 0 {ed}/
302 318 402 444 /
"""
@ -411,7 +411,7 @@ def make_pendf(filename, pendf='pendf', **kwargs):
def make_ace(filename, temperatures=None, acer=True, xsdir=None,
output_dir=None, pendf=False, error=0.001, broadr=True,
heatr=True, gaspr=True, purr=True, evaluation=None,
smoothing=True, **kwargs):
smoothing=True, displacement_energy=None, **kwargs):
"""Generate incident neutron ACE file from an ENDF file
File names can be passed to
@ -463,6 +463,9 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
indicates which evaluation should be used.
smoothing : bool, optional
If the smoothing option (ACER card 6) is on (True) or off (False).
displacement_energy : float, optional
Threshold displacement energy in [eV]. Used in HEATR to calculate
damage-energy cross section. When None, use NJOY defaults.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.run`
@ -515,6 +518,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
# heatr
if heatr:
ed = displacement_energy if displacement_energy is not None else 0
nheatr_in = nlast
nheatr_local = nheatr_in + 1
tapeout[nheatr_local] = (output_dir / "heatr_local") if heatr is True \

3
openmc/deplete/microxs.py
View file

@ -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

185
openmc/deplete/r2s.py
View file

@ -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

35
openmc/deplete/results.py
View file

@ -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

56
openmc/lib/plot.py
View file

@ -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

32
openmc/material.py
View file

@ -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())

20
openmc/model/model.py
View file

@ -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))

14
openmc/plots.py
View file

@ -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
View file

@ -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:

1
pyproject.toml
View file

@ -40,6 +40,7 @@ dependencies = [
[project.optional-dependencies]
depletion-mpi = ["mpi4py"]
docs = [
"breathe",
"sphinx",
"sphinxcontrib-katex",
"sphinx-numfig",

10
src/boundary_condition.cpp
View file

@ -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);

11
src/distribution.cpp
View file

@ -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];

8
src/distribution_spatial.cpp
View file

@ -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

37
src/geometry.cpp
View file

@ -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;
}
//==============================================================================

33
src/lattice.cpp
View file

@ -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};
}

50
src/math_functions.cpp
View file

@ -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

22
src/mesh.cpp
View file

@ -6,7 +6,8 @@
#define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers
#include <cmath> // for ceil
#include <cstddef> // for size_t
#include <numeric> // for accumulate
#include <limits>
#include <numeric> // for accumulate
#include <string>
#ifdef _MSC_VER
@ -1080,13 +1081,26 @@ int StructuredMesh::get_bin(Position r) const
int StructuredMesh::n_bins() const
{
return std::accumulate(
shape_.begin(), shape_.begin() + n_dimension_, 1, std::multiplies<>());
// Bin indices are stored as 32-bit ints in the tally system.
int64_t n = 1;
for (int i = 0; i < n_dimension_; ++i)
n *= shape_[i];
if (n > std::numeric_limits<int>::max()) {
fatal_error(fmt::format(
"Mesh {} has too many bins ({}) for 32-bit tally indexing", id_, n));
}
return static_cast<int>(n);
}
int StructuredMesh::n_surface_bins() const
{
return 4 * n_dimension_ * n_bins();
// Surface bin indices are stored as 32-bit ints in the tally system.
int64_t n = static_cast<int64_t>(n_bins()) * 4 * n_dimension_;
if (n > std::numeric_limits<int>::max()) {
fatal_error(fmt::format(
"Mesh {} has too many surface bins ({}) for tally indexing", id_, n));
}
return static_cast<int>(n);
}
tensor::Tensor<double> StructuredMesh::count_sites(

9
src/output.cpp
View file

@ -621,8 +621,15 @@ void write_tallies()
if (model::tallies.empty())
return;
// Tag tallies.out written during the forward solve of an adjoint run
const char* forward =
(FlatSourceDomain::solve_ == RandomRaySolve::FORWARD_FOR_ADJOINT)
? "forward."
: "";
// Set filename for tallies_out
std::string filename = fmt::format("{}tallies.out", settings::path_output);
std::string filename =
fmt::format("{}tallies.{}out", settings::path_output, forward);
// Open the tallies.out file.
std::ofstream tallies_out;

28
src/particle.cpp
View file

@ -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;

4
src/physics.cpp
View file

@ -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
View file

@ -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);

72
src/random_lcg.cpp
View file

@ -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;
}
//==============================================================================

20
src/random_ray/flat_source_domain.cpp
View file

@ -30,7 +30,8 @@ namespace openmc {
RandomRayVolumeEstimator FlatSourceDomain::volume_estimator_ {
RandomRayVolumeEstimator::HYBRID};
bool FlatSourceDomain::volume_normalized_flux_tallies_ {false};
bool FlatSourceDomain::adjoint_ {false};
bool FlatSourceDomain::adjoint_requested_ {false};
RandomRaySolve FlatSourceDomain::solve_ {RandomRaySolve::FORWARD};
bool FlatSourceDomain::fw_cadis_local_ {false};
double FlatSourceDomain::diagonal_stabilization_rho_ {1.0};
std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>
@ -556,7 +557,7 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const
// If we are in adjoint mode of a fixed source problem, the external
// source is already normalized, such that all resulting fluxes are
// also normalized.
if (adjoint_) {
if (solve_ == RandomRaySolve::ADJOINT) {
return 1.0;
}
@ -724,7 +725,7 @@ void FlatSourceDomain::random_ray_tally()
for (int i = 0; i < model::tallies.size(); i++) {
Tally& tally {*model::tallies[i]};
#pragma omp parallel for
for (int bin = 0; bin < tally.n_filter_bins(); bin++) {
for (int64_t bin = 0; bin < tally.n_filter_bins(); bin++) {
for (int score_idx = 0; score_idx < tally.n_scores(); score_idx++) {
auto score_type = tally.scores_[score_idx];
if (score_type == SCORE_FLUX) {
@ -795,6 +796,12 @@ void FlatSourceDomain::output_to_vtk() const
double z_delta = width.z / Nz;
std::string filename = openmc_plot->path_plot();
// Tag plots written during the forward solve of an adjoint run
if (solve_ == RandomRaySolve::FORWARD_FOR_ADJOINT) {
auto dot = filename.find_last_of('.');
filename = filename.substr(0, dot) + ".forward" + filename.substr(dot);
}
// Perform sanity checks on file size
uint64_t bytes = Nx * Ny * Nz * (negroups_ + 1 + 1 + 1) * sizeof(float);
write_message(5, "Processing plot {}: {}... (Estimated size is {} MB)",
@ -1002,9 +1009,10 @@ void FlatSourceDomain::output_to_vtk() const
void FlatSourceDomain::apply_external_source_to_source_region(
int src_idx, SourceRegionHandle& srh)
{
auto s = (adjoint_ && !model::adjoint_sources.empty())
? model::adjoint_sources[src_idx].get()
: model::external_sources[src_idx].get();
auto s =
(solve_ == RandomRaySolve::ADJOINT && !model::adjoint_sources.empty())
? model::adjoint_sources[src_idx].get()
: model::external_sources[src_idx].get();
auto is = dynamic_cast<IndependentSource*>(s);
auto discrete = dynamic_cast<Discrete*>(is->energy());
double strength_factor = is->strength();

87
src/random_ray/random_ray_simulation.cpp
View file

@ -189,7 +189,7 @@ void validate_random_ray_inputs()
// Validate adjoint sources
///////////////////////////////////////////////////////////////////
if (FlatSourceDomain::adjoint_ && !model::adjoint_sources.empty()) {
if (FlatSourceDomain::adjoint_requested_ && !model::adjoint_sources.empty()) {
for (int i = 0; i < model::adjoint_sources.size(); i++) {
Source* s = model::adjoint_sources[i].get();
@ -289,7 +289,8 @@ void openmc_finalize_random_ray()
{
FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::HYBRID;
FlatSourceDomain::volume_normalized_flux_tallies_ = false;
FlatSourceDomain::adjoint_ = false;
FlatSourceDomain::adjoint_requested_ = false;
FlatSourceDomain::solve_ = RandomRaySolve::FORWARD;
FlatSourceDomain::fw_cadis_local_ = false;
FlatSourceDomain::fw_cadis_local_targets_.clear();
FlatSourceDomain::mesh_domain_map_.clear();
@ -356,19 +357,17 @@ void RandomRaySimulation::prepare_local_fixed_sources_adjoint()
}
}
void RandomRaySimulation::prepare_adjoint_simulation(bool fw_adjoint)
void RandomRaySimulation::prepare_adjoint_simulation(bool from_forward)
{
reset_timers();
if (mpi::master)
header("ADJOINT FLUX SOLVE", 3);
if (fw_adjoint) {
// Forward simulation has already been run;
// Configure the domain for adjoint simulation and
// re-initialize OpenMC general data structures
FlatSourceDomain::adjoint_ = true;
if (from_forward) {
// The forward solve has already run. Re-initialize OpenMC's general data
// structures for the adjoint solve and derive the adjoint source from the
// forward flux.
openmc_simulation_init();
prepare_fw_fixed_sources_adjoint();
@ -603,7 +602,8 @@ void RandomRaySimulation::print_results_random_ray(
}
fmt::print(" Volume Estimator Type = {}\n", estimator);
std::string adjoint_true = (FlatSourceDomain::adjoint_) ? "ON" : "OFF";
std::string adjoint_true =
(FlatSourceDomain::solve_ == RandomRaySolve::ADJOINT) ? "ON" : "OFF";
fmt::print(" Adjoint Flux Mode = {}\n", adjoint_true);
std::string shape;
@ -675,60 +675,49 @@ void RandomRaySimulation::print_results_random_ray(
void openmc_run_random_ray()
{
//////////////////////////////////////////////////////////
// Run forward simulation
//////////////////////////////////////////////////////////
using namespace openmc;
// Check if adjoint calculation is needed, and if local adjoint source(s)
// are present. If an adjoint calculation is needed and no sources are
// specified, we will run a forward calculation first to calculate adjoint
// sources for global variance reduction, then perform an adjoint
// calculation later.
bool adjoint_needed = openmc::FlatSourceDomain::adjoint_;
bool fw_adjoint = openmc::model::adjoint_sources.empty() && adjoint_needed;
// Determine which solves to run. If adjoint results are requested and no
// user-defined adjoint source is present, an initial forward solve is needed
// to construct the adjoint source from the forward flux (FW-CADIS). If the
// user has defined an adjoint source, the forward solve is skipped and only
// the adjoint solve is run.
const bool run_adjoint = FlatSourceDomain::adjoint_requested_;
const bool have_adjoint_source = !model::adjoint_sources.empty();
const bool run_forward = !(run_adjoint && have_adjoint_source);
// If we're going to do an adjoint simulation with forward-weighted adjoint
// sources afterwards, report that this is the initial forward flux solve.
if (!adjoint_needed || fw_adjoint) {
// Configure the domain for forward simulation
openmc::FlatSourceDomain::adjoint_ = false;
if (adjoint_needed && openmc::mpi::master)
openmc::header("FORWARD FLUX SOLVE", 3);
// Set the initial solve type
if (!run_forward) {
FlatSourceDomain::solve_ = RandomRaySolve::ADJOINT;
} else if (run_adjoint) {
FlatSourceDomain::solve_ = RandomRaySolve::FORWARD_FOR_ADJOINT;
} else {
// Configure domain for adjoint simulation (later)
openmc::FlatSourceDomain::adjoint_ = true;
FlatSourceDomain::solve_ = RandomRaySolve::FORWARD;
}
// Initialize OpenMC general data structures
openmc_simulation_init();
// Validate that inputs meet requirements for random ray mode
if (openmc::mpi::master)
openmc::validate_random_ray_inputs();
if (mpi::master)
validate_random_ray_inputs();
// Initialize Random Ray Simulation Object
openmc::RandomRaySimulation sim;
RandomRaySimulation sim;
if (!adjoint_needed || fw_adjoint) {
// Initialize fixed sources, if present
// Run the forward solve
if (run_forward) {
// When an adjoint solve follows, report this as the initial forward solve
if (run_adjoint && mpi::master)
header("FORWARD FLUX SOLVE", 3);
sim.apply_fixed_sources_and_mesh_domains();
// Execute random ray simulation
sim.simulate();
}
//////////////////////////////////////////////////////////
// Run adjoint simulation (if enabled)
//////////////////////////////////////////////////////////
if (!adjoint_needed) {
return;
// Run the adjoint solve
if (run_adjoint) {
FlatSourceDomain::solve_ = RandomRaySolve::ADJOINT;
sim.prepare_adjoint_simulation(run_forward);
sim.simulate();
}
// Setup for adjoint simulation
sim.prepare_adjoint_simulation(fw_adjoint);
// Execute random ray simulation
sim.simulate();
}

3
src/settings.cpp
View file

@ -322,7 +322,7 @@ void get_run_parameters(pugi::xml_node node_base)
get_node_value_bool(random_ray_node, "volume_normalized_flux_tallies");
}
if (check_for_node(random_ray_node, "adjoint")) {
FlatSourceDomain::adjoint_ =
FlatSourceDomain::adjoint_requested_ =
get_node_value_bool(random_ray_node, "adjoint");
}
if (check_for_node(random_ray_node, "sample_method")) {
@ -1247,6 +1247,7 @@ void read_settings_xml(pugi::xml_node root)
// read weight windows from file
if (check_for_node(root, "weight_windows_file")) {
weight_windows_file = get_node_value(root, "weight_windows_file");
weight_windows_on = true;
}
// read settings for weight windows value, this will override

169
src/simulation.cpp
View file

@ -12,10 +12,12 @@
#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"
#include "openmc/random_lcg.h"
#include "openmc/random_ray/flat_source_domain.h"
#include "openmc/settings.h"
#include "openmc/source.h"
#include "openmc/state_point.h"
@ -40,6 +42,7 @@
#include <algorithm>
#include <cmath>
#include <numeric>
#include <string>
//==============================================================================
@ -200,9 +203,12 @@ int openmc_simulation_finalize()
if (settings::output_tallies && mpi::master)
write_tallies();
// If weight window generators are present in this simulation,
// write a weight windows file
if (variance_reduction::weight_windows_generators.size() > 0) {
// If weight window generators are present in this simulation, write a
// weight windows file. This is skipped during the forward solve of an
// adjoint (FW-CADIS) run, where only the adjoint-derived weight windows
// are meaningful.
if (variance_reduction::weight_windows_generators.size() > 0 &&
FlatSourceDomain::solve_ != RandomRaySolve::FORWARD_FOR_ADJOINT) {
openmc_weight_windows_export();
}
@ -350,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
//==============================================================================
@ -882,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);
}
}
}
@ -914,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;
@ -947,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
@ -979,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);
@ -996,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
@ -1071,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;
@ -1141,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;

483
src/source.cpp
View file

@ -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));
}
@ -523,9 +530,16 @@ void FileSource::load_sites_from_file(const std::string& path)
file_close(file_id);
}
// Make sure particles in source file have valid types
// Make sure particles in source file have valid types. If any particle is a
// photon, electron, or positron, enable photon transport so that the
// appropriate cross sections are loaded.
for (const auto& site : this->sites_) {
validate_particle_type(site.particle, "FileSource");
if (site.particle == ParticleType::photon() ||
site.particle == ParticleType::electron() ||
site.particle == ParticleType::positron()) {
settings::photon_transport = true;
}
}
}
@ -666,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
//==============================================================================

11
src/state_point.cpp
View file

@ -22,6 +22,7 @@
#include "openmc/nuclide.h"
#include "openmc/output.h"
#include "openmc/particle_type.h"
#include "openmc/random_ray/flat_source_domain.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/derivative.h"
@ -46,9 +47,15 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
// Determine width for zero padding
int w = std::to_string(settings::n_max_batches).size();
// Tag statepoints written during the forward solve of an adjoint run
const char* forward =
(FlatSourceDomain::solve_ == RandomRaySolve::FORWARD_FOR_ADJOINT)
? "forward."
: "";
// Set filename for state point
filename_ = fmt::format("{0}statepoint.{1:0{2}}.h5", settings::path_output,
simulation::current_batch, w);
filename_ = fmt::format("{0}statepoint.{3}{1:0{2}}.h5",
settings::path_output, simulation::current_batch, w, forward);
}
// If a file name was specified, ensure it has .h5 file extension

6
src/tallies/tally.cpp
View file

@ -512,7 +512,7 @@ void Tally::set_strides()
// longest stride.
auto n = filters_.size();
strides_.resize(n, 0);
int stride = 1;
int64_t stride = 1;
for (int i = n - 1; i >= 0; --i) {
strides_[i] = stride;
stride *= model::tally_filters[filters_[i]]->n_bins();
@ -887,7 +887,7 @@ void Tally::accumulate()
if (higher_moments_) {
#pragma omp parallel for
// filter bins (specific cell, energy bins)
for (int i = 0; i < results_.shape(0); ++i) {
for (int64_t i = 0; i < results_.shape(0); ++i) {
// score bins (flux, total reaction rate, fission reaction rate, etc.)
for (int j = 0; j < results_.shape(1); ++j) {
double val = results_(i, j, TallyResult::VALUE) * norm;
@ -902,7 +902,7 @@ void Tally::accumulate()
} else {
#pragma omp parallel for
// filter bins (specific cell, energy bins)
for (int i = 0; i < results_.shape(0); ++i) {
for (int64_t i = 0; i < results_.shape(0); ++i) {
// score bins (flux, total reaction rate, fission reaction rate, etc.)
for (int j = 0; j < results_.shape(1); ++j) {
double val = results_(i, j, TallyResult::VALUE) * norm;

52
src/tallies/tally_scoring.cpp
View file

@ -158,7 +158,7 @@ void score_fission_delayed_dg(int i_tally, int d_bin, double score,
dg_match.bins_[i_bin] = d_bin;
// Determine the filter scoring index
auto filter_index = 0;
int64_t filter_index = 0;
double filter_weight = 1.;
for (auto i = 0; i < tally.filters().size(); ++i) {
auto i_filt = tally.filters(i);
@ -449,7 +449,7 @@ void score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin)
(score_bin == SCORE_PROMPT_NU_FISSION && g == 0)) {
// Find the filter scoring index for this filter combination
int filter_index = 0;
int64_t filter_index = 0;
double filter_weight = 1.0;
for (auto j = 0; j < tally.filters().size(); ++j) {
auto i_filt = tally.filters(j);
@ -497,7 +497,7 @@ void score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin)
} else {
// Find the filter index and weight for this filter combination
int filter_index = 0;
int64_t filter_index = 0;
double filter_weight = 1.;
for (auto j = 0; j < tally.filters().size(); ++j) {
auto i_filt = tally.filters(j);
@ -578,8 +578,8 @@ double get_nuclide_xs(const Particle& p, int i_nuclide, int score_bin)
//! collision estimator.
void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
int filter_index, double filter_weight, int i_nuclide, double atom_density,
double flux)
int64_t filter_index, double filter_weight, int i_nuclide,
double atom_density, double flux)
{
Tally& tally {*model::tallies[i_tally]};
@ -1112,8 +1112,8 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
//! is not used for analog tallies.
void score_general_ce_analog(Particle& p, int i_tally, int start_index,
int filter_index, double filter_weight, int i_nuclide, double atom_density,
double flux)
int64_t filter_index, double filter_weight, int i_nuclide,
double atom_density, double flux)
{
Tally& tally {*model::tallies[i_tally]};
@ -1615,8 +1615,8 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
//! argument is really just used for filter weights.
void score_general_mg(Particle& p, int i_tally, int start_index,
int filter_index, double filter_weight, int i_nuclide, double atom_density,
double flux)
int64_t filter_index, double filter_weight, int i_nuclide,
double atom_density, double flux)
{
auto& tally {*model::tallies[i_tally]};
@ -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

4
src/tallies/trigger.cpp
View file

@ -28,7 +28,7 @@ KTrigger keff_trigger;
//==============================================================================
std::pair<double, double> get_tally_uncertainty(
int i_tally, int score_index, int filter_index)
int i_tally, int score_index, int64_t filter_index)
{
const auto& tally {model::tallies[i_tally]};
@ -71,7 +71,7 @@ void check_tally_triggers(double& ratio, int& tally_id, int& score)
continue;
const auto& results = t.results_;
for (auto filter_index = 0; filter_index < results.shape(0);
for (int64_t filter_index = 0; filter_index < results.shape(0);
++filter_index) {
// Compute the tally uncertainty metrics.
auto uncert_pair =

7
src/weight_windows.cpp
View file

@ -791,7 +791,7 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
if (method_string == "magic") {
method_ = WeightWindowUpdateMethod::MAGIC;
if (settings::solver_type == SolverType::RANDOM_RAY &&
FlatSourceDomain::adjoint_) {
FlatSourceDomain::adjoint_requested_) {
fatal_error("Random ray weight window generation with MAGIC cannot be "
"done in adjoint mode.");
}
@ -800,7 +800,7 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
if (settings::solver_type != SolverType::RANDOM_RAY) {
fatal_error("FW-CADIS can only be run in random ray solver mode.");
}
FlatSourceDomain::adjoint_ = true;
FlatSourceDomain::adjoint_requested_ = true;
if (check_for_node(node, "targets")) {
FlatSourceDomain::fw_cadis_local_ = true;
targets_ = get_node_array<size_t>(node, "targets");
@ -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();

1
tests/cpp_unit_tests/CMakeLists.txt
View file

@ -6,6 +6,7 @@ set(TEST_NAMES
test_math
test_mcpl_stat_sum
test_mesh
test_ray
test_region
test_tensor
# Add additional unit test files here

25
tests/cpp_unit_tests/test_distribution.cpp
View file

@ -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};

36
tests/cpp_unit_tests/test_math.cpp
View file

@ -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));
}

117
tests/cpp_unit_tests/test_ray.cpp Normal file
View file

@ -0,0 +1,117 @@
#include <cmath>
#include <utility>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <pugixml.hpp>
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/ray.h"
#include "openmc/settings.h"
#include "openmc/surface.h"
namespace {
using Catch::Matchers::WithinAbs;
constexpr double DISTANCE_TOLERANCE = 1.0e-12;
class SphericalShellFixture {
public:
SphericalShellFixture()
: run_mode_(openmc::settings::run_mode),
root_universe_(openmc::model::root_universe),
n_coord_levels_(openmc::model::n_coord_levels)
{
openmc::settings::run_mode = openmc::RunMode::PLOTTING;
constexpr auto geometry = R"(
<geometry>
<surface id="1" type="sphere" coeffs="0 0 0 1"/>
<surface id="2" type="sphere" coeffs="0 0 0 2"
boundary="vacuum"/>
<cell id="1" material="void" region="-1"/>
<cell id="2" material="void" region="1 -2"/>
</geometry>
)";
auto result = document_.load_string(geometry);
REQUIRE(result);
openmc::read_geometry_xml(document_.document_element());
openmc::finalize_geometry();
openmc::finalize_cell_densities();
}
~SphericalShellFixture()
{
openmc::free_memory_geometry();
openmc::free_memory_surfaces();
openmc::model::universe_level_counts.clear();
openmc::model::root_universe = root_universe_;
openmc::model::n_coord_levels = n_coord_levels_;
openmc::settings::run_mode = run_mode_;
}
private:
pugi::xml_document document_;
openmc::RunMode run_mode_;
int root_universe_;
int n_coord_levels_;
};
class RecordingRay : public openmc::Ray {
public:
using Ray::Ray;
void on_intersection() override
{
intersections.emplace_back(
traversal_distance_, boundary().surface_index() + 1);
}
openmc::vector<std::pair<double, int>> intersections;
};
} // namespace
TEST_CASE_METHOD(SphericalShellFixture, "Trace a single chord through a sphere")
{
constexpr double impact_parameter = 1.5;
const double half_chord =
std::sqrt(4.0 - impact_parameter * impact_parameter);
RecordingRay ray({-3.0, impact_parameter, 0.0}, {1.0, 0.0, 0.0});
ray.trace();
REQUIRE(ray.intersections.size() == 2);
CHECK(ray.intersections[0].second == 2);
CHECK_THAT(ray.intersections[0].first, WithinAbs(0.0, DISTANCE_TOLERANCE));
CHECK(ray.intersections[1].second == 2);
CHECK_THAT(ray.intersections[1].first,
WithinAbs(2.0 * half_chord, DISTANCE_TOLERANCE));
}
TEST_CASE_METHOD(
SphericalShellFixture, "Trace both segments of a spherical shell")
{
constexpr double impact_parameter = 0.5;
const double outer = std::sqrt(4.0 - impact_parameter * impact_parameter);
const double inner = std::sqrt(1.0 - impact_parameter * impact_parameter);
RecordingRay ray({-3.0, impact_parameter, 0.0}, {1.0, 0.0, 0.0});
ray.trace();
REQUIRE(ray.intersections.size() == 4);
CHECK(ray.intersections[0].second == 2);
CHECK_THAT(ray.intersections[0].first, WithinAbs(0.0, DISTANCE_TOLERANCE));
CHECK(ray.intersections[1].second == 1);
CHECK_THAT(
ray.intersections[1].first, WithinAbs(outer - inner, DISTANCE_TOLERANCE));
CHECK(ray.intersections[2].second == 1);
CHECK_THAT(
ray.intersections[2].first, WithinAbs(outer + inner, DISTANCE_TOLERANCE));
CHECK(ray.intersections[3].second == 2);
CHECK_THAT(
ray.intersections[3].first, WithinAbs(2.0 * outer, DISTANCE_TOLERANCE));
}

27
tests/cpp_unit_tests/test_tally.cpp
View file

@ -1,3 +1,4 @@
#include "openmc/tallies/filter_energy.h"
#include "openmc/tallies/tally.h"
#include <catch2/catch_test_macros.hpp>
@ -46,10 +47,34 @@ TEST_CASE("Test add/set_filter")
REQUIRE(model::filter_map[cell_filter->id()] == tally->filters(0));
REQUIRE(model::filter_map[particle_filter->id()] == tally->filters(1));
// set filters with a duplicate filter, should only add the filter to the tally once
// set filters with a duplicate filter, should only add the filter to the
// tally once
filters = {cell_filter, cell_filter};
tally->set_filters(filters);
REQUIRE(tally->filters().size() == 1);
REQUIRE(model::filter_map[cell_filter->id()] == tally->filters(0));
}
// Regression test for 64-bit tally filter-bin counts (mesh x groups > 2^31).
TEST_CASE("Tally filter-bin count does not overflow 32 bits")
{
// Two energy filters whose bin counts multiply to 2.5e9, above INT32_MAX.
constexpr int64_t bins_per_filter = 50000;
// Only the bin count matters here, so the edge values are an arbitrary ramp.
std::vector<double> edges(bins_per_filter + 1);
for (int64_t i = 0; i < bins_per_filter + 1; ++i)
edges[i] = static_cast<double>(i);
Tally* tally = Tally::create();
for (int i = 0; i < 2; ++i) {
Filter* filter = Filter::create("energy");
dynamic_cast<EnergyFilter*>(filter)->set_bins(edges);
tally->add_filter(filter);
}
tally->set_strides();
// set_strides() previously accumulated this product in a 32-bit int.
REQUIRE(tally->n_filter_bins() == bins_per_filter * bins_per_filter);
REQUIRE(tally->n_filter_bins() > 2147483647);
}

0
tests/regression_tests/lattice_large_pitch/__init__.py Normal file
View file

51
tests/regression_tests/lattice_large_pitch/inputs_true.dat Normal file
View 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
View 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
View 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()

55
tests/regression_tests/weightwindows/test.py
View file

@ -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]

3
tests/testing_harness.py
View file

@ -143,7 +143,8 @@ class TestHarness:
def _cleanup(self):
"""Delete statepoints, tally, and test files."""
output = glob.glob('statepoint.*.h5')
output += ['tallies.out', 'results_test.dat', 'summary.h5']
output += ['tallies.out', 'tallies.forward.out']
output += ['results_test.dat', 'summary.h5']
output += glob.glob('volume_*.h5')
for f in output:
if os.path.exists(f):

22
tests/unit_tests/test_data_decay.py
View file

@ -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)

27
tests/unit_tests/test_data_dose.py
View file

@ -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

33
tests/unit_tests/test_data_misc.py
View file

@ -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

31
tests/unit_tests/test_deplete_resultslist.py
View file

@ -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")

16
tests/unit_tests/test_lib.py
View file

@ -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)

31
tests/unit_tests/test_material.py
View file

@ -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'

23
tests/unit_tests/test_model.py
View file

@ -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
View 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])

120
tests/unit_tests/test_r2s.py
View file

@ -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
View 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

24
tests/unit_tests/test_source_file.py
View file

@ -145,3 +145,27 @@ def test_source_file_transport(run_in_tmpdir):
# Try running OpenMC
model.run()
def test_source_file_photon_transport(run_in_tmpdir):
# Create a source file containing a photon. Note that photon_transport is
# not explicitly enabled in the settings -- it should be turned on
# automatically because the source file contains a photon.
particle = openmc.SourceParticle(E=1.0e6, particle='photon')
openmc.write_source_file([particle], 'photon_source.h5')
# Create simple model to use the photon source file
model = openmc.Model()
al = openmc.Material()
al.add_element('Al', 1.0)
al.set_density('g/cm3', 2.7)
sph = openmc.Sphere(r=10.0, boundary_type='vacuum')
cell = openmc.Cell(fill=al, region=-sph)
model.geometry = openmc.Geometry([cell])
model.settings.source = openmc.FileSource(path='photon_source.h5')
model.settings.particles = 10
model.settings.batches = 3
model.settings.run_mode = 'fixed source'
# Running OpenMC should succeed
model.run()

263
tests/unit_tests/test_source_tokamak.py Normal file
View 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)

44
tests/unit_tests/test_surface_flux.py
View file

@ -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