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

Merge branch 'develop' into xs-change

Browse source
This commit is contained in:
David Foster 2026-06-25 11:37:25 +01:00 • committed by GitHub
commit 670fee26cb
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
166 changed files with 6490 additions and 1418 deletions
Download patch file Download diff file Expand all files Collapse all files

2
.github/workflows/ci.yml vendored
View file

@ -177,7 +177,7 @@ jobs:
- name: Setup tmate debug session
continue-on-error: true
if: ${{ contains(env.COMMIT_MESSAGE, '[gha-debug]') }}
if: ${{ failure() && contains(env.COMMIT_MESSAGE, '[gha-debug]') }}
uses: mxschmitt/action-tmate@v3
timeout-minutes: 10

24
CMakeLists.txt
View file

@ -611,9 +611,7 @@ add_custom_command(TARGET libopenmc POST_BUILD
#===============================================================================
# Install executable, scripts, manpage, license
#===============================================================================
configure_file(cmake/OpenMCConfig.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" @ONLY)
configure_file(cmake/OpenMCConfigVersion.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake" @ONLY)
include(CMakePackageConfigHelpers)
set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC)
install(TARGETS openmc libopenmc
@ -627,10 +625,24 @@ install(EXPORT openmc-targets
NAMESPACE OpenMC::
DESTINATION ${INSTALL_CONFIGDIR})
configure_package_config_file(
"cmake/OpenMCConfig.cmake.in"
"${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake"
INSTALL_DESTINATION ${INSTALL_CONFIGDIR}
)
write_basic_package_version_file(
"${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake"
VERSION ${OPENMC_VERSION}
COMPATIBILITY AnyNewerVersion
)
install(FILES
"${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake"
"${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake"
DESTINATION ${INSTALL_CONFIGDIR})
"${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake"
"${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake"
DESTINATION "${INSTALL_CONFIGDIR}"
)
install(FILES man/man1/openmc.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1)
install(FILES LICENSE DESTINATION "${CMAKE_INSTALL_DOCDIR}" RENAME copyright)
install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})

34
cmake/OpenMCConfig.cmake.in
View file

@ -1,12 +1,18 @@
get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY)
@PACKAGE_INIT@
# Compute the install prefix from this file's location
get_filename_component(_OPENMC_PREFIX "${OpenMC_CMAKE_DIR}/../../.." ABSOLUTE)
include("${CMAKE_CURRENT_LIST_DIR}/OpenMCConfigVersion.cmake")
include(CMakeFindDependencyMacro)
# Explicitly calculate prefix if it was not generated above
if(NOT DEFINED PACKAGE_PREFIX_DIR)
get_filename_component(PACKAGE_PREFIX_DIR "${CMAKE_CURRENT_LIST_DIR}/../../.." ABSOLUTE)
endif()
find_dependency(fmt CONFIG REQUIRED HINTS ${PACKAGE_PREFIX_DIR})
find_dependency(pugixml CONFIG REQUIRED HINTS ${PACKAGE_PREFIX_DIR})
find_package(fmt CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
find_package(pugixml CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
if(@OPENMC_USE_DAGMC@)
find_package(DAGMC REQUIRED HINTS @DAGMC_DIR@)
find_dependency(DAGMC REQUIRED HINTS @DAGMC_DIR@)
endif()
if(@OPENMC_USE_LIBMESH@)
@ -16,20 +22,24 @@ if(@OPENMC_USE_LIBMESH@)
pkg_check_modules(LIBMESH REQUIRED @LIBMESH_PC_FILE@>=1.7.0 IMPORTED_TARGET)
endif()
find_package(PNG)
if(NOT TARGET OpenMC::libopenmc)
include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake")
if("@PNG_FOUND@")
find_dependency(PNG)
endif()
if(@OPENMC_USE_MPI@)
find_package(MPI REQUIRED)
find_dependency(MPI REQUIRED)
endif()
if(@OPENMC_USE_OPENMP@)
find_package(OpenMP REQUIRED)
find_dependency(OpenMP REQUIRED)
endif()
if(@OPENMC_USE_UWUW@ AND NOT ${DAGMC_BUILD_UWUW})
message(FATAL_ERROR "UWUW is enabled in OpenMC but the DAGMC installation discovered was not configured with UWUW.")
endif()
include("${CMAKE_CURRENT_LIST_DIR}/OpenMCTargets.cmake")
if(NOT OpenMC_FIND_QUIETLY)
message(STATUS "Found OpenMC: ${PACKAGE_VERSION} (found in ${PACKAGE_PREFIX_DIR})")
endif()

11
cmake/OpenMCConfigVersion.cmake.in
View file

@ -1,11 +0,0 @@
set(PACKAGE_VERSION "@OPENMC_VERSION@")
# Check whether the requested PACKAGE_FIND_VERSION is compatible
if("${PACKAGE_VERSION}" VERSION_LESS "${PACKAGE_FIND_VERSION}")
set(PACKAGE_VERSION_COMPATIBLE FALSE)
else()
set(PACKAGE_VERSION_COMPATIBLE TRUE)
if ("${PACKAGE_VERSION}" VERSION_EQUAL "${PACKAGE_FIND_VERSION}")
set(PACKAGE_VERSION_EXACT TRUE)
endif()
endif()

4
docs/source/io_formats/collision_track.rst
View file

@ -10,7 +10,7 @@ may also be written after each batch when multiple files are requested
(``collision_track.N.h5``) or when the run is performed in parallel. The file
contains the information needed to reconstruct each recorded collision.
The current revision of the collision track file format is 1.1.
The current revision of the collision track file format is 1.2.
**/**
@ -33,7 +33,7 @@ The current revision of the collision track file format is 1.1.
- ``event_mt`` (*int*) -- ENDF MT number identifying the reaction.
- ``delayed_group`` (*int*) -- Delayed neutron group index (non-zero for delayed events).
- ``cell_id`` (*int*) -- ID of the cell in which the collision occurred.
- ``nuclide_id`` (*int*) -- ZA identifier of the nuclide (ZZZAAAM format).
- ``nuclide_id`` (*int*) -- PDG number of the nuclide (100ZZZAAAM).
- ``material_id`` (*int*) -- ID of the material containing the collision site.
- ``universe_id`` (*int*) -- ID of the universe containing the collision site.
- ``n_collision`` (*int*) -- Collision counter for the particle history.

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

@ -406,24 +406,55 @@ Each ``<dagmc_universe>`` element can have the following attributes or sub-eleme
*Default*: None
:material_overrides:
This element contains information on material overrides to be applied to the
DAGMC universe. It has the following attributes and sub-elements:
:cell:
Zero or more ``<cell>`` sub-elements may appear to override properties of
individual DAGMC volumes. Each ``<cell>`` element supports the following
attributes and sub-elements:
:cell:
Material override information for a single cell. It contains the following
attributes and sub-elements:
:id:
The integer cell ID in the DAGMC geometry to override. Required.
:id:
The cell ID in the DAGMC geometry for which the material override will
apply.
:name:
An optional string label for the cell.
:materials:
A list of material IDs that will apply to instances of the cell. If the
list contains only one ID, it will replace the original material
assignment of all instances of the DAGMC cell. If the list contains more
than one material, each material ID of the list will be assigned to the
various instances of the DAGMC cell.
*Default*: None
:material:
The material ID to assign to this cell. Use ``void`` for vacuum. Multiple
space-separated IDs may be given to specify a distribmat (distributed
material) assignment. Required.
:temperature:
Temperature(s) in [K] to assign to the cell. Must be ≥ 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
material fill. Multiple space-separated values may be given.
*Default*: None
:volume:
Volume of the cell in [cm³].
.. note:: DAGMC can compute cell volumes exactly from the triangulated
mesh surfaces. Specifying a manual volume risks inconsistency
with that capability.
*Default*: None
The following standard ``<cell>`` attributes are **not** supported inside
``<dagmc_universe>`` and will raise an error if present: ``region``,
``fill``, ``universe``, ``translation``, ``rotation``.
.. deprecated::
The ``<material_overrides>`` sub-element (containing ``<cell_override>``
children with ``<material_ids>``) is deprecated. A deprecation warning is
emitted and the overrides are converted to the ``<cell>`` format at parse
time. It is an error to specify both ``<material_overrides>`` and
``<cell>`` sub-elements on the same ``<dagmc_universe>``.
*Default*: None

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

@ -98,6 +98,11 @@ sub-elements:
A list of strings representing the nuclide, to define specific
define specific target nuclide collisions to be banked.
.. note::
Electron and positron collision-track events are not associated with
a specific nuclide. If a ``nuclides`` entry is specified, these events
are omitted.
*Default*: None
:reactions:
@ -606,30 +611,30 @@ found in the :ref:`random ray user guide <random_ray>`.
*Default*: None
:adjoint_source:
Specifies an adjoint fixed source for adjoint transport simulations, and
follows the format for :ref:`source_element`. The distributions which make
up the adjoint source are subject to the same restrictions as forward
Specifies an adjoint fixed source for adjoint transport simulations, and
follows the format for :ref:`source_element`. The distributions which make
up the adjoint source are subject to the same restrictions as forward
fixed sources in Random Ray mode.
*Default*: None
:adjoint:
Specifies whether to perform adjoint transport. The default is 'False',
Specifies whether to perform adjoint transport. The default is 'False',
corresponding to forward transport.
*Default*: None
:volume_estimator:
Specifies choice of volume estimator for the random ray solver. Options
Specifies choice of volume estimator for the random ray solver. Options
are 'naive', 'simulation_averaged', or 'hybrid'. The default is 'hybrid'.
*Default*: None
:volume_normalized_flux_tallies:
Specifies whether to normalize flux tallies by volume (bool). The
default is 'False'. When enabled, flux tallies will be reported in units
of cm/cm^3. When disabled, flux tallies will be reported in units of cm
(i.e., total distance traveled by neutrons in the spatial tally
Specifies whether to normalize flux tallies by volume (bool). The
default is 'False'. When enabled, flux tallies will be reported in units
of cm/cm^3. When disabled, flux tallies will be reported in units of cm
(i.e., total distance traveled by neutrons in the spatial tally
region).
*Default*: None
@ -741,14 +746,16 @@ pseudo-random number generator.
*Default*: 1
--------------------
``<stride>`` Element
--------------------
-----------------------------------
``<shared_secondary_bank>`` Element
-----------------------------------
The ``stride`` element is used to specify how many random numbers are allocated
for each source particle history.
*Default*: 152,917
The ``shared_secondary_bank`` element indicates whether to use a shared
secondary particle bank. When enabled, secondary particles are collected into
a global bank, sorted for reproducibility, and load-balanced across MPI ranks
between generations. If not specified, the shared secondary bank is enabled
automatically for fixed-source simulations with weight windows active, and
disabled otherwise.
.. _source_element:
@ -1058,17 +1065,19 @@ variable and whose sub-elements/attributes are as follows:
:type:
The type of the distribution. Valid options are "uniform", "discrete",
"tabular", "maxwell", "watt", and "mixture". The "uniform" option produces
variates sampled from a uniform distribution over a finite interval. The
"discrete" option produces random variates that can assume a finite number
of values (i.e., a distribution characterized by a probability mass function).
The "tabular" option produces random variates sampled from a tabulated
distribution where the density function is either a histogram or
"tabular", "maxwell", "watt", "mixture", and "decay_spectrum". The "uniform"
option produces variates sampled from a uniform distribution over a finite
interval. The "discrete" option produces random variates that can assume a
finite number of values (i.e., a distribution characterized by a probability
mass function). The "tabular" option produces random variates sampled from a
tabulated distribution where the density function is either a histogram or
linearly-interpolated between tabulated points. The "watt" option produces
random variates is sampled from a Watt fission spectrum (only used for
energies). The "maxwell" option produce variates sampled from a Maxwell
fission spectrum (only used for energies). The "mixture" option produces samples
from univariate sub-distributions with given probabilities.
fission spectrum (only used for energies). The "mixture" option produces
samples from univariate sub-distributions with given probabilities. The
"decay_spectrum" option produces photon energies sampled from decay photon
spectra in a depletion chain (only used for energies).
*Default*: None
@ -1086,6 +1095,10 @@ variable and whose sub-elements/attributes are as follows:
:math:`(x,p)` pairs defining the discrete/tabular distribution. All :math:`x`
points are given first followed by corresponding :math:`p` points.
For a "decay_spectrum" distribution, ``parameters`` gives the atom densities
in [atom/b-cm] for the nuclides listed in the ``nuclides`` element, in the
same order.
For a "watt" distribution, ``parameters`` should be given as two real numbers
:math:`a` and :math:`b` that parameterize the distribution :math:`p(x) dx = c
e^{-x/a} \sinh \sqrt{b \, x} dx`.
@ -1115,6 +1128,21 @@ variable and whose sub-elements/attributes are as follows:
This sub-element of a ``pair`` element provides information on the
corresponding univariate distribution.
:volume:
For a "decay_spectrum" distribution, this attribute specifies the source
region volume in cm\ :sup:`3`. It is used together with atom densities to
determine the absolute photon emission rate. When a source uses a
"decay_spectrum" energy distribution, the source strength is set from this
emission rate.
:nuclides:
For a "decay_spectrum" distribution, this element specifies a
whitespace-separated list of nuclide names contributing to the decay photon
source. The atom densities for these nuclides are given by the ``parameters``
element in the same order. Nuclides are resolved against the depletion chain,
and nuclides without decay photon spectra do not contribute to the
distribution.
:bias:
This optional element specifies a biased distribution for importance sampling.
For continuous distributions, the ``bias`` element should contain another
@ -1135,23 +1163,6 @@ based on constraints.
*Default*: 0.05
-------------------------
``<state_point>`` Element
-------------------------
The ``<state_point>`` element indicates at what batches a state point file
should be written. A state point file can be used to restart a run or to get
tally results at any batch. The default behavior when using this tag is to
write out the source bank in the state_point file. This behavior can be
customized by using the ``<source_point>`` element. This element has the
following attributes/sub-elements:
:batches:
A list of integers separated by spaces indicating at what batches a state
point file should be written.
*Default*: Last batch only
--------------------------
``<source_point>`` Element
--------------------------
@ -1201,6 +1212,32 @@ attributes/sub-elements:
*Default*: false
-------------------------
``<state_point>`` Element
-------------------------
The ``<state_point>`` element indicates at what batches a state point file
should be written. A state point file can be used to restart a run or to get
tally results at any batch. The default behavior when using this tag is to
write out the source bank in the state_point file. This behavior can be
customized by using the ``<source_point>`` element. This element has the
following attributes/sub-elements:
:batches:
A list of integers separated by spaces indicating at what batches a state
point file should be written.
*Default*: Last batch only
--------------------
``<stride>`` Element
--------------------
The ``stride`` element is used to specify how many random numbers are allocated
for each source particle history.
*Default*: 152,917
------------------------------
``<surf_source_read>`` Element
------------------------------
@ -1300,7 +1337,7 @@ The ``<surface_grazing_cutoff>`` element specifies the surface flux cosine cutof
``<surface_grazing_ratio>`` Element
-----------------------------------
The ``<surface_grazing_ratio>`` element specifies the surface flux cosine
The ``<surface_grazing_ratio>`` element specifies the surface flux cosine
substitution ratio.
*Default*: 0.5
@ -1725,11 +1762,11 @@ mesh-based weight windows.
The ratio of the lower to upper weight window bounds.
*Default*: 5.0
For FW-CADIS:
:targets:
A sequence of IDs corresponding to the tallies which cover phase
A sequence of IDs corresponding to the tallies which cover phase
space regions of interest for local variance reduction.
*Default*: None

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

@ -22,6 +22,7 @@ Univariate Probability Distributions
openmc.stats.Legendre
openmc.stats.Mixture
openmc.stats.Normal
openmc.stats.DecaySpectrum
.. autosummary::
:toctree: generated

39
docs/source/usersguide/depletion.rst
View file

@ -449,3 +449,42 @@ to transfer xenon from one material to another, you'd use::
...
integrator.add_transfer_rate(mat1, ['Xe'], 0.1, destination_material=mat2)
Comparing to Other Codes
========================
Comparing depletion results from OpenMC with those from another code, such as
MCNP or Serpent, requires more than constructing equivalent transport models.
At each depletion step, differences in the transport solution, nuclear data,
reaction rate normalization, and numerical integration can all affect the
result. Small differences can also accumulate over successive depletion steps.
For a meaningful comparison, align as many of the following inputs and methods
as possible:
- Geometry and material definitions and associated physical properties such as
temperature
- Neutron cross section library (e.g., ENDF/B-VIII.0)
- Treatment of thermal scattering and unresolved resonance probability tables
- Neutron reactions accounted for in the depletion chain
- Decay data in the depletion chain
- Isomeric branching ratios for reactions in the depletion chain
- Fission product yields in the depletion chain
- Fission product yield interpolation method
(``CoupledOperator(fission_yield_mode=...)``)
- Reaction rate normalization, including fission Q values
(``CoupledOperator(fission_q=...)``)
- Depletion integration method (``PredictorIntegrator``, ``CECMIntegrator``,
etc.) and time-step sizes
When comparing to codes that use ACE format cross sections, it is recommended to
directly convert the ACE files to HDF5 format using functionality from the
:mod:`openmc.data` module (see :ref:`create_xs_library`). Some of the
LANL-distributed ACE libraries used with MCNP have also been converted to HDF5
format and are available for download at https://openmc.org/data.
Even after these choices have been aligned, exact agreement should not be
expected. Codes may use different approximations or numerical methods that
cannot be configured identically. When investigating a discrepancy, first
compare transport results and one-group reaction rates at the initial time, then
compare changes over subsequent timesteps.

10
docs/source/usersguide/processing.rst
View file

@ -97,7 +97,15 @@ VTK Mesh File Generation
------------------------
VTK files of OpenMC meshes can be created using the
:meth:`openmc.Mesh.write_data_to_vtk` method. Data can be applied to the
:meth:`openmc.Mesh.write_data_to_vtk` method. This method supports several VTK
formats depending on the mesh type. Structured meshes
(:class:`~openmc.RegularMesh`, :class:`~openmc.RectilinearMesh`,
:class:`~openmc.CylindricalMesh`, and :class:`~openmc.SphericalMesh`) can be
exported to legacy VTK format (``.vtk``). The :class:`~openmc.UnstructuredMesh`
class supports VTK unstructured grid formats (``.vtu``) as well as an HDF5-based
format (``.vtkhdf``) that does not require the ``vtk`` module to write.
Data can be applied to the
elements of the resulting mesh from mesh filter objects. This data can be
provided either as a flat array or, in the case of structured meshes
(:class:`~openmc.RegularMesh`, :class:`~openmc.RectilinearMesh`,

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

@ -792,6 +792,11 @@ collision_track.h5 file at the end of the simulation. The file contains
300 recorded collisions that occurred in materials with IDs 1 or 2, involving
fission or (n,2n) reactions on the nuclides U-238 or O-16, within cells
with IDs 5 and 12.
.. note::
Electron and positron collision-track events are not associated with a
specific nuclide. If a ``nuclides`` entry is specified, these events are omitted.
The file can be read using :func:`openmc.read_collision_track_file`.
The example below shows how to extract the data from the collision_track
feature and displays the fields stored in the file:

8
include/openmc/bank.h
View file

@ -34,18 +34,24 @@ extern vector<vector<double>> ifp_fission_lifetime_bank;
extern vector<int64_t> progeny_per_particle;
extern SharedArray<SourceSite> shared_secondary_bank_read;
extern SharedArray<SourceSite> shared_secondary_bank_write;
} // namespace simulation
//==============================================================================
// Non-member functions
//==============================================================================
void sort_fission_bank();
void sort_bank(SharedArray<SourceSite>& bank, bool is_fission_bank);
void free_memory_bank();
void init_fission_bank(int64_t max);
int64_t synchronize_global_secondary_bank(
SharedArray<SourceSite>& shared_secondary_bank);
} // namespace openmc
#endif // OPENMC_BANK_H

5
include/openmc/capi.h
View file

@ -123,8 +123,13 @@ int openmc_new_filter(const char* type, int32_t* index);
int openmc_next_batch(int* status);
int openmc_nuclide_name(int index, const char** name);
int openmc_plot_geometry();
// Deprecated; use openmc_slice_data.
int openmc_id_map(const void* slice, int32_t* data_out);
// Deprecated; use openmc_slice_data.
int openmc_property_map(const void* slice, double* data_out);
int openmc_slice_data(const double origin[3], const double u_span[3],
const double v_span[3], const size_t pixels[2], bool show_overlaps, int level,
int32_t filter_index, int32_t* geom_data, double* property_data);
int openmc_get_plot_index(int32_t id, int32_t* index);
int openmc_plot_get_id(int32_t index, int32_t* id);
int openmc_plot_set_id(int32_t index, int32_t id);

24
include/openmc/cell.h
View file

@ -145,6 +145,30 @@ private:
bool simple_; //!< Does the region contain only intersections?
};
//==============================================================================
// XML parsing helpers for <cell> nodes
//==============================================================================
//! Parse material IDs from a <cell> XML node.
//! \param node XML node containing a "material" attribute or child element
//! \param cell_id Cell ID used in error messages
//! \return Vector of material IDs (MATERIAL_VOID for "void")
vector<int32_t> parse_cell_material_xml(pugi::xml_node node, int32_t cell_id);
//! Parse temperatures in [K] from a <cell> XML node.
//! Validates that all values are non-negative and the list is non-empty.
//! \param node XML node containing a "temperature" attribute or child element
//! \param cell_id Cell ID used in error messages
//! \return Vector of temperatures in [K]
vector<double> parse_cell_temperature_xml(pugi::xml_node node, int32_t cell_id);
//! Parse densities in [g/cm³] from a <cell> XML node.
//! Validates that all values are positive and the list is non-empty.
//! \param node XML node containing a "density" attribute or child element
//! \param cell_id Cell ID used in error messages
//! \return Vector of densities in [g/cm³]
vector<double> parse_cell_density_xml(pugi::xml_node node, int32_t cell_id);
//==============================================================================
class Cell {

2
include/openmc/chain.h
View file

@ -101,6 +101,8 @@ extern vector<unique_ptr<ChainNuclide>> chain_nuclides;
void read_chain_file_xml();
void free_memory_chain();
} // namespace openmc
#endif // OPENMC_CHAIN_H

2
include/openmc/constants.h
View file

@ -35,7 +35,7 @@ constexpr array<int, 2> VERSION_VOXEL {2, 0};
constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
constexpr array<int, 2> VERSION_PROPERTIES {1, 1};
constexpr array<int, 2> VERSION_WEIGHT_WINDOWS {1, 0};
constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 1};
constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 2};
// ============================================================================
// ADJUSTABLE PARAMETERS

19
include/openmc/dagmc.h
View file

@ -94,6 +94,10 @@ private:
class DAGUniverse : public Universe {
public:
using MaterialOverrides = std::unordered_map<int32_t, vector<int32_t>>;
using TemperatureOverrides = std::unordered_map<int32_t, vector<double>>;
using DensityOverrides = std::unordered_map<int32_t, vector<double>>;
explicit DAGUniverse(pugi::xml_node node);
//! Create a new DAGMC universe
@ -112,6 +116,9 @@ public:
//! Initialize the DAGMC accel. data structures, indices, material
//! assignments, etc.
void initialize();
void initialize(const MaterialOverrides& material_overrides,
const TemperatureOverrides& temperature_overrides,
const DensityOverrides& density_overrides = {});
//! Reads UWUW materials and returns an ID map
void read_uwuw_materials();
@ -146,7 +153,8 @@ public:
//! Assign a material overriding normal assignement to a cell
//! \param[in] c The OpenMC cell to which the material is assigned
void override_assign_material(std::unique_ptr<DAGCell>& c) const;
void override_assign_material(std::unique_ptr<DAGCell>& c,
const MaterialOverrides& material_overrides) const;
//! Return the index into the model cells vector for a given DAGMC volume
//! handle in the universe
@ -187,7 +195,9 @@ private:
void set_id(); //!< Deduce the universe id from model::universes
void init_dagmc(); //!< Create and initialise DAGMC pointer
void init_metadata(); //!< Create and initialise dagmcMetaData pointer
void init_geometry(); //!< Create cells and surfaces from DAGMC entities
void init_geometry(const MaterialOverrides& material_overrides,
const TemperatureOverrides& temperature_overrides,
const DensityOverrides& density_overrides);
std::string
filename_; //!< Name of the DAGMC file used to create this universe
@ -201,11 +211,6 @@ private:
//!< generate new material IDs for the universe
bool has_graveyard_; //!< Indicates if the DAGMC geometry has a "graveyard"
//!< volume
std::unordered_map<int32_t, vector<int32_t>>
material_overrides_; //!< Map of material overrides
//!< keys correspond to the DAGMCCell id
//!< values are a list of material ids used
//!< for the override
};
//==============================================================================

65
include/openmc/distribution.h
View file

@ -407,6 +407,71 @@ private:
double integral_; //!< Integral of distribution
};
//==============================================================================
// DecaySpectrum — non-owning mixture of decay photon distributions
//==============================================================================
//! Energy distribution formed by mixing multiple decay photon spectra.
//!
//! Unlike the general Mixture distribution, this class holds non-owning
//! pointers to the component distributions (which live in
//! data::chain_nuclides). Each component is weighted by the activity
//! (atoms * decay_constant) of the corresponding nuclide.
class DecaySpectrum : public Distribution {
public:
//============================================================================
// Types, aliases
struct Sample {
double energy;
double weight;
int parent_nuclide;
};
//============================================================================
// Constructors
//! Construct from an XML node containing nuclide names and atom densities.
//!
//! Reads child ``<nuclide>`` elements with ``name`` and ``density``
//! attributes, resolves them against the loaded depletion chain, and
//! constructs the mixed distribution.
explicit DecaySpectrum(pugi::xml_node node);
//============================================================================
// Methods
//! Sample a value from the distribution and return the parent nuclide index
//! \param seed Pseudorandom number seed pointer
//! \return (Sampled energy, sample weight, chain nuclide index)
Sample sample_with_parent(uint64_t* seed) const;
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return (sampled value, sample weight)
std::pair<double, double> sample(uint64_t* seed) const override;
double integral() const override;
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
//! Initialize decay spectrum sampling data
//! \param nuclide_indices Indices of decay photon emitters in
//! data::chain_nuclides
//! \param atoms Number of atoms for each component.
void init(vector<int> nuclide_indices, const vector<double>& atoms);
vector<int> nuclide_indices_; //!< Indices of emitting nuclides in the chain
DiscreteIndex di_; //!< Discrete index for component selection
double integral_; //!< Total photon emission rate
};
} // namespace openmc
#endif // OPENMC_DISTRIBUTION_H

13
include/openmc/event.h
View file

@ -112,6 +112,19 @@ void process_collision_events();
//! \param n_particles The number of particles in the particle buffer
void process_death_events(int64_t n_particles);
//! Process event queues until all are empty. Each iteration processes the
//! longest queue first to maximize vectorization efficiency.
void process_transport_events();
//! Initialize secondary particles from a shared secondary bank for
//! event-based transport
//
//! \param n_particles The number of particles to initialize
//! \param offset The offset index in the shared secondary bank
//! \param shared_secondary_bank The shared secondary bank to read from
void process_init_secondary_events(int64_t n_particles, int64_t offset,
const SharedArray<SourceSite>& shared_secondary_bank);
} // namespace openmc
#endif // OPENMC_EVENT_H

2
include/openmc/nuclide.h
View file

@ -88,6 +88,8 @@ public:
void rebuild_derived_xs();
std::vector<double> get_xs(int MT, int T_index) const;
std::vector<double> get_energy_grid(int T_index) const;
//! Return a ParticleType object representing this nuclide
ParticleType particle_type() const { return {Z_, A_, metastable_}; }
//============================================================================
// Data members

3
include/openmc/particle.h
View file

@ -71,7 +71,8 @@ public:
void event_advance();
void event_cross_surface();
void event_collide();
void event_revive_from_secondary();
void event_revive_from_secondary(const SourceSite& site);
void event_check_limit_and_revive();
void event_death();
//! pulse-height recording

34
include/openmc/particle_data.h
View file

@ -50,8 +50,11 @@ struct SourceSite {
// Extra attributes that don't show up in source written to file
int parent_nuclide {-1};
int64_t parent_id;
int64_t progeny_id;
int64_t parent_id {0};
int64_t progeny_id {0};
double wgt_born {1.0};
double wgt_ww_born {-1.0};
int64_t n_split {0};
};
struct CollisionTrackSite {
@ -533,14 +536,14 @@ private:
uint64_t seeds_[N_STREAMS];
int stream_;
vector<SourceSite> secondary_bank_;
vector<SourceSite> local_secondary_bank_;
// Keep track of how many secondary particles were created in the collision
// and what the starting index is in the secondary bank for this particle
int n_secondaries_ {0};
int secondary_bank_index_ {0};
int64_t current_work_;
int64_t current_work_ {0};
vector<double> flux_derivs_;
@ -563,7 +566,9 @@ private:
int n_event_ {0};
int n_split_ {0};
int64_t n_tracks_ {0}; //!< number of tracks in this particle history
int64_t n_split_ {0};
double ww_factor_ {0.0};
int64_t n_progeny_ {0};
@ -693,12 +698,14 @@ public:
int& stream() { return stream_; }
// secondary particle bank
SourceSite& secondary_bank(int i) { return secondary_bank_[i]; }
const SourceSite& secondary_bank(int i) const { return secondary_bank_[i]; }
decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; }
decltype(secondary_bank_) const& secondary_bank() const
SourceSite& local_secondary_bank(int i) { return local_secondary_bank_[i]; }
const SourceSite& local_secondary_bank(int i) const
{
return secondary_bank_;
return local_secondary_bank_[i];
}
decltype(local_secondary_bank_)& local_secondary_bank()
{
return local_secondary_bank_;
}
// Number of secondaries created in a collision
@ -747,13 +754,16 @@ public:
int& n_event() { return n_event_; }
// Number of times variance reduction has caused a particle split
int n_split() const { return n_split_; }
int& n_split() { return n_split_; }
int64_t n_split() const { return n_split_; }
int64_t& n_split() { return n_split_; }
// Particle-specific factor for on-the-fly weight window adjustment
double ww_factor() const { return ww_factor_; }
double& ww_factor() { return ww_factor_; }
// Number of tracks in this particle history
int64_t& n_tracks() { return n_tracks_; }
// Number of progeny produced by this particle
int64_t& n_progeny() { return n_progeny_; }

114
include/openmc/plot.h
View file

@ -18,6 +18,8 @@
#include "openmc/position.h"
#include "openmc/random_lcg.h"
#include "openmc/ray.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/filter_match.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -148,10 +150,11 @@ public:
struct IdData {
// Constructor
IdData(size_t h_res, size_t v_res);
IdData(size_t h_res, size_t v_res, bool include_filter = false);
// Methods
void set_value(size_t y, size_t x, const GeometryState& p, int level);
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);
// Members
@ -160,16 +163,34 @@ struct IdData {
struct PropertyData {
// Constructor
PropertyData(size_t h_res, size_t v_res);
PropertyData(size_t h_res, size_t v_res, bool include_filter = false);
// Methods
void set_value(size_t y, size_t x, const GeometryState& p, int level);
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);
// Members
tensor::Tensor<double> data_; //!< 2D array of temperature & density data
};
struct RasterData {
// Constructor
RasterData(size_t h_res, size_t v_res, bool include_filter = false);
// 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);
// Members
tensor::Tensor<int32_t>
id_data_; //!< [v_res, h_res, 3 or 4]: cell, instance, mat, [filter_bin]
tensor::Tensor<double>
property_data_; //!< [v_res, h_res, 2]: temperature, density
bool include_filter_; //!< Whether filter bin index is included
};
//===============================================================================
// Plot class
//===============================================================================
@ -177,7 +198,7 @@ struct PropertyData {
class SlicePlotBase {
public:
template<class T>
T get_map() const;
T get_map(int32_t filter_index = -1) const;
enum class PlotBasis { xy = 1, xz = 2, yz = 3 };
@ -188,70 +209,65 @@ public:
// Members
public:
Position origin_; //!< Plot origin in geometry
Position width_; //!< Plot width in geometry
PlotBasis basis_; //!< Plot basis (XY/XZ/YZ)
array<size_t, 3> pixels_; //!< Plot size in pixels
bool slice_color_overlaps_; //!< Show overlapping cells?
int slice_level_ {-1}; //!< Plot universe level
Position origin_; //!< Plot origin in geometry
Direction u_span_; //!< Full-width span vector in geometry
Direction v_span_; //!< Full-height span vector in geometry
array<size_t, 3> pixels_; //!< Plot size in pixels
bool show_overlaps_; //!< Show overlapping cells?
int slice_level_ {-1}; //!< Plot universe level
private:
};
template<class T>
T SlicePlotBase::get_map() const
T SlicePlotBase::get_map(int32_t filter_index) const
{
size_t width = pixels_[0];
size_t height = pixels_[1];
// get pixel size
double in_pixel = (width_[0]) / static_cast<double>(width);
double out_pixel = (width_[1]) / static_cast<double>(height);
// size data array
T data(width, height);
// setup basis indices and initial position centered on pixel
int in_i, out_i;
Position xyz = origin_;
switch (basis_) {
case PlotBasis::xy:
in_i = 0;
out_i = 1;
break;
case PlotBasis::xz:
in_i = 0;
out_i = 2;
break;
case PlotBasis::yz:
in_i = 1;
out_i = 2;
break;
default:
UNREACHABLE();
// Determine if filter is being used
bool include_filter = (filter_index >= 0);
Filter* filter = nullptr;
if (include_filter) {
filter = model::tally_filters[filter_index].get();
}
// set initial position
xyz[in_i] = origin_[in_i] - width_[0] / 2. + in_pixel / 2.;
xyz[out_i] = origin_[out_i] + width_[1] / 2. - out_pixel / 2.;
// size data array
T data(width, height, include_filter);
// arbitrary direction
Direction dir = {1. / std::sqrt(2.), 1. / std::sqrt(2.), 0.0};
// compute pixel steps and top-left pixel center
Direction u_step = u_span_ / static_cast<double>(width);
Direction v_step = v_span_ / static_cast<double>(height);
Position start =
origin_ - 0.5 * u_span_ + 0.5 * v_span_ + 0.5 * u_step - 0.5 * v_step;
// Validate that span vectors define a valid plane
Position cross = u_span_.cross(v_span_);
if (cross.norm() == 0.0) {
fatal_error("Slice span vectors are invalid (zero area).");
}
// Use an arbitrary direction that is not aligned with any coordinate axis.
// The direction has no physical meaning for plotting but is used by
// Surface::sense() to break ties when a pixel is coincident with a surface.
Direction dir = {1.0 / std::sqrt(2.0), 1.0 / std::sqrt(2.0), 0.0};
#pragma omp parallel
{
GeometryState p;
p.r() = xyz;
Particle p;
p.r() = start;
p.u() = dir;
p.coord(0).universe() = model::root_universe;
int level = slice_level_;
int j {};
FilterMatch match;
#pragma omp for
for (int y = 0; y < height; y++) {
p.r()[out_i] = xyz[out_i] - out_pixel * y;
Position row = start - v_step * static_cast<double>(y);
for (int x = 0; x < width; x++) {
p.r()[in_i] = xyz[in_i] + in_pixel * x;
p.r() = row + u_step * static_cast<double>(x);
p.n_coord() = 1;
// local variables
bool found_cell = exhaustive_find_cell(p);
@ -260,9 +276,9 @@ T SlicePlotBase::get_map() const
j = level;
}
if (found_cell) {
data.set_value(y, x, p, j);
data.set_value(y, x, p, j, filter, &match);
}
if (slice_color_overlaps_ && check_cell_overlap(p, false)) {
if (show_overlaps_ && check_cell_overlap(p, false)) {
data.set_overlap(y, x);
}
} // inner for
@ -297,6 +313,8 @@ public:
void print_info() const override;
PlotType type_; //!< Plot type (Slice/Voxel)
Position width_; //!< Axis-aligned width from plot.xml
PlotBasis basis_; //!< Basis from plot.xml for slice plots
int meshlines_width_; //!< Width of lines added to the plot
int index_meshlines_mesh_ {-1}; //!< Index of the mesh to draw on the plot
RGBColor meshlines_color_; //!< Color of meshlines on the plot

3
include/openmc/settings.h
View file

@ -98,7 +98,8 @@ extern bool uniform_source_sampling; //!< sample sources uniformly?
extern bool ufs_on; //!< uniform fission site method on?
extern bool urr_ptables_on; //!< use unresolved resonance prob. tables?
extern bool use_decay_photons; //!< use decay photons for D1S
extern "C" bool weight_windows_on; //!< are weight windows are enabled?
extern bool use_shared_secondary_bank; //!< Use shared bank for secondaries
extern "C" bool weight_windows_on; //!< are weight windows are enabled?
extern bool weight_window_checkpoint_surface; //!< enable weight window check
//!< upon surface crossing?
extern bool weight_window_checkpoint_collision; //!< enable weight window check

30
include/openmc/shared_array.h
View file

@ -4,6 +4,8 @@
//! \file shared_array.h
//! \brief Shared array data structure
#include <algorithm> // for copy_n
#include "openmc/memory.h"
namespace openmc {
@ -30,14 +32,12 @@ public:
//! Default constructor.
SharedArray() = default;
//! Construct a zero size container with space to hold capacity number of
//! elements.
//! Construct a container with `size` elements and capacity equal to `size`.
//
//! \param capacity The number of elements for the container to allocate
//! space for
SharedArray(int64_t capacity) : capacity_(capacity)
//! \param size The number of elements to allocate and initialize
SharedArray(int64_t size) : size_(size), capacity_(size)
{
data_ = make_unique<T[]>(capacity);
data_ = make_unique<T[]>(size);
}
//==========================================================================
@ -97,8 +97,26 @@ public:
capacity_ = 0;
}
//! Push back an element to the array, with capacity and reallocation behavior
//! as if this were a vector. This does not perform any thread safety checks.
//! If the size exceeds the capacity, then the capacity will double just as
//! with a vector. Data will be reallocated and moved to a new pointer and
//! copied in before the new item is appended. Old data will be freed.
void thread_unsafe_append(const T& value)
{
if (size_ == capacity_) {
int64_t new_capacity = capacity_ == 0 ? 8 : 2 * capacity_;
unique_ptr<T[]> new_data = make_unique<T[]>(new_capacity);
std::copy_n(data_.get(), size_, new_data.get());
data_ = std::move(new_data);
capacity_ = new_capacity;
}
data_[size_++] = value;
}
//! Return the number of elements in the container
int64_t size() { return size_; }
int64_t size() const { return size_; }
//! Resize the container to contain a specified number of elements. This is
//! useful in cases where the container is written to in a non-thread safe

33
include/openmc/simulation.h
View file

@ -49,6 +49,9 @@ extern const RegularMesh* ufs_mesh;
extern vector<double> k_generation;
extern vector<int64_t> work_index;
extern int64_t
simulation_tracks_completed; //!< Number of tracks completed on this rank
} // namespace simulation
//==============================================================================
@ -59,7 +62,7 @@ extern vector<int64_t> work_index;
void allocate_banks();
//! Determine number of particles to transport per process
void calculate_work();
void calculate_work(int64_t n_particles);
//! Initialize nuclear data before a simulation
void initialize_data();
@ -70,8 +73,9 @@ void initialize_batch();
//! Initialize a fission generation
void initialize_generation();
//! Full initialization of a particle history
void initialize_history(Particle& p, int64_t index_source);
//! Full initialization of a particle track
void initialize_particle_track(
Particle& p, int64_t index_source, bool is_secondary);
//! Finalize a batch
//!
@ -92,16 +96,35 @@ void broadcast_results();
void free_memory_simulation();
//! Simulate a single particle history (and all generated secondary particles,
//! if enabled), from birth to death
//! Compute unique particle ID from a 1-based source index
//! \param index_source 1-based source index within this rank's work
//! \return globally unique particle ID
int64_t compute_particle_id(int64_t index_source);
//! Compute the transport RNG seed from a particle ID
//! \param particle_id the particle's globally unique ID
//! \return seed value passed to init_particle_seeds()
int64_t compute_transport_seed(int64_t particle_id);
//! Simulate a single particle history from birth to death, inclusive of any
//! secondary particles. In shared secondary mode, only a single track is
//! transported and secondaries are deposited into a shared bank instead.
void transport_history_based_single_particle(Particle& p);
//! Simulate all particle histories using history-based parallelism
void transport_history_based();
//! Simulate all particles using history-based parallelism, with a shared
//! secondary bank
void transport_history_based_shared_secondary();
//! Simulate all particle histories using event-based parallelism
void transport_event_based();
//! Simulate all particles using event-based parallelism, with a shared
//! secondary bank
void transport_event_based_shared_secondary();
} // namespace openmc
#endif // OPENMC_SIMULATION_H

230
openmc/dagmc.py
View file

@ -36,12 +36,6 @@ class DAGMCUniverse(openmc.UniverseBase):
auto_mat_ids : bool
Set IDs automatically on initialization (True) or report overlaps in ID
space between OpenMC and UWUW materials (False)
material_overrides : dict, optional
A dictionary of material overrides. The keys are material name strings
and the values are Iterables of openmc.Material objects. If a material
name is found in the DAGMC file, the material will be replaced with the
openmc.Material object in the value.
Attributes
----------
id : int
@ -78,15 +72,6 @@ class DAGMCUniverse(openmc.UniverseBase):
The number of surfaces in the model.
.. versionadded:: 0.13.2
material_overrides : dict
A dictionary of material overrides. Keys are cell IDs; values are
iterables of :class:`openmc.Material` objects. The material assignment
of each DAGMC cell ID key will be replaced with the
:class:`~openmc.Material` object in the value. If the value contains
multiple :class:`~openmc.Material` objects, each Material in the list
will be assigned to the corresponding instance of the cell.
.. versionadded:: 0.15.1
"""
def __init__(self,
@ -94,16 +79,12 @@ class DAGMCUniverse(openmc.UniverseBase):
universe_id=None,
name='',
auto_geom_ids=False,
auto_mat_ids=False,
material_overrides=None):
auto_mat_ids=False):
super().__init__(universe_id, name)
# Initialize class attributes
self.filename = filename
self.auto_geom_ids = auto_geom_ids
self.auto_mat_ids = auto_mat_ids
self._material_overrides = {}
if material_overrides is not None:
self.material_overrides = material_overrides
def __repr__(self):
string = super().__repr__()
@ -130,47 +111,17 @@ class DAGMCUniverse(openmc.UniverseBase):
@property
def material_overrides(self):
return self._material_overrides
raise AttributeError(
"DAGMCUniverse.material_overrides has been removed. Use "
"DAGMCCell objects added via add_cell() to manage per-cell "
"material assignments.")
@material_overrides.setter
def material_overrides(self, val):
cv.check_type('material overrides', val, Mapping)
for key, value in val.items():
self.add_material_override(key, value)
def replace_material_assignment(self, material_name: str, material: openmc.Material):
"""Replace a material assignment within the DAGMC universe.
Replace the material assignment of all cells filled with a material in
the DAGMC universe. The universe must be synchronized in an initialized
Model (see :meth:`~openmc.DAGMCUniverse.sync_dagmc_cells`) before
calling this method.
.. versionadded:: 0.15.1
Parameters
----------
material_name : str
Material name to replace
material : openmc.Material
Material to replace the material_name with
"""
if material_name not in self.material_names:
raise ValueError(
f"No material with name '{material_name}' found in the DAGMC universe")
if not self.cells:
raise RuntimeError("This DAGMC universe has not been synchronized "
"on an initialized Model.")
for cell in self.cells.values():
if cell.fill is None:
continue
if isinstance(cell.fill, openmc.Iterable):
cell.fill = list(map(lambda x: material if x.name == material_name else x, cell.fill))
else:
cell.fill = material if cell.fill.name == material_name else cell.fill
raise AttributeError(
"DAGMCUniverse.material_overrides has been removed. Use "
"DAGMCCell objects added via add_cell() to manage per-cell "
"material assignments.")
def add_material_override(self, key, overrides=None):
"""Add a material override to the universe.
@ -201,7 +152,10 @@ class DAGMCUniverse(openmc.UniverseBase):
if key not in self.cells:
raise ValueError(f"Cell ID '{key}' not found in DAGMC universe")
self._material_overrides[key] = overrides
if len(overrides) == 1:
self.cells[key].fill = overrides[0]
else:
self.cells[key].fill = list(overrides)
@property
def auto_geom_ids(self):
@ -290,12 +244,6 @@ class DAGMCUniverse(openmc.UniverseBase):
memo.add(self)
# Ensure that the material overrides are up-to-date
for cell in self.cells.values():
if cell.fill is None:
continue
self.add_material_override(cell, cell.fill)
# Set xml element values
dagmc_element = ET.Element('dagmc_universe')
dagmc_element.set('id', str(self.id))
@ -307,17 +255,10 @@ class DAGMCUniverse(openmc.UniverseBase):
if self.auto_mat_ids:
dagmc_element.set('auto_mat_ids', 'true')
dagmc_element.set('filename', str(self.filename))
if self._material_overrides:
mat_element = ET.Element('material_overrides')
for key in self._material_overrides:
cell_overrides = ET.Element('cell_override')
cell_overrides.set("id", str(key))
material_element = ET.Element('material_ids')
material_element.text = ' '.join(
str(t.id) for t in self._material_overrides[key])
cell_overrides.append(material_element)
mat_element.append(cell_overrides)
dagmc_element.append(mat_element)
if self.cells:
for cell in self.cells.values():
cell_element = cell.create_xml_subelement(xml_element, memo)
dagmc_element.append(cell_element)
xml_element.append(dagmc_element)
def bounding_region(
@ -442,7 +383,7 @@ class DAGMCUniverse(openmc.UniverseBase):
return out
@classmethod
def from_xml_element(cls, elem, mats = None):
def from_xml_element(cls, elem, mats=None):
"""Generate DAGMC universe from XML element
Parameters
@ -471,21 +412,56 @@ class DAGMCUniverse(openmc.UniverseBase):
out.auto_geom_ids = bool(get_text(elem, "auto_geom_ids"))
out.auto_mat_ids = bool(get_text(elem, "auto_mat_ids"))
el_mat_override = elem.find('material_overrides')
if el_mat_override is not None:
if mats is None:
raise ValueError("Material overrides found in DAGMC universe "
"but no materials were provided to populate "
"the mapping.")
out._material_overrides = {}
for elem in el_mat_override.findall('cell_override'):
cell_id = int(get_text(elem, 'id'))
mat_ids = get_elem_list(elem, "material_ids", str) or []
mat_objs = [mats[mat_id] for mat_id in mat_ids]
out._material_overrides[cell_id] = mat_objs
has_overrides = elem.find('material_overrides') is not None
has_cells = elem.find('cell') is not None
if has_overrides and has_cells:
raise ValueError(
"DAGMCUniverse cannot specify both <material_overrides> and "
"<cell> sub-elements. Use <cell> elements only.")
if has_overrides:
warnings.warn(
"DAGMCUniverse <material_overrides> is deprecated and will be "
"removed in a future version. Use nested <cell> elements "
"instead.", DeprecationWarning, stacklevel=2)
out._parse_legacy_material_overrides(elem, mats)
elif has_cells:
out._parse_cell_overrides(elem, mats)
return out
def _parse_legacy_material_overrides(self, elem, mats):
"""Parse the deprecated <material_overrides> XML format and populate
the universe with equivalent DAGMCCell objects."""
if mats is None:
raise ValueError(
"DAGMC material overrides found but no materials were "
"provided to populate the mapping.")
mo_elem = elem.find('material_overrides')
for co_elem in mo_elem.findall('cell_override'):
cell_id = int(get_text(co_elem, 'id'))
mat_ids = co_elem.find('material_ids').text.split()
fill_objs = [mats[mid] for mid in mat_ids]
fill = fill_objs[0] if len(fill_objs) == 1 else fill_objs
if cell_id in self.cells:
raise ValueError(
f"Duplicate DAGMC cell override specified for cell {cell_id}.")
self.add_cell(DAGMCCell(cell_id=cell_id, fill=fill))
def _parse_cell_overrides(self, elem, mats):
if mats is None:
raise ValueError("DAGMC cell overrides found in DAGMC universe but "
"no materials were provided to populate the "
"mapping.")
for cell_elem in elem.findall('cell'):
cell_id = int(get_text(cell_elem, 'id'))
if cell_id in self.cells:
raise ValueError(
f"Duplicate DAGMC cell override specified for cell {cell_id}.")
DAGMCCell.from_xml_element(cell_elem, mats, self)
def _partial_deepcopy(self):
"""Clone all of the openmc.DAGMCUniverse object's attributes except for
its cells, as they are copied within the clone function. This should
@ -565,7 +541,13 @@ class DAGMCUniverse(openmc.UniverseBase):
fill = [mats_per_id[mat.id] for mat in dag_cell.fill if mat]
else:
fill = mats_per_id[dag_cell.fill.id] if dag_cell.fill else None
self.add_cell(openmc.DAGMCCell(cell_id=dag_cell_id, fill=fill))
name = dag_cell.name
if dag_cell_id in self._cells:
self._cells[dag_cell_id].name = name
self._cells[dag_cell_id].fill = fill
else:
self.add_cell(
openmc.DAGMCCell(cell_id=dag_cell_id, name=name, fill=fill))
@add_plot_params
def plot(self, *args, **kwargs):
@ -594,6 +576,14 @@ class DAGMCCell(openmc.Cell):
DAG_parent_universe : int
The parent universe of the cell.
Notes
-----
DAGMC geometries are composed of triangulated surfaces, which means cell
volumes can in principle be computed exactly (e.g. via mesh-based
integration). Manually specifying :attr:`volume` overrides any such
calculation and may introduce inconsistencies if the value does not
accurately reflect the true geometric volume.
"""
def __init__(self, cell_id=None, name='', fill=None):
super().__init__(cell_id, name, fill, None)
@ -625,8 +615,62 @@ class DAGMCCell(openmc.Cell):
raise TypeError("plot is not available for DAGMC cells.")
def create_xml_subelement(self, xml_element, memo=None):
raise TypeError("create_xml_subelement is not available for DAGMC cells.")
if self.fill_type not in ('void', 'material', 'distribmat'):
raise TypeError("DAGMC cell overrides currently only support "
"material fills.")
if self.temperature is not None and self.fill_type not in (
'material', 'distribmat'
):
raise TypeError("DAGMC cell temperature overrides require a "
"material fill.")
if self.density is not None and self.fill_type not in ('material', 'distribmat'):
raise TypeError("DAGMC cell density overrides require a "
"material fill.")
if any(getattr(self, attr) is not None for attr in ('translation', 'rotation')):
raise TypeError("DAGMC cell overrides do not support translation "
"or rotation.")
return super().create_xml_subelement(xml_element, memo)
@classmethod
def from_xml_element(cls, elem, surfaces, materials, get_universe):
raise TypeError("from_xml_element is not available for DAGMC cells.")
def from_xml_element(cls, elem, mats, universe):
"""Generate a DAGMCCell from an XML <cell> override element.
Parameters
----------
elem : lxml.etree._Element
`<cell>` element containing a DAGMC cell property override
mats : dict
Dictionary mapping material ID strings to :class:`openmc.Material`
instances
universe : DAGMCUniverse
Universe to add the parsed cell to.
Returns
-------
DAGMCCell
DAGMCCell instance
"""
if not isinstance(universe, DAGMCUniverse):
raise TypeError(
f"universe must be a DAGMCUniverse instance, "
f"got {type(universe).__name__}.")
cell_id = int(get_text(elem, 'id'))
# Validate attributes that are unsupported for DAGMC cell overrides
for tag in ('region', 'fill', 'universe'):
if get_text(elem, tag) is not None:
raise ValueError(
f"DAGMC cell {cell_id} override cannot specify '{tag}'.")
for tag in ('translation', 'rotation'):
if get_text(elem, tag) is not None:
raise ValueError(
f"DAGMC cell {cell_id} override does not support "
f"'{tag}'.")
if get_elem_list(elem, 'material', str) is None:
raise ValueError(
f"DAGMC cell {cell_id} must specify a material override.")
return super().from_xml_element(
elem, surfaces={}, materials=mats,
get_universe=lambda _: universe)

7
openmc/data/angle_distribution.py
View file

@ -10,8 +10,8 @@ from openmc.mixin import EqualityMixin
from openmc.stats import Univariate, Tabular, Uniform, Legendre
from .function import INTERPOLATION_SCHEME
from .data import EV_PER_MEV
from .endf import get_head_record, get_cont_record, get_tab1_record, \
get_list_record, get_tab2_record
from .endf import as_evaluation, get_head_record, get_cont_record, \
get_tab1_record, get_list_record, get_tab2_record
class AngleDistribution(EqualityMixin):
@ -213,7 +213,7 @@ class AngleDistribution(EqualityMixin):
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF evaluation
mt : int
The MT value of the reaction to get angular distributions for
@ -224,6 +224,7 @@ class AngleDistribution(EqualityMixin):
Angular distribution
"""
ev = as_evaluation(ev)
file_obj = StringIO(ev.section[4, mt])
# Read HEAD record

25
openmc/data/decay.py
View file

@ -14,7 +14,8 @@ from openmc.mixin import EqualityMixin
from openmc.stats import Discrete, Tabular, Univariate, combine_distributions
from .data import gnds_name, zam
from .function import INTERPOLATION_SCHEME
from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
from .endf import (
as_evaluation, get_head_record, get_list_record, get_tab1_record)
# Gives name and (change in A, change in Z) resulting from decay
@ -75,7 +76,7 @@ class FissionProductYields(EqualityMixin):
Parameters
----------
ev_or_filename : str of openmc.data.endf.Evaluation
ev_or_filename : str, openmc.data.endf.Evaluation, or endf.Material
ENDF fission product yield evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -133,11 +134,7 @@ class FissionProductYields(EqualityMixin):
return energies, data
# Get evaluation if str is passed
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
ev = as_evaluation(ev_or_filename)
# Assign basic nuclide properties
self.nuclide = {
@ -164,7 +161,7 @@ class FissionProductYields(EqualityMixin):
Parameters
----------
ev_or_filename : str or openmc.data.endf.Evaluation
ev_or_filename : str, openmc.data.endf.Evaluation, or endf.Material
ENDF fission product yield evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -292,7 +289,7 @@ class Decay(EqualityMixin):
Parameters
----------
ev_or_filename : str of openmc.data.endf.Evaluation
ev_or_filename : str, openmc.data.endf.Evaluation, or endf.Material
ENDF radioactive decay data evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -323,11 +320,7 @@ class Decay(EqualityMixin):
"""
def __init__(self, ev_or_filename):
# Get evaluation if str is passed
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
ev = as_evaluation(ev_or_filename)
file_obj = StringIO(ev.section[8, 457])
@ -486,7 +479,7 @@ class Decay(EqualityMixin):
Parameters
----------
ev_or_filename : str or openmc.data.endf.Evaluation
ev_or_filename : str, openmc.data.endf.Evaluation, or endf.Material
ENDF radioactive decay data evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -649,5 +642,3 @@ def decay_energy(nuclide: str):
warn(f"Chain file '{chain_file}' does not have any decay energy.")
return _DECAY_ENERGY.get(nuclide, 0.0)

34
openmc/data/endf.py
View file

@ -12,7 +12,12 @@ import re
from .data import gnds_name
from .function import Tabulated1D
from endf.material import _LIBRARY, _SUBLIBRARY, get_materials as get_evaluations
from endf.material import (
Material,
_LIBRARY,
_SUBLIBRARY,
get_materials as get_evaluations,
)
from endf.incident_neutron import SUM_RULES
from endf.records import (
float_endf,
@ -44,7 +49,7 @@ class Evaluation:
Parameters
----------
filename_or_obj : str or file-like
filename_or_obj : str, file-like, or endf.Material
Path to ENDF file to read or an open file positioned at the start of an
ENDF material
@ -64,17 +69,25 @@ class Evaluation:
"""
def __init__(self, filename_or_obj):
self.section = {}
self.info = {}
self.target = {}
self.projectile = {}
self.reaction_list = []
if isinstance(filename_or_obj, Material):
self.section = dict(filename_or_obj.section_text)
self.section_data = filename_or_obj.section_data
self.material = filename_or_obj.MAT
self._read_header()
return
if isinstance(filename_or_obj, (str, PurePath)):
fh = open(str(filename_or_obj), 'r')
need_to_close = True
else:
fh = filename_or_obj
need_to_close = False
self.section = {}
self.info = {}
self.target = {}
self.projectile = {}
self.reaction_list = []
# Skip TPID record. Evaluators sometimes put in TPID records that are
# ill-formated because they lack MF/MT values or put them in the wrong
@ -199,3 +212,10 @@ class Evaluation:
self.target['mass_number'],
self.target['isomeric_state'])
def as_evaluation(ev_or_filename):
"""Return an object supporting OpenMC's legacy Evaluation interface."""
if isinstance(ev_or_filename, Evaluation):
return ev_or_filename
else:
return Evaluation(ev_or_filename)

7
openmc/data/fission_energy.py
View file

@ -5,7 +5,8 @@ from io import StringIO
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from .data import EV_PER_MEV
from .endf import get_cont_record, get_list_record, get_tab1_record, Evaluation
from .endf import (
as_evaluation, get_cont_record, get_list_record, get_tab1_record)
from .function import Function1D, Tabulated1D, Polynomial, sum_functions
@ -195,7 +196,7 @@ class FissionEnergyRelease(EqualityMixin):
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF evaluation
incident_neutron : openmc.data.IncidentNeutron
Corresponding incident neutron dataset
@ -206,7 +207,7 @@ class FissionEnergyRelease(EqualityMixin):
Fission energy release data
"""
cv.check_type('evaluation', ev, Evaluation)
ev = as_evaluation(ev)
# Check to make sure this ENDF file matches the expected isomer.
if ev.target['atomic_number'] != incident_neutron.atomic_number:

10
openmc/data/neutron.py
View file

@ -13,7 +13,8 @@ from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from .ace import Library, Table, get_table, get_metadata
from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV, gnds_name
from .endf import (
Evaluation, SUM_RULES, get_head_record, get_tab1_record, get_evaluations)
Evaluation, SUM_RULES, as_evaluation, get_head_record, get_tab1_record,
get_evaluations)
from .fission_energy import FissionEnergyRelease
from .function import Tabulated1D, Sum, ResonancesWithBackground
from .njoy import make_ace, make_pendf
@ -652,7 +653,7 @@ class IncidentNeutron(EqualityMixin):
Parameters
----------
ev_or_filename : openmc.data.endf.Evaluation or str
ev_or_filename : openmc.data.endf.Evaluation, endf.Material, or str
ENDF evaluation to read from. If given as a string, it is assumed to
be the filename for the ENDF file.
@ -666,10 +667,7 @@ class IncidentNeutron(EqualityMixin):
Incident neutron continuous-energy data
"""
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
ev = as_evaluation(ev_or_filename)
atomic_number = ev.target['atomic_number']
mass_number = ev.target['mass_number']

22
openmc/data/photon.py
View file

@ -15,7 +15,8 @@ from openmc.mixin import EqualityMixin
from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from .ace import Table, get_metadata, get_table
from .data import ATOMIC_SYMBOL, EV_PER_MEV
from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record
from .endf import (
as_evaluation, get_head_record, get_tab1_record, get_list_record)
from .function import Tabulated1D
@ -272,7 +273,7 @@ class AtomicRelaxation(EqualityMixin):
Parameters
----------
ev_or_filename : str or openmc.data.endf.Evaluation
ev_or_filename : str, openmc.data.endf.Evaluation, or endf.Material
ENDF atomic relaxation evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -282,10 +283,7 @@ class AtomicRelaxation(EqualityMixin):
Atomic relaxation data
"""
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
ev = as_evaluation(ev_or_filename)
# Atomic relaxation data is always MF=28, MT=533
if (28, 533) not in ev.section:
@ -606,10 +604,10 @@ class IncidentPhoton(EqualityMixin):
Parameters
----------
photoatomic : str or openmc.data.endf.Evaluation
photoatomic : str, openmc.data.endf.Evaluation, or endf.Material
ENDF photoatomic data evaluation to read from. If given as a string,
it is assumed to be the filename for the ENDF file.
relaxation : str or openmc.data.endf.Evaluation, optional
relaxation : str, openmc.data.endf.Evaluation, or endf.Material, optional
ENDF atomic relaxation data evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -619,10 +617,7 @@ class IncidentPhoton(EqualityMixin):
Photon interaction data
"""
if isinstance(photoatomic, Evaluation):
ev = photoatomic
else:
ev = Evaluation(photoatomic)
ev = as_evaluation(photoatomic)
Z = ev.target['atomic_number']
data = cls(Z)
@ -1071,7 +1066,7 @@ class PhotonReaction(EqualityMixin):
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF photo-atomic interaction data evaluation
mt : int
The MT value of the reaction to get data for
@ -1082,6 +1077,7 @@ class PhotonReaction(EqualityMixin):
Photon reaction data
"""
ev = as_evaluation(ev)
rx = cls(mt)
# Read photon cross section

7
openmc/data/reaction.py
View file

@ -13,8 +13,8 @@ from .angle_distribution import AngleDistribution
from .angle_energy import AngleEnergy
from .correlated import CorrelatedAngleEnergy
from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV
from .endf import get_head_record, get_tab1_record, get_list_record, \
get_tab2_record, get_cont_record
from .endf import as_evaluation, get_head_record, get_tab1_record, \
get_list_record, get_tab2_record, get_cont_record
from .energy_distribution import EnergyDistribution, LevelInelastic, \
DiscretePhoton
from .function import Tabulated1D, Polynomial
@ -1151,7 +1151,7 @@ class Reaction(EqualityMixin):
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF evaluation
mt : int
The MT value of the reaction to get data for
@ -1162,6 +1162,7 @@ class Reaction(EqualityMixin):
Reaction data
"""
ev = as_evaluation(ev)
rx = Reaction(mt)
# Integrated cross section

7
openmc/data/resonance.py
View file

@ -7,7 +7,9 @@ import pandas as pd
import openmc.checkvalue as cv
from .data import NEUTRON_MASS
from .endf import get_head_record, get_cont_record, get_tab1_record, get_list_record
from .endf import (
as_evaluation, get_head_record, get_cont_record, get_tab1_record,
get_list_record)
try:
from .reconstruct import wave_number, penetration_shift, reconstruct_mlbw, \
reconstruct_slbw, reconstruct_rm
@ -77,7 +79,7 @@ class Resonances:
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF evaluation
Returns
@ -86,6 +88,7 @@ class Resonances:
Resonance data
"""
ev = as_evaluation(ev)
file_obj = io.StringIO(ev.section[2, 151])
# Determine whether discrete or continuous representation

3
openmc/data/resonance_covariance.py
View file

@ -74,7 +74,7 @@ class ResonanceCovariances(Resonances):
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF evaluation
resonances : openmc.data.Resonance object
openmc.data.Resonanance object generated from the same evaluation
@ -86,6 +86,7 @@ class ResonanceCovariances(Resonances):
Resonance covariance data
"""
ev = endf.as_evaluation(ev)
file_obj = io.StringIO(ev.section[32, 151])
# Determine whether discrete or continuous representation

7
openmc/data/thermal.py
View file

@ -1042,7 +1042,7 @@ class ThermalScattering(EqualityMixin):
Parameters
----------
ev_or_filename : openmc.data.endf.Evaluation or str
ev_or_filename : openmc.data.endf.Evaluation, endf.Material, or str
ENDF evaluation to read from. If given as a string, it is assumed to
be the filename for the ENDF file.
divide_incoherent_elastic : bool
@ -1056,10 +1056,7 @@ class ThermalScattering(EqualityMixin):
Thermal scattering data
"""
if isinstance(ev_or_filename, endf.Evaluation):
ev = ev_or_filename
else:
ev = endf.Evaluation(ev_or_filename)
ev = endf.as_evaluation(ev_or_filename)
# Read incoherent inelastic data
assert (7, 4) in ev.section, 'No MF=7, MT=4 found in thermal scattering'

12
openmc/deplete/chain.py
View file

@ -319,16 +319,16 @@ class Chain:
String arguments in ``decay_files``, ``fpy_files``, and
``neutron_files`` will be treated as file names to be read.
Alternatively, :class:`openmc.data.endf.Evaluation` instances
can be included in these arguments.
Alternatively, :class:`openmc.data.endf.Evaluation` or
``endf.Material`` instances can be included in these arguments.
Parameters
----------
decay_files : list of str or openmc.data.endf.Evaluation
decay_files : list of str, openmc.data.endf.Evaluation, or endf.Material
List of ENDF decay sub-library files
fpy_files : list of str or openmc.data.endf.Evaluation
fpy_files : list of str, openmc.data.endf.Evaluation, or endf.Material
List of ENDF neutron-induced fission product yield sub-library files
neutron_files : list of str or openmc.data.endf.Evaluation
neutron_files : list of str, openmc.data.endf.Evaluation, or endf.Material
List of ENDF neutron reaction sub-library files
reactions : iterable of str, optional
Transmutation reactions to include in the depletion chain, e.g.,
@ -363,7 +363,7 @@ class Chain:
print('Processing neutron sub-library files...')
reactions = {}
for f in neutron_files:
evaluation = openmc.data.endf.Evaluation(f)
evaluation = openmc.data.endf.as_evaluation(f)
name = evaluation.gnds_name
reactions[name] = {}
for mf, mt, nc, mod in evaluation.reaction_list:

11
openmc/deplete/d1s.py
View file

@ -186,8 +186,6 @@ def apply_time_correction(
# Apply TCF, broadcasting to the correct dimensions
tcf.shape = (1, -1, 1, 1, 1)
new_tally._sum = tally_sum * tcf
new_tally._sum_sq = tally_sum_sq * (tcf*tcf)
new_tally._mean = tally_mean * tcf
new_tally._std_dev = tally_std_dev * tcf
@ -196,6 +194,8 @@ def apply_time_correction(
if sum_nuclides:
# Sum over parent nuclides (note that when combining different bins for
# parent nuclide, we can't work directly on sum_sq)
new_tally._sum = None
new_tally._sum_sq = None
new_tally._mean = new_tally.mean.sum(axis=1).reshape(shape)
new_tally._std_dev = np.linalg.norm(new_tally.std_dev, axis=1).reshape(shape)
new_tally._derived = True
@ -203,9 +203,10 @@ def apply_time_correction(
# Remove ParentNuclideFilter
new_tally.filters.pop(i_filter)
else:
# Change shape back to (filter combinations, nuclides, scores)
new_tally._sum.shape = shape
new_tally._sum_sq.shape = shape
# Apply TCF and change shape back to (filter combinations, nuclides,
# scores)
new_tally._sum = (tally_sum * tcf).reshape(shape)
new_tally._sum_sq = (tally_sum_sq * (tcf*tcf)).reshape(shape)
new_tally._mean.shape = shape
new_tally._std_dev.shape = shape

4
openmc/deplete/independent_operator.py
View file

@ -339,8 +339,12 @@ class IndependentOperator(OpenMCOperator):
for i_nuc in nuc_index:
nuc = self.nuc_ind_map[i_nuc]
if nuc not in xs._index_nuc:
continue
for i_rx in react_index:
rx = self.rx_ind_map[i_rx]
if rx not in xs._index_rx:
continue
# Determine reaction rate by multiplying xs in [b] by flux
# in [n-cm/src] to give [(reactions/src)*b-cm/atom]

63
openmc/deplete/microxs.py
View file

@ -84,7 +84,12 @@ def get_microxs_and_flux(
reactions listed in the depletion chain file are used.
energies : iterable of float or str
Energy group boundaries in [eV] or the name of the group structure.
If left as None energies will default to [0.0, 100e6]
If left as None, no energy filter is applied to the flux tally. When
`reaction_rate_mode` is "direct", these boundaries define the output
flux and microscopic cross section energy group structure. When
`reaction_rate_mode` is "flux", these boundaries define the multigroup
flux tally used to collapse continuous-energy cross sections; returned
fluxes and microscopic cross sections are one-group.
reaction_rate_mode : {"direct", "flux"}, optional
The "direct" method tallies reaction rates directly (per energy
group). The "flux" method tallies a multigroup flux spectrum and then
@ -110,7 +115,9 @@ def get_microxs_and_flux(
reaction_rate_opts : dict, optional
When `reaction_rate_mode="flux"`, allows selecting a subset of
nuclide/reaction pairs to be computed via direct reaction-rate tallies
(per energy group). Supported keys: "nuclides", "reactions".
over one energy bin spanning the full `energies` range. Supported keys:
"nuclides", "reactions". If "reactions" are specified without
"nuclides", all selected nuclides are used.
Returns
-------
@ -139,10 +146,14 @@ def get_microxs_and_flux(
nuclides = [nuc.name for nuc in chain.nuclides
if nuc.name in nuclides_with_data]
# Set up the reaction rate and flux tallies
# Set up the reaction rate and flux tallies. When energies are omitted, no
# energy filter is needed for the transport calculation. A one-group energy
# range is still needed later if flux collapse is requested.
collapse_energies = energies
if energies is None:
energies = [0.0, 100.0e6]
if isinstance(energies, str):
energy_filter = None
collapse_energies = [0.0, 100.0e6]
elif isinstance(energies, str):
energy_filter = openmc.EnergyFilter.from_group_structure(energies)
else:
energy_filter = openmc.EnergyFilter(energies)
@ -172,8 +183,11 @@ def get_microxs_and_flux(
rr_reactions = list(reactions)
elif reaction_rate_mode == 'flux' and reaction_rate_opts:
opts = reaction_rate_opts or {}
rr_nuclides = list(opts.get('nuclides', []))
rr_reactions = list(opts.get('reactions', []))
if rr_reactions:
rr_nuclides = list(opts.get('nuclides', nuclides))
else:
rr_nuclides = list(opts.get('nuclides', []))
# Keep only requested pairs within overall sets
if rr_nuclides:
rr_nuclides = [n for n in rr_nuclides if n in set(nuclides)]
@ -182,7 +196,7 @@ def get_microxs_and_flux(
# Use 1-group energy filter for RR in flux mode
has_rr = bool(rr_nuclides and rr_reactions)
if has_rr and reaction_rate_mode == 'flux':
if has_rr and reaction_rate_mode == 'flux' and energy_filter is not None:
rr_energy_filter = openmc.EnergyFilter(
[energy_filter.values[0], energy_filter.values[-1]])
else:
@ -194,14 +208,18 @@ def get_microxs_and_flux(
model.tallies = []
for i, domain_filter in enumerate(domain_filters):
flux_tally = openmc.Tally(name=f'MicroXS flux {i}')
flux_tally.filters = [domain_filter, energy_filter]
flux_tally.filters = [domain_filter]
if energy_filter is not None:
flux_tally.filters.append(energy_filter)
flux_tally.scores = ['flux']
model.tallies.append(flux_tally)
flux_tallies.append(flux_tally)
if has_rr:
rr_tally = openmc.Tally(name=f'MicroXS RR {i}')
rr_tally.filters = [domain_filter, rr_energy_filter]
rr_tally.filters = [domain_filter]
if rr_energy_filter is not None:
rr_tally.filters.append(rr_energy_filter)
rr_tally.nuclides = rr_nuclides
rr_tally.multiply_density = False
rr_tally.scores = rr_reactions
@ -255,8 +273,12 @@ def get_microxs_and_flux(
all_flux_arrays = []
for flux_tally in flux_tallies:
# Get flux values and make energy groups last dimension
flux = flux_tally.get_reshaped_data() # (domains, groups, 1, 1)
flux = np.moveaxis(flux, 1, -1) # (domains, 1, 1, groups)
flux = flux_tally.get_reshaped_data()
if energy_filter is None:
flux = flux[..., np.newaxis] # (domains, 1, 1, groups)
else:
# (domains, groups, 1, 1) -> (domains, 1, 1, groups)
flux = np.moveaxis(flux, 1, -1)
all_flux_arrays.append(flux)
fluxes.extend(flux.squeeze((1, 2)))
@ -266,8 +288,15 @@ def get_microxs_and_flux(
for flux_arr, rr_tally in zip(all_flux_arrays, rr_tallies):
flux = flux_arr
# Get reaction rates and make energy groups last dimension
reaction_rates = rr_tally.get_reshaped_data() # (domains, groups, nuclides, reactions)
reaction_rates = np.moveaxis(reaction_rates, 1, -1) # (domains, nuclides, reactions, groups)
reaction_rates = rr_tally.get_reshaped_data()
if rr_energy_filter is None:
# (domains, nuclides, reactions) ->
# (domains, nuclides, reactions, groups)
reaction_rates = reaction_rates[..., np.newaxis]
else:
# (domains, groups, nuclides, reactions) ->
# (domains, nuclides, reactions, groups)
reaction_rates = np.moveaxis(reaction_rates, 1, -1)
# If RR is 1-group, sum flux over groups
if reaction_rate_mode == "flux":
@ -279,16 +308,20 @@ def get_microxs_and_flux(
direct_micros.extend(
MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs)
# If using flux mode, compute flux-collapsed microscopic XS
if reaction_rate_mode == 'flux':
# Compute flux-collapsed microscopic XS
flux_micros = [MicroXS.from_multigroup_flux(
energies=energies,
energies=collapse_energies,
multigroup_flux=flux_i,
chain_file=chain_file,
nuclides=nuclides,
reactions=reactions
) for flux_i in fluxes]
# We need to return one-group fluxes to match the microscopic cross
# sections, which are always one-group by virtue of the collapse
fluxes = [flux.sum(keepdims=True) for flux in fluxes]
# Decide which micros to use and merge if needed
if reaction_rate_mode == 'flux' and rr_tallies:
micros = [m1.merge(m2) for m1, m2 in zip(flux_micros, direct_micros)]

244
openmc/deplete/r2s.py
View file

@ -1,5 +1,6 @@
from __future__ import annotations
from collections.abc import Sequence
from contextlib import nullcontext
import copy
from datetime import datetime
import json
@ -8,6 +9,7 @@ from pathlib import Path
import numpy as np
import openmc
from . import IndependentOperator, PredictorIntegrator
from .chain import Chain
from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5
from .results import Results
from ..checkvalue import PathLike
@ -149,7 +151,7 @@ class R2SManager:
photon_time_indices: Sequence[int] | None = None,
output_dir: PathLike | None = None,
bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
chain_file: PathLike | None = None,
chain_file: PathLike | Chain | None = None,
micro_kwargs: dict | None = None,
mat_vol_kwargs: dict | None = None,
run_kwargs: dict | None = None,
@ -187,9 +189,10 @@ class R2SManager:
Dictionary mapping cell IDs to bounding boxes used for spatial
source sampling in cell-based R2S calculations. Required if method
is 'cell-based'.
chain_file : PathLike, optional
Path to the depletion chain XML file to use during activation. If
not provided, the default configured chain file will be used.
chain_file : PathLike or openmc.deplete.Chain, optional
Path to the depletion chain XML file or depletion chain object to
use during activation. If not provided, the default configured
chain file will be used.
micro_kwargs : dict, optional
Additional keyword arguments passed to
:func:`openmc.deplete.get_microxs_and_flux` during the neutron
@ -216,6 +219,8 @@ class R2SManager:
# consistency (different ranks may have slightly different times)
stamp = datetime.now().strftime('%Y-%m-%dT%H-%M-%S')
output_dir = Path(comm.bcast(f'r2s_{stamp}'))
else:
output_dir = Path(output_dir)
# Set run_kwargs for the neutron transport step
if micro_kwargs is None:
@ -226,22 +231,35 @@ class R2SManager:
operator_kwargs = {}
run_kwargs.setdefault('output', False)
micro_kwargs.setdefault('run_kwargs', run_kwargs)
# If a chain file is provided, prefer it for steps 1 and 2
if chain_file is not None:
micro_kwargs.setdefault('chain_file', chain_file)
operator_kwargs.setdefault('chain_file', chain_file)
self.step1_neutron_transport(
output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs
)
self.step2_activation(
timesteps, source_rates, timestep_units, output_dir / 'activation',
operator_kwargs=operator_kwargs
)
self.step3_photon_transport(
photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
# DecaySpectrum distributions are resolved in the C++ solver using
# OPENMC_CHAIN_FILE. If a Chain object was passed, write an XML
# representation alongside the R2S outputs.
if isinstance(chain_file, Chain):
output_dir.mkdir(parents=True, exist_ok=True)
chain_path = output_dir / 'chain.xml'
if comm.rank == 0:
chain_file.export_to_xml(chain_path)
comm.barrier()
else:
chain_path = chain_file
chain_context = (
openmc.config.patch('chain_file', chain_path)
if chain_path is not None else nullcontext()
)
with chain_context:
self.step1_neutron_transport(
output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs
)
self.step2_activation(
timesteps, source_rates, timestep_units,
output_dir / 'activation', operator_kwargs=operator_kwargs
)
self.step3_photon_transport(
photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
)
return output_dir
@ -516,45 +534,30 @@ class R2SManager:
if different_photon_model:
photon_cells = self.photon_model.geometry.get_all_cells()
for time_index in time_indices:
# Create decay photon source
if self.method == 'mesh-based':
self.photon_model.settings.source = \
self.get_decay_photon_source_mesh(time_index)
else:
sources = []
results = self.results['depletion_results']
for cell, original_mat in zip(self.domains, self.results['activation_materials']):
# Skip if the cell is not in the photon model or the
# material has changed
if different_photon_model:
if cell.id not in photon_cells or \
# Determine eligible work items upfront (independent of time index).
if self.method == 'mesh-based':
work_items = self._get_mesh_work_items()
else:
work_items = []
for cell, original_mat in zip(
self.domains, self.results['activation_materials']):
if different_photon_model:
if cell.id not in photon_cells or \
cell.fill.id != photon_cells[cell.id].fill.id:
continue
continue
work_items.append((cell, original_mat, bounding_boxes[cell.id]))
# Get bounding box for the cell
bounding_box = bounding_boxes[cell.id]
# Get activated material composition
activated_mat = results[time_index].get_material(str(original_mat.id))
# Create decay photon source source
space = openmc.stats.Box(*bounding_box)
energy = activated_mat.get_decay_photon_energy()
strength = energy.integral() if energy is not None else 0.0
source = openmc.IndependentSource(
space=space,
energy=energy,
particle='photon',
strength=strength,
constraints={'domains': [cell]}
)
sources.append(source)
self.photon_model.settings.source = sources
# Ensure photon transport is enabled in settings
self.photon_model.settings.photon_transport = True
for time_index in time_indices:
# Convert time_index (which may be negative) to a normal index
if time_index < 0:
time_index = len(self.results['depletion_results']) + time_index
time_index += len(self.results['depletion_results'])
# Build decay photon sources and assign to the photon model
sources = self._create_photon_sources(time_index, work_items)
self.photon_model.settings.source = sources
# Run photon transport calculation
photon_dir = Path(output_dir) / f'time_{time_index}'
@ -567,58 +570,30 @@ class R2SManager:
sp.tallies[tally.id] for tally in self.photon_model.tallies
]
def get_decay_photon_source_mesh(
self,
time_index: int = -1
) -> list[openmc.IndependentSource]:
"""Create decay photon source for a mesh-based calculation.
def _get_mesh_work_items(self):
"""Enumerate mesh-based work items across all meshes.
For each mesh element-material combination across all meshes, an
:class:`~openmc.IndependentSource` is created with a
:class:`~openmc.stats.Box` spatial distribution based on the bounding
box of the material within the mesh element. A material constraint is
also applied so that sampled source sites are limited to the correct
region.
When the photon transport model is different from the neutron model, the
photon MeshMaterialVolumes is used to determine whether an (element,
material) combination exists in the photon model.
Parameters
----------
time_index : int, optional
Time index for the decay photon source. Default is -1 (last time).
Returns a list of (index_mat, mat_id, bbox) tuples for each eligible
mesh element--material combination, where index_mat is the index into
the activation materials list, mat_id is the material ID, and bbox is
the bounding box for that mesh element--material combination.
Returns
-------
list of openmc.IndependentSource
A list of IndependentSource objects for the decay photons, one for
each mesh element-material combination with non-zero source strength.
list of tuple
Each tuple is (index_mat, mat_id, bbox).
"""
mat_dict = self.neutron_model._get_all_materials()
# List to hold all sources
sources = []
# Index in the overall list of activated materials
index_mat = 0
# Get various results from previous steps
mmv_list = self.results['mesh_material_volumes']
materials = self.results['activation_materials']
results = self.results['depletion_results']
photon_mmv_list = self.results.get('mesh_material_volumes_photon')
work_items = []
index_mat = 0
for mesh_idx, mat_vols in enumerate(mmv_list):
photon_mat_vols = photon_mmv_list[mesh_idx] \
if photon_mmv_list is not None else None
# Total number of mesh elements for this mesh
n_elements = mat_vols.num_elements
for index_elem in range(n_elements):
# Determine which materials exist in the photon model for this element
if photon_mat_vols is not None:
photon_materials = {
mat_id
@ -626,36 +601,77 @@ class R2SManager:
if mat_id is not None
}
for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
# Skip void volume
for mat_id, _, bbox in mat_vols.by_element(
index_elem, include_bboxes=True):
if mat_id is None:
continue
# Skip if this material doesn't exist in photon model
if photon_mat_vols is not None and mat_id not in photon_materials:
if photon_mat_vols is not None \
and mat_id not in photon_materials:
index_mat += 1
continue
# Get activated material composition
original_mat = materials[index_mat]
activated_mat = results[time_index].get_material(str(original_mat.id))
# Create decay photon source
energy = activated_mat.get_decay_photon_energy()
if energy is not None:
strength = energy.integral()
space = openmc.stats.Box(*bbox)
sources.append(openmc.IndependentSource(
space=space,
energy=energy,
particle='photon',
strength=strength,
constraints={'domains': [mat_dict[mat_id]]}
))
# Increment index of activated material
work_items.append((index_mat, mat_id, bbox))
index_mat += 1
return work_items
def _create_photon_sources(self, time_index, work_items):
"""Create decay photon sources for a set of regions.
Builds :class:`openmc.IndependentSource` objects with
:class:`openmc.stats.DecaySpectrum` energy distributions that will be
serialized to XML and resolved against the depletion chain by the C++
solver.
Parameters
----------
time_index : int
Index into depletion results.
work_items : list of tuple
For mesh-based: list of (index_mat, mat_id, bbox).
For cell-based: list of (cell, original_mat, bbox).
Returns
-------
list of openmc.IndependentSource
Photon sources for each activated region.
"""
step_result = self.results['depletion_results'][time_index]
materials = self.results['activation_materials']
mesh_based = self.method == 'mesh-based'
if mesh_based:
mat_dict = self.neutron_model._get_all_materials()
sources = []
for item in work_items:
if mesh_based:
index_mat, domain_id, bbox = item
original_mat = materials[index_mat]
domain = mat_dict[domain_id]
else:
cell, original_mat, bbox = item
domain = cell
activated_mat = step_result.get_material(str(original_mat.id))
nuclides = activated_mat.get_nuclide_atom_densities()
if not nuclides:
continue
# Eliminate nuclides with zero density
nuclides = {nuclide: density for nuclide, density in nuclides.items()
if density > 0}
energy = openmc.stats.DecaySpectrum(nuclides, activated_mat.volume)
energy.clip(inplace=True)
if not energy.nuclides:
continue
sources.append(openmc.IndependentSource(
space=openmc.stats.Box(bbox.lower_left, bbox.upper_right),
energy=energy,
particle='photon',
constraints={'domains': [domain]},
))
return sources
def load_results(self, path: PathLike):

4
openmc/examples.py
View file

@ -150,7 +150,7 @@ def pwr_core() -> openmc.Model:
rpv_steel.add_nuclide('Ni60', 0.0026776, 'wo')
rpv_steel.add_nuclide('Mn55', 0.01, 'wo')
rpv_steel.add_nuclide('Cr52', 0.002092475, 'wo')
rpv_steel.add_nuclide('C0', 0.0025, 'wo')
rpv_steel.add_element('C', 0.0025, 'wo')
rpv_steel.add_nuclide('Cu63', 0.0013696, 'wo')
lower_rad_ref = openmc.Material(6, name='Lower radial reflector')
@ -1618,4 +1618,4 @@ def random_ray_three_region_cube_with_detectors() -> openmc.Model:
model.settings = settings
model.tallies = tallies
return model
return model

5
openmc/lib/core.py
View file

@ -32,7 +32,10 @@ class _SourceSite(Structure):
('particle', c_int32),
('parent_nuclide', c_int),
('parent_id', c_int64),
('progeny_id', c_int64)]
('progeny_id', c_int64),
('wgt_born', c_double),
('wgt_ww_born', c_double),
('n_split', c_int64)]
# Define input type for numpy arrays that will be passed into C++ functions
# Must be an int or double array, with single dimension that is contiguous

358
openmc/lib/plot.py
View file

@ -9,6 +9,7 @@ from .core import _FortranObjectWithID
from .error import _error_handler
import numpy as np
import warnings
class _Position(Structure):
@ -51,218 +52,209 @@ class _Position(Structure):
return f"({self.x}, {self.y}, {self.z})"
class _PlotBase(Structure):
"""A structure defining a 2-D geometry slice with underlying c-types
def _extract_slice_data_args(plot):
"""Convert a legacy plot-like object into slice_data keyword arguments."""
try:
kwargs = {
'origin': tuple(plot.origin),
'width': (plot.width, plot.height),
'basis': plot.basis,
'pixels': (plot.h_res, plot.v_res),
'show_overlaps': getattr(plot, 'color_overlaps', False),
'level': getattr(plot, 'level', -1),
}
except AttributeError as exc:
raise TypeError(
"plot must be a legacy plot-like object with origin, width, "
"height, basis, h_res, and v_res attributes."
) from exc
return kwargs
C-Type Attributes
-----------------
origin_ : openmc.lib.plot._Position
A position defining the origin of the plot.
width_ : openmc.lib.plot._Position
The width of the plot along the x, y, and z axes, respectively
basis_ : c_int
The axes basis of the plot view.
pixels_ : c_size_t[3]
The resolution of the plot in the horizontal and vertical dimensions
color_overlaps_ : c_bool
Whether to assign unique IDs (-3) to overlapping regions.
level_ : c_int
The universe level for the plot view
Attributes
_dll.openmc_slice_data.argtypes = [
POINTER(c_double * 3), # origin
POINTER(c_double * 3), # u_span
POINTER(c_double * 3), # v_span
POINTER(c_size_t * 2), # pixels
c_bool, # show_overlaps
c_int, # level
c_int32, # filter_index
POINTER(c_int32), # geom_data
POINTER(c_double), # property_data (can be None)
]
_dll.openmc_slice_data.restype = c_int
_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,
include_properties=True):
"""Generate a 2D raster of geometry and property data for plotting.
Parameters
----------
origin : tuple or list of ndarray
Origin (center) of the plot
width : float
The horizontal dimension of the plot in geometry units (cm)
height : float
The vertical dimension of the plot in geometry units (cm)
basis : string
One of {'xy', 'xz', 'yz'} indicating the horizontal and vertical
axes of the plot.
h_res : int
The horizontal resolution of the plot in pixels
v_res : int
The vertical resolution of the plot in pixels
level : int
The universe level for the plot (default: -1 -> all universes shown)
origin : sequence of float
Center position of the plot [x, y, z]
width : sequence of float
Width of the plot [horizontal, vertical]. Mutually exclusive with
u_span/v_span.
basis : {'xy', 'xz', 'yz'} or int
Plot basis. Ignored if u_span/v_span are provided.
u_span : sequence of float, optional
Full-width span vector for the horizontal axis (3 values). Mutually
exclusive with width.
v_span : sequence of float, optional
Full-height span vector for the vertical axis (3 values). Mutually
exclusive with width.
pixels : sequence of int
Number of pixels [horizontal, vertical]
show_overlaps : bool, optional
Whether to detect overlapping cells
level : int, optional
Universe level (-1 for deepest)
filter : openmc.lib.Filter, optional
Filter for bin index lookup
include_properties : bool, optional
Whether to compute temperature/density
Returns
-------
geom_data : numpy.ndarray
Array of shape (v_res, h_res, 3) or (v_res, h_res, 4) with int32 dtype.
Contains [cell_id, cell_instance, material_id] when no filter is provided,
or [cell_id, cell_instance, material_id, filter_bin] when a filter is provided.
property_data : numpy.ndarray or None
Array of shape (v_res, h_res, 2) with float64 dtype containing
[temperature, density], or None if include_properties=False
"""
_fields_ = [('origin_', _Position),
('width_', _Position),
('basis_', c_int),
('pixels_', 3*c_size_t),
('color_overlaps_', c_bool),
('level_', c_int)]
if pixels is None:
raise ValueError("pixels must be specified.")
if len(pixels) != 2:
raise ValueError("pixels must be a length-2 sequence.")
def __init__(self):
self.level_ = -1
self.basis_ = 1
self.color_overlaps_ = False
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.")
@property
def origin(self):
return self.origin_
@origin.setter
def origin(self, origin):
self.origin_.x = origin[0]
self.origin_.y = origin[1]
self.origin_.z = origin[2]
@property
def width(self):
return self.width_.x
@width.setter
def width(self, width):
self.width_.x = width
@property
def height(self):
return self.width_.y
@height.setter
def height(self, height):
self.width_.y = height
@property
def basis(self):
if self.basis_ == 1:
return 'xy'
elif self.basis_ == 2:
return 'xz'
elif self.basis_ == 3:
return 'yz'
raise ValueError(f"Plot basis {self.basis_} is invalid")
@basis.setter
def basis(self, basis):
if u_span is not None or v_span is not None:
if u_span is None or v_span is None:
raise ValueError("Both u_span and v_span must be provided.")
u_span = np.asarray(u_span, dtype=float)
v_span = np.asarray(v_span, dtype=float)
if u_span.shape != (3,) or v_span.shape != (3,):
raise ValueError("u_span and v_span must be length-3 sequences.")
u_norm = np.linalg.norm(u_span)
v_norm = np.linalg.norm(v_span)
if u_norm == 0.0 or v_norm == 0.0:
raise ValueError("u_span and v_span must be non-zero vectors.")
dot = float(np.dot(u_span, v_span))
ortho_tol = 1.0e-10 * u_norm * v_norm
if abs(dot) > ortho_tol:
raise ValueError("u_span and v_span must be orthogonal.")
else:
if width is None:
raise ValueError("width must be provided when u_span/v_span are not set.")
if len(width) != 2:
raise ValueError("width must be a length-2 sequence.")
basis_map = {'xy': 1, 'xz': 2, 'yz': 3}
if isinstance(basis, str):
valid_bases = ('xy', 'xz', 'yz')
basis = basis.lower()
if basis not in valid_bases:
if basis not in basis_map:
raise ValueError(f"{basis} is not a valid plot basis.")
if basis == 'xy':
self.basis_ = 1
elif basis == 'xz':
self.basis_ = 2
elif basis == 'yz':
self.basis_ = 3
return
if isinstance(basis, int):
valid_bases = (1, 2, 3)
if basis not in valid_bases:
basis = basis_map[basis]
elif isinstance(basis, int):
if basis not in basis_map.values():
raise ValueError(f"{basis} is not a valid plot basis.")
self.basis_ = basis
return
else:
raise ValueError(f"{basis} is not a valid plot basis.")
raise ValueError(f"{basis} of type {type(basis)} is an invalid plot basis")
if basis == 1:
u_span = np.array([width[0], 0.0, 0.0], dtype=float)
v_span = np.array([0.0, width[1], 0.0], dtype=float)
elif basis == 2:
u_span = np.array([width[0], 0.0, 0.0], dtype=float)
v_span = np.array([0.0, 0.0, width[1]], dtype=float)
else:
u_span = np.array([0.0, width[0], 0.0], dtype=float)
v_span = np.array([0.0, 0.0, width[1]], dtype=float)
@property
def h_res(self):
return self.pixels_[0]
origin = np.asarray(origin, dtype=float)
if origin.shape != (3,):
raise ValueError("origin must be a length-3 sequence.")
@h_res.setter
def h_res(self, h_res):
self.pixels_[0] = h_res
# Prepare ctypes arrays
origin_arr = (c_double * 3)(*origin)
u_span_arr = (c_double * 3)(*u_span)
v_span_arr = (c_double * 3)(*v_span)
pixels_arr = (c_size_t * 2)(*pixels)
@property
def v_res(self):
return self.pixels_[1]
# Get internal filter index from filter ID if filter is provided
if filter is not None:
filter_index = c_int32()
_dll.openmc_get_filter_index(filter.id, filter_index)
filter_index = filter_index.value
else:
filter_index = -1
@v_res.setter
def v_res(self, v_res):
self.pixels_[1] = v_res
# Allocate output arrays with dynamic size based on filter
n_geom_fields = 4 if filter is not None else 3
geom_data = np.zeros((pixels[1], pixels[0], n_geom_fields), dtype=np.int32)
if include_properties:
property_data = np.zeros((pixels[1], pixels[0], 2), dtype=np.float64)
prop_ptr = property_data.ctypes.data_as(POINTER(c_double))
else:
property_data = None
prop_ptr = None
@property
def level(self):
return int(self.level_)
_dll.openmc_slice_data(
origin_arr,
u_span_arr,
v_span_arr,
pixels_arr,
show_overlaps,
level,
filter_index,
geom_data.ctypes.data_as(POINTER(c_int32)),
prop_ptr
)
@level.setter
def level(self, level):
self.level_ = level
@property
def color_overlaps(self):
return self.color_overlaps_
@color_overlaps.setter
def color_overlaps(self, color_overlaps):
self.color_overlaps_ = color_overlaps
def __repr__(self):
out_str = ["-----",
"Plot:",
"-----",
f"Origin: {self.origin}",
f"Width: {self.width}",
f"Height: {self.height}",
f"Basis: {self.basis}",
f"HRes: {self.h_res}",
f"VRes: {self.v_res}",
f"Color Overlaps: {self.color_overlaps}",
f"Level: {self.level}"]
return '\n'.join(out_str)
_dll.openmc_id_map.argtypes = [POINTER(_PlotBase), POINTER(c_int32)]
_dll.openmc_id_map.restype = c_int
_dll.openmc_id_map.errcheck = _error_handler
return geom_data, property_data
def id_map(plot):
"""Deprecated compatibility wrapper for geometry ID maps.
This function is kept for compatibility and will be removed in a future
release. Use `slice_data(..., include_properties=False)` instead.
"""
Generate a 2-D map of cell and material IDs. Used for in-memory image
generation.
warnings.warn(
"openmc.lib.id_map is deprecated and will be removed in a future "
"release; use openmc.lib.slice_data(..., include_properties=False).",
FutureWarning,
)
Parameters
----------
plot : openmc.lib.plot._PlotBase
Object describing the slice of the model to be generated
Returns
-------
id_map : numpy.ndarray
A NumPy array with shape (vertical pixels, horizontal pixels, 3) of
OpenMC property ids with dtype int32. The last dimension of the array
contains, in order, cell IDs, cell instances, and material IDs.
"""
img_data = np.zeros((plot.v_res, plot.h_res, 3),
dtype=np.dtype('int32'))
_dll.openmc_id_map(plot, img_data.ctypes.data_as(POINTER(c_int32)))
return img_data
_dll.openmc_property_map.argtypes = [POINTER(_PlotBase), POINTER(c_double)]
_dll.openmc_property_map.restype = c_int
_dll.openmc_property_map.errcheck = _error_handler
kwargs = _extract_slice_data_args(plot)
geom_data, _ = slice_data(include_properties=False, **kwargs)
return geom_data[:, :, :3]
def property_map(plot):
"""Deprecated compatibility wrapper for temperature/density maps.
This function is kept for compatibility and will be removed in a future
release. Use `slice_data(..., include_properties=True)` instead.
"""
Generate a 2-D map of cell temperatures and material densities. Used for
in-memory image generation.
warnings.warn(
"openmc.lib.property_map is deprecated and will be removed in a "
"future release; use openmc.lib.slice_data(..., "
"include_properties=True).",
FutureWarning,
)
Parameters
----------
plot : openmc.lib.plot._PlotBase
Object describing the slice of the model to be generated
Returns
-------
property_map : numpy.ndarray
A NumPy array with shape (vertical pixels, horizontal pixels, 2) of
OpenMC property ids with dtype float
"""
prop_data = np.zeros((plot.v_res, plot.h_res, 2))
_dll.openmc_property_map(plot, prop_data.ctypes.data_as(POINTER(c_double)))
kwargs = _extract_slice_data_args(plot)
_, prop_data = slice_data(include_properties=True, **kwargs)
return prop_data
_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

18
openmc/material.py
View file

@ -296,6 +296,8 @@ class Material(IDManagerMixin):
mass += nuc.percent
# Compute and return the molar mass
if moles == 0.0:
raise ValueError("Material has no nuclides; cannot compute molar mass")
return mass / moles
@property
@ -1418,6 +1420,9 @@ class Material(IDManagerMixin):
if volume is None:
volume = self.volume
if units in {'Bq', 'Ci'} and volume is None:
raise ValueError(f"Volume must be set in order to compute activity in '{units}'.")
if units == 'Bq':
multiplier = volume
elif units == 'Bq/cm3':
@ -1474,6 +1479,8 @@ class Material(IDManagerMixin):
if units == 'W':
multiplier = volume if volume is not None else self.volume
if multiplier is None:
raise ValueError("Volume must be set in order to compute total decay heat.")
elif units == 'W/cm3':
multiplier = 1
elif units == 'W/m3':
@ -2282,7 +2289,7 @@ class Materials(cv.CheckedList):
multigroup_fluxes: Sequence[Sequence[float]]
Energy-dependent multigroup flux values, where each sublist corresponds
to a specific material. Will be normalized so that it sums to 1.
energy_group_structures': Sequence[Sequence[float] | str]
energy_group_structures: Sequence[Sequence[float] | str]
Energy group boundaries in [eV] or the name of the group structure.
timesteps : iterable of float or iterable of tuple
Array of timesteps. Note that values are not cumulative. The units are
@ -2317,6 +2324,11 @@ class Materials(cv.CheckedList):
for mat in self:
mat.depletable = True
if len(multigroup_fluxes) != len(self):
raise ValueError("multigroup_fluxes length must match number of materials")
if len(energy_group_structures) != len(self):
raise ValueError("energy_group_structures length must match number of materials")
chain = _get_chain(chain_file)
# Create MicroXS objects for all materials
@ -2327,6 +2339,10 @@ class Materials(cv.CheckedList):
for material, flux, energy in zip(
self, multigroup_fluxes, energy_group_structures
):
if material.volume is None:
raise ValueError(
f"Material {material.id} has no volume; cannot deplete"
)
temperature = material.temperature or 293.6
micro_xs = openmc.deplete.MicroXS.from_multigroup_flux(
energies=energy,

225
openmc/mesh.py
View file

@ -3,7 +3,7 @@ import warnings
from abc import ABC, abstractmethod
from collections.abc import Iterable, Sequence, Mapping
from functools import wraps
from math import pi, sqrt, atan2
import math
from numbers import Integral, Real
from pathlib import Path
from typing import Protocol
@ -11,12 +11,10 @@ from typing import Protocol
import h5py
import lxml.etree as ET
import numpy as np
from pathlib import Path
import openmc
import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from openmc.utility_funcs import change_directory
from .bounding_box import BoundingBox
from ._xml import get_elem_list, get_text
from .mixin import IDManagerMixin
@ -291,7 +289,7 @@ class MeshBase(IDManagerMixin, ABC):
"""
mesh_type = 'regular' if 'type' not in group.keys() else group['type'][()].decode()
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
mesh_name = '' if not 'name' in group else group['name'][()].decode()
mesh_name = '' if 'name' not in group else group['name'][()].decode()
if mesh_type == 'regular':
return RegularMesh.from_hdf5(group, mesh_id, mesh_name)
@ -1903,7 +1901,7 @@ class CylindricalMesh(StructuredMesh):
self,
r_grid: Sequence[float],
z_grid: Sequence[float],
phi_grid: Sequence[float] = (0, 2*pi),
phi_grid: Sequence[float] = (0, 2*math.pi),
origin: Sequence[float] = (0., 0., 0.),
mesh_id: int | None = None,
name: str = '',
@ -1960,7 +1958,7 @@ class CylindricalMesh(StructuredMesh):
cv.check_length('mesh phi_grid', grid, 2)
cv.check_increasing('mesh phi_grid', grid)
grid = np.asarray(grid, dtype=float)
if np.any((grid < 0.0) | (grid > 2*pi)):
if np.any((grid < 0.0) | (grid > 2*math.pi)):
raise ValueError("phi_grid values must be in [0, 2π].")
self._phi_grid = grid
@ -2028,9 +2026,9 @@ class CylindricalMesh(StructuredMesh):
return string
def get_indices_at_coords(
self,
coords: Sequence[float]
) -> tuple[int, int, int]:
self,
coords: Sequence[float]
) -> tuple[int, int, int]:
"""Finds the index of the mesh element at the specified coordinates.
.. versionadded:: 0.15.0
@ -2046,7 +2044,7 @@ class CylindricalMesh(StructuredMesh):
The r, phi, z indices
"""
r_value_from_origin = sqrt((coords[0]-self.origin[0])**2 + (coords[1]-self.origin[1])**2)
r_value_from_origin = math.hypot(coords[0]-self.origin[0], coords[1]-self.origin[1])
if r_value_from_origin < self.r_grid[0] or r_value_from_origin > self.r_grid[-1]:
raise ValueError(
@ -2070,13 +2068,13 @@ class CylindricalMesh(StructuredMesh):
delta_x = coords[0] - self.origin[0]
delta_y = coords[1] - self.origin[1]
# atan2 returns values in -pi to +pi range
phi_value = atan2(delta_y, delta_x)
phi_value = math.atan2(delta_y, delta_x)
if delta_x < 0 and delta_y < 0:
# returned phi_value anticlockwise and negative
phi_value += 2 * pi
phi_value += 2 * math.pi
if delta_x > 0 and delta_y < 0:
# returned phi_value anticlockwise and negative
phi_value += 2 * pi
phi_value += 2 * math.pi
phi_grid_values = np.array(self.phi_grid)
@ -2111,7 +2109,7 @@ class CylindricalMesh(StructuredMesh):
dimension: Sequence[int] = (10, 10, 10),
mesh_id: int | None = None,
name: str = '',
phi_grid_bounds: Sequence[float] = (0.0, 2*pi),
phi_grid_bounds: Sequence[float] = (0.0, 2*math.pi),
enclose_domain: bool = False,
) -> CylindricalMesh:
"""Create CylindricalMesh from a bounding box.
@ -2339,8 +2337,8 @@ class SphericalMesh(StructuredMesh):
def __init__(
self,
r_grid: Sequence[float],
phi_grid: Sequence[float] = (0, 2*pi),
theta_grid: Sequence[float] = (0, pi),
phi_grid: Sequence[float] = (0, 2*math.pi),
theta_grid: Sequence[float] = (0, math.pi),
origin: Sequence[float] = (0., 0., 0.),
mesh_id: int | None = None,
name: str = '',
@ -2397,7 +2395,7 @@ class SphericalMesh(StructuredMesh):
cv.check_length('mesh theta_grid', grid, 2)
cv.check_increasing('mesh theta_grid', grid)
grid = np.asarray(grid, dtype=float)
if np.any((grid < 0.0) | (grid > pi)):
if np.any((grid < 0.0) | (grid > math.pi)):
raise ValueError("theta_grid values must be in [0, π].")
self._theta_grid = grid
@ -2411,7 +2409,7 @@ class SphericalMesh(StructuredMesh):
cv.check_length('mesh phi_grid', grid, 2)
cv.check_increasing('mesh phi_grid', grid)
grid = np.asarray(grid, dtype=float)
if np.any((grid < 0.0) | (grid > 2*pi)):
if np.any((grid < 0.0) | (grid > 2*math.pi)):
raise ValueError("phi_grid values must be in [0, 2π].")
self._phi_grid = grid
@ -2483,8 +2481,8 @@ class SphericalMesh(StructuredMesh):
dimension: Sequence[int] = (10, 10, 10),
mesh_id: int | None = None,
name: str = '',
phi_grid_bounds: Sequence[float] = (0.0, 2*pi),
theta_grid_bounds: Sequence[float] = (0.0, pi),
phi_grid_bounds: Sequence[float] = (0.0, 2*math.pi),
theta_grid_bounds: Sequence[float] = (0.0, math.pi),
enclose_domain: bool = False,
) -> SphericalMesh:
"""Create SphericalMesh from a bounding box.
@ -2645,10 +2643,82 @@ class SphericalMesh(StructuredMesh):
arr[..., 2] = z + origin[2]
return arr
def get_indices_at_coords(self, coords: Sequence[float]) -> tuple:
raise NotImplementedError(
"get_indices_at_coords is not yet implemented for SphericalMesh"
)
def get_indices_at_coords(
self,
coords: Sequence[float]
) -> tuple[int, int, int]:
"""Find the mesh cell indices containing the specified coordinates.
.. versionadded:: 0.15.4
Parameters
----------
coords : Sequence[float]
Cartesian coordinates of the point as (x, y, z).
Returns
-------
tuple[int, int, int]
The r, theta, phi indices.
Raises
------
ValueError
If the coordinates fall outside the mesh grid boundaries.
"""
dx = coords[0] - self.origin[0]
dy = coords[1] - self.origin[1]
dz = coords[2] - self.origin[2]
r_value = math.hypot(dx, dy, dz)
if r_value < self.r_grid[0] or r_value > self.r_grid[-1]:
raise ValueError(
f'The r value {r_value} computed from the specified '
f'coordinates is outside the r grid range '
f'[{self.r_grid[0]}, {self.r_grid[-1]}].'
)
r_index = int(min(
np.searchsorted(self.r_grid, r_value, side='right') - 1,
len(self.r_grid) - 2
))
if r_value == 0.0:
theta_value = 0.0
phi_value = 0.0
else:
theta_value = math.acos(dz / r_value)
phi_value = math.atan2(dy, dx)
if phi_value < 0:
phi_value += 2 * math.pi
if theta_value < self.theta_grid[0] or theta_value > self.theta_grid[-1]:
raise ValueError(
f'The theta value {theta_value} computed from the specified '
f'coordinates is outside the theta grid range '
f'[{self.theta_grid[0]}, {self.theta_grid[-1]}].'
)
theta_index = int(min(
np.searchsorted(self.theta_grid, theta_value, side='right') - 1,
len(self.theta_grid) - 2
))
if phi_value < self.phi_grid[0] or phi_value > self.phi_grid[-1]:
raise ValueError(
f'The phi value {phi_value} computed from the specified '
f'coordinates is outside the phi grid range '
f'[{self.phi_grid[0]}, {self.phi_grid[-1]}].'
)
phi_index = int(min(
np.searchsorted(self.phi_grid, phi_value, side='right') - 1,
len(self.phi_grid) - 2
))
return (r_index, theta_index, phi_index)
def require_statepoint_data(func):
@ -2738,7 +2808,8 @@ class UnstructuredMesh(MeshBase):
_UNSUPPORTED_ELEM = -1
_LINEAR_TET = 0
_LINEAR_HEX = 1
_VTK_TETRA = 10
_VTK_TET = 10
_VTK_HEX = 12
def __init__(self, filename: PathLike, library: str, mesh_id: int | None = None,
name: str = '', length_multiplier: float = 1.0,
@ -3049,7 +3120,9 @@ class UnstructuredMesh(MeshBase):
n_skipped += 1
continue
else:
raise RuntimeError(f"Invalid element type {elem_type} found")
raise RuntimeError(
f"Invalid element type {elem_type} found in mesh {self.id}"
)
for i, c in enumerate(conn):
if c == -1:
@ -3106,35 +3179,42 @@ class UnstructuredMesh(MeshBase):
datasets: dict | None = None,
volume_normalization: bool = True,
):
def append_dataset(dset, array):
"""Convenience function to append data to an HDF5 dataset"""
origLen = dset.shape[0]
dset.resize(origLen + array.shape[0], axis=0)
dset[origLen:] = array
# This writer supports linear tetrahedra and linear hexahedra elements
conn_list = [] # flattened connectivity ids
cell_sizes = [] # number of points per cell
vtk_types = [] # VTK cell types per cell (uint8)
n_skipped = 0
if self.library != "moab":
raise NotImplementedError("VTKHDF output is only supported for MOAB meshes")
# the self.connectivity contains arrays of length 8 to support hex
# elements as well, in the case of tetrahedra mesh elements, the
# last 4 values are -1 and are removed
trimmed_connectivity = []
for cell in self.connectivity:
# Find the index of the first -1 value, if any
first_negative_index = np.where(cell == -1)[0]
if first_negative_index.size > 0:
# Slice the array up to the first -1 value
trimmed_connectivity.append(cell[: first_negative_index[0]])
for conn, etype in zip(self.connectivity, self.element_types):
if etype == self._LINEAR_TET:
ids = conn[:4]
vtk_types.append(self._VTK_TET)
elif etype == self._LINEAR_HEX:
ids = conn[:8]
vtk_types.append(self._VTK_HEX)
elif etype == self._UNSUPPORTED_ELEM:
n_skipped += 1
continue
else:
# No -1 values, append the whole cell
trimmed_connectivity.append(cell)
trimmed_connectivity = np.array(trimmed_connectivity, dtype="int32").flatten()
raise RuntimeError(
f"Invalid element type {etype} found in mesh {self.id}"
)
conn_list.extend(ids.tolist())
cell_sizes.append(len(ids))
# MOAB meshes supports tet elements only so we know it has 4 points per cell
points_per_cell = 4
if n_skipped > 0:
warnings.warn(
f"{n_skipped} elements were not written because "
"they are not of type linear tet/hex"
)
# offsets are the indices of the first point of each cell in the array of points
offsets = np.arange(0, self.n_elements * points_per_cell + 1, points_per_cell)
connectivity = np.asarray(conn_list, dtype=np.int64)
# Offsets must be length (numCells + 1) with a leading 0 and
# cumulative end-indices thereafter, per VTK's layout
cell_sizes_arr = np.asarray(cell_sizes, dtype=np.int64)
offsets = np.zeros(cell_sizes_arr.size + 1, dtype=np.int64)
np.cumsum(cell_sizes_arr, out=offsets[1:])
for name, data in datasets.items():
if data.shape != self.dimension:
@ -3156,42 +3236,27 @@ class UnstructuredMesh(MeshBase):
dtype=h5py.string_dtype("ascii", len(ascii_type)),
)
# create hdf5 file structure
root.create_dataset("NumberOfPoints", (0,), maxshape=(None,), dtype="i8")
root.create_dataset("Types", (0,), maxshape=(None,), dtype="uint8")
root.create_dataset("Points", (0, 3), maxshape=(None, 3), dtype="f")
root.create_dataset(
"NumberOfConnectivityIds", (0,), maxshape=(None,), dtype="i8"
)
root.create_dataset("NumberOfCells", (0,), maxshape=(None,), dtype="i8")
root.create_dataset("Offsets", (0,), maxshape=(None,), dtype="i8")
root.create_dataset("Connectivity", (0,), maxshape=(None,), dtype="i8")
# Create HDF5 file structure compliant with VTKHDF UnstructuredGrid
n_points = int(len(self.vertices))
n_cells = int(len(cell_sizes))
n_conn_ids = int(len(connectivity))
append_dataset(root["NumberOfPoints"], np.array([len(self.vertices)]))
append_dataset(root["Points"], self.vertices)
append_dataset(
root["NumberOfConnectivityIds"],
np.array([len(trimmed_connectivity)]),
)
append_dataset(root["Connectivity"], trimmed_connectivity)
append_dataset(root["NumberOfCells"], np.array([self.n_elements]))
append_dataset(root["Offsets"], offsets)
append_dataset(
root["Types"], np.full(self.n_elements, self._VTK_TETRA, dtype="uint8")
)
root.create_dataset("NumberOfPoints", data=(n_points,), dtype="i8")
root.create_dataset("NumberOfCells", data=(n_cells,), dtype="i8")
root.create_dataset("NumberOfConnectivityIds", data=(n_conn_ids,), dtype="i8")
root.create_dataset("Points", data=self.vertices.astype(np.float64, copy=False), dtype="f8")
root.create_dataset("Types", data=np.asarray(vtk_types, dtype=np.uint8), dtype="uint8")
root.create_dataset("Offsets", data=offsets.astype("i8"), dtype="i8")
root.create_dataset("Connectivity", data=connectivity.astype("i8"), dtype="i8")
cell_data_group = root.create_group("CellData")
for name, data in datasets.items():
cell_data_group.create_dataset(
name, (0,), maxshape=(None,), dtype="float64", chunks=True
)
if volume_normalization:
data /= self.volumes
append_dataset(cell_data_group[name], data)
cell_data_group.create_dataset(
name, data=data, dtype="float64", chunks=True
)
@classmethod
def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):

8
openmc/mgxs/library.py
View file

@ -1290,6 +1290,14 @@ class Library:
'are ignored since multiplicity or nu-scatter matrices '\
'were not tallied for ' + xsdata_name
warn(msg, RuntimeWarning)
if 'scatter matrix' in self.mgxs_types:
scatt_mgxs = self.get_mgxs(domain, 'scatter matrix')
elif 'consistent scatter matrix' in self.mgxs_types:
scatt_mgxs = self.get_mgxs(domain, 'consistent scatter matrix')
else:
raise ValueError(f'No scatter matrix found for {xsdata_name}.')
xsdata.set_scatter_matrix_mgxs(scatt_mgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],

201
openmc/model/model.py
View file

@ -265,22 +265,22 @@ class Model:
denom_tally = openmc.Tally(name='IFP denominator')
denom_tally.scores = ['ifp-denominator']
self.tallies.append(denom_tally)
# TODO: This should also be incorporated into lower-level calls in
# TODO: This should also be incorporated into lower-level calls in
# settings.py, but it requires information about the tallies currently
# on the active Model
def _assign_fw_cadis_tally_IDs(self):
# Verify that all tallies assigned as targets on WeightWindowGenerators
# exist within model.tallies. If this is the case, convert the .targets
# Verify that all tallies assigned as targets on WeightWindowGenerators
# exist within model.tallies. If this is the case, convert the .targets
# attribute of each WeightWindowGenerator to a sequence of tally IDs.
if len(self.settings.weight_window_generators) == 0:
return
# List of valid tally IDs
reference_tally_ids = np.asarray([tal.id for tal in self.tallies])
for wwg in self.settings.weight_window_generators:
# Only proceeds if the "targets" attribute is an openmc.Tallies,
# Only proceeds if the "targets" attribute is an openmc.Tallies,
# which means it hasn't been checked against model.tallies.
if isinstance(wwg.targets, openmc.Tallies):
id_vec = []
@ -291,8 +291,8 @@ class Model:
if tal == reference_tal:
id_next = reference_tal.id
break
if id_next == None:
if id_next is None:
raise RuntimeError(
f'Local FW-CADIS target tally {tal.id} not found on model.tallies!')
else:
@ -467,6 +467,8 @@ class Model:
This method iterates over all DAGMC universes in the geometry and
synchronizes their cells with the current material assignments. Requires
that the model has been initialized via :meth:`Model.init_lib`.
Synchronized DAGMC cells can then be edited and exported as nested
`<cell>` overrides inside each `<dagmc_universe>` element.
.. versionadded:: 0.15.1
@ -1125,28 +1127,148 @@ class Model:
array contains cell IDs, cell instances, and material IDs (in that
order).
"""
ids, _ = self.slice_data(
origin=origin,
width=width,
pixels=pixels,
basis=basis,
show_overlaps=color_overlaps,
level=-1,
include_properties=False,
**init_kwargs,
)
return ids
def slice_data(
self,
origin: Sequence[float] | None = None,
width: Sequence[float] | None = None,
pixels: int | Sequence[int] = 40000,
basis: str = 'xy',
u_span: Sequence[float] | None = None,
v_span: Sequence[float] | None = None,
show_overlaps: bool = False,
level: int = -1,
filter: openmc.Filter | None = None,
include_properties: bool = True,
**init_kwargs
) -> tuple[np.ndarray, np.ndarray | None]:
"""Generate geometry and property data for a 2D plot slice.
This method combines the functionality of :meth:`id_map` and property
mapping into a single call, avoiding duplicate geometry lookups. It also
supports filter bin index lookup for tally visualization.
.. versionadded:: 0.16.0
Parameters
----------
origin : Sequence[float], optional
Origin of the plot. If unspecified, this argument defaults to the
center of the bounding box if the bounding box does not contain inf
values for the provided basis, otherwise (0.0, 0.0, 0.0).
width : Sequence[float], optional
Width of the plot. If unspecified, this argument defaults to the
width of the bounding box if the bounding box does not contain inf
values for the provided basis, otherwise (10.0, 10.0).
pixels : int | Sequence[int], optional
If an iterable of ints is provided then this directly sets the
number of pixels to use in each basis direction. If a single int is
provided then this sets the total number of pixels in the plot and
the number of pixels in each basis direction is calculated from this
total and the image aspect ratio based on the width argument.
basis : {'xy', 'yz', 'xz'}, optional
Basis of the plot.
u_span : Sequence[float], optional
Full-width span vector for an oriented slice (3 values). Mutually
exclusive with width.
v_span : Sequence[float], optional
Full-height span vector for an oriented slice (3 values). Mutually
exclusive with width.
show_overlaps : bool, optional
Whether to identify and assign unique IDs (-3) to overlapping
regions. If False, overlapping regions will be assigned the ID of
the lowest-numbered cell that occupies that region. Defaults to
False.
level : int, optional
Universe level to plot (-1 for deepest). Defaults to -1.
filter : openmc.Filter, optional
If provided, the information for each pixel also includes an index
in the filter corresponding to the pixel position.
include_properties : bool, optional
Whether to include temperature/density data. Defaults to True.
**init_kwargs
Keyword arguments passed to :meth:`Model.init_lib`.
Returns
-------
geom_data : numpy.ndarray
Shape (v_res, h_res, 3) or (v_res, h_res, 4) int32 array. Contains
[cell_id, cell_instance, material_id] when no filter, or [cell_id,
cell_instance, material_id, filter_bin] with filter.
property_data : numpy.ndarray or None
Shape (v_res, h_res, 2) float64 array with [temperature, density],
or None if include_properties=False.
"""
import openmc.lib
origin, width, pixels = self._set_plot_defaults(
origin, width, pixels, basis)
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.")
# initialize the openmc.lib.plot._PlotBase object
plot_obj = openmc.lib.plot._PlotBase()
plot_obj.origin = origin
plot_obj.width = width[0]
plot_obj.height = width[1]
plot_obj.h_res = pixels[0]
plot_obj.v_res = pixels[1]
plot_obj.basis = basis
plot_obj.color_overlaps = color_overlaps
if u_span is not None or v_span is not None:
if u_span is None or v_span is None:
raise ValueError("Both u_span and v_span must be provided.")
if origin is None:
origin = (0.0, 0.0, 0.0)
if isinstance(pixels, int):
u_norm = np.linalg.norm(u_span)
v_norm = np.linalg.norm(v_span)
aspect_ratio = u_norm / v_norm
pixels_y = math.sqrt(pixels / aspect_ratio)
pixels = (int(pixels / pixels_y), int(pixels_y))
else:
origin, width, pixels = self._set_plot_defaults(
origin, width, pixels, basis)
# Silence output by default. Also set arguments to start in volume
# calculation mode to avoid loading cross sections
init_kwargs.setdefault('output', False)
init_kwargs.setdefault('args', ['-c'])
# If filter does not already appear in the model, temporarily add a
# tally with the filter
original_length = len(self.tallies)
if filter is not None:
filter_ids = {f.id for t in self.tallies for f in t.filters}
if filter.id not in filter_ids:
# Create temporary tally while preserving ID assignment
next_id = openmc.Tally.next_id
temp_tally = openmc.Tally()
temp_tally.filters = [filter]
temp_tally.scores = ['flux']
self.tallies.append(temp_tally)
openmc.Tally.used_ids.remove(temp_tally.id)
openmc.Tally.next_id = next_id
with openmc.lib.TemporarySession(self, **init_kwargs):
return openmc.lib.id_map(plot_obj)
geom_data, property_data = openmc.lib.slice_data(
origin=origin,
width=width,
basis=basis,
u_span=u_span,
v_span=v_span,
pixels=pixels,
show_overlaps=show_overlaps,
level=level,
filter=filter,
include_properties=include_properties,
)
# If filter was temporarily added, remove it
if len(self.tallies) > original_length:
self.tallies.pop()
return geom_data, property_data
@add_plot_params
def plot(
@ -1214,13 +1336,14 @@ class Model:
"openmc.config before plotting.")
break
# Get ID map from the C API
id_map = self.id_map(
# Get plot IDs from the C API
id_map, _ = self.slice_data(
origin=origin,
width=width,
pixels=pixels,
basis=basis,
color_overlaps=show_overlaps
show_overlaps=show_overlaps,
include_properties=False,
)
# Generate colors if not provided
@ -1746,10 +1869,10 @@ class Model:
@staticmethod
def _auto_generate_mgxs_lib(
model: openmc.model.model,
model: openmc.model.Model,
groups: openmc.mgxs.EnergyGroups,
correction: str | none,
directory: pathlike,
correction: str | None,
directory: PathLike,
) -> openmc.mgxs.Library:
"""
Automatically generate a multi-group cross section libray from a model
@ -1956,7 +2079,7 @@ class Model:
# Set materials on the model
model.materials = [material]
if temperature != None:
if temperature is not None:
model.materials[-1].temperature = temperature
# Settings
@ -1983,7 +2106,7 @@ class Model:
mgxs_lib = Model._auto_generate_mgxs_lib(
model, groups, correction, directory)
if temperature != None:
if temperature is not None:
return mgxs_lib.get_xsdata(domain=material, xsdata_name=name,
temperature=temperature)
else:
@ -2056,12 +2179,12 @@ class Model:
)
temp_settings = {}
if temperature_settings == None:
if temperature_settings is None:
temp_settings = self.settings.temperature
else:
temp_settings = temperature_settings
if temperatures == None:
if temperatures is None:
mgxs_sets = []
for material in self.materials:
xs_data = Model._isothermal_infinite_media_mgxs(
@ -2234,7 +2357,7 @@ class Model:
model = openmc.Model()
model.geometry = stoch_geom
if temperature != None:
if temperature is not None:
for material in model.geometry.get_all_materials().values():
material.temperature = temperature
@ -2258,7 +2381,7 @@ class Model:
model, groups, correction, directory)
# Fetch all of the isothermal results.
if temperature != None:
if temperature is not None:
return {
mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name,
temperature=temperature)
@ -2344,12 +2467,12 @@ class Model:
)
temp_settings = {}
if temperature_settings == None:
if temperature_settings is None:
temp_settings = self.settings.temperature
else:
temp_settings = temperature_settings
if temperatures == None:
if temperatures is None:
mgxs_sets = Model._isothermal_stochastic_slab_mgxs(
geo,
groups,
@ -2442,7 +2565,7 @@ class Model:
model = copy.deepcopy(input_model)
model.tallies = openmc.Tallies()
if temperature != None:
if temperature is not None:
for material in model.geometry.get_all_materials().values():
material.temperature = temperature
@ -2458,7 +2581,7 @@ class Model:
model, groups, correction, directory)
# Fetch all of the isothermal results.
if temperature != None:
if temperature is not None:
return {
mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name,
temperature=temperature)
@ -2513,12 +2636,12 @@ class Model:
entries in openmc.Settings.temperature_settings.
"""
temp_settings = {}
if temperature_settings == None:
if temperature_settings is None:
temp_settings = self.settings.temperature
else:
temp_settings = temperature_settings
if temperatures == None:
if temperatures is None:
mgxs_sets = Model._isothermal_materialwise_mgxs(
self,
groups,

43
openmc/settings.py
View file

@ -67,6 +67,10 @@ class Settings:
(ex: ["(n,fission)", 2, "(n,2n)"] ). (list of str or int)
:deposited_E_threshold: Number to define the minimum deposited energy during
per collision to trigger banking. (float)
create_delayed_neutrons : bool
Whether delayed neutrons are created in fission.
.. versionadded:: 0.13.3
create_fission_neutrons : bool
Indicate whether fission neutrons should be created or not.
cutoff : dict
@ -256,8 +260,15 @@ class Settings:
The type of calculation to perform (default is 'eigenvalue')
seed : int
Seed for the linear congruential pseudorandom number generator
stride : int
Number of random numbers allocated for each source particle history
shared_secondary_bank : bool
Whether to use a shared secondary particle bank. When enabled,
secondary particles are collected into a global bank, sorted for
reproducibility, and load-balanced across MPI ranks between
generations. If not specified, the shared secondary bank is
enabled automatically for fixed-source simulations with weight
windows active, and disabled otherwise.
.. versionadded:: 0.15.4
source : Iterable of openmc.SourceBase
Distribution of source sites in space, angle, and energy
source_rejection_fraction : float
@ -277,6 +288,8 @@ class Settings:
Options for writing state points. Acceptable keys are:
:batches: list of batches at which to write statepoint files
stride : int
Number of random numbers allocated for each source particle history
surf_source_read : dict
Options for reading surface source points. Acceptable keys are:
@ -373,10 +386,6 @@ class Settings:
.. versionadded:: 0.14.0
create_delayed_neutrons : bool
Whether delayed neutrons are created in fission.
.. versionadded:: 0.13.3
weight_windows_on : bool
Whether weight windows are enabled
@ -482,6 +491,7 @@ class Settings:
self._weight_window_generators = cv.CheckedList(
WeightWindowGenerator, 'weight window generators')
self._weight_windows_on = None
self._shared_secondary_bank = None
self._weight_windows_file = None
self._weight_window_checkpoints = {}
self._max_history_splits = None
@ -1306,6 +1316,15 @@ class Settings:
cv.check_type('weight windows on', value, bool)
self._weight_windows_on = value
@property
def shared_secondary_bank(self) -> bool:
return self._shared_secondary_bank
@shared_secondary_bank.setter
def shared_secondary_bank(self, value: bool):
cv.check_type('shared secondary bank', value, bool)
self._shared_secondary_bank = value
@property
def weight_window_checkpoints(self) -> dict:
return self._weight_window_checkpoints
@ -1924,6 +1943,11 @@ class Settings:
elem = ET.SubElement(root, "weight_windows_on")
elem.text = str(self._weight_windows_on).lower()
def _create_shared_secondary_bank_subelement(self, root):
if self._shared_secondary_bank is not None:
elem = ET.SubElement(root, "shared_secondary_bank")
elem.text = str(self._shared_secondary_bank).lower()
def _create_weight_window_generators_subelement(self, root, mesh_memo=None):
if not self.weight_window_generators:
return
@ -2430,6 +2454,11 @@ class Settings:
if text is not None:
self.weight_windows_on = text in ('true', '1')
def _shared_secondary_bank_from_xml_element(self, root):
text = get_text(root, 'shared_secondary_bank')
if text is not None:
self.shared_secondary_bank = text in ('true', '1')
def _weight_windows_file_from_xml_element(self, root):
text = get_text(root, 'weight_windows_file')
if text is not None:
@ -2597,6 +2626,7 @@ class Settings:
self._create_write_initial_source_subelement(element)
self._create_weight_windows_subelement(element, mesh_memo)
self._create_weight_windows_on_subelement(element)
self._create_shared_secondary_bank_subelement(element)
self._create_weight_window_generators_subelement(element, mesh_memo)
self._create_weight_windows_file_element(element)
self._create_weight_window_checkpoints_subelement(element)
@ -2714,6 +2744,7 @@ class Settings:
settings._write_initial_source_from_xml_element(elem)
settings._weight_windows_from_xml_element(elem, meshes)
settings._weight_windows_on_from_xml_element(elem)
settings._shared_secondary_bank_from_xml_element(elem)
settings._weight_windows_file_from_xml_element(elem)
settings._weight_window_generators_from_xml_element(elem, meshes)
settings._weight_window_checkpoints_from_xml_element(elem)

309
openmc/stats/univariate.py
View file

@ -2,8 +2,10 @@ from __future__ import annotations
from abc import ABC, abstractmethod
from collections import defaultdict
from collections.abc import Iterable, Sequence
from functools import cache
from math import sqrt, pi, exp, log
from numbers import Real
from pathlib import Path
from warnings import warn
import lxml.etree as ET
@ -14,6 +16,7 @@ import scipy
import openmc.checkvalue as cv
from openmc.data import atomic_mass, NEUTRON_MASS
import openmc.data
from .._xml import get_elem_list, get_text
from ..mixin import EqualityMixin
@ -124,6 +127,8 @@ class Univariate(EqualityMixin, ABC):
return Legendre.from_xml_element(elem)
elif distribution == 'mixture':
return Mixture.from_xml_element(elem)
elif distribution == 'decay_spectrum':
return DecaySpectrum.from_xml_element(elem)
@abstractmethod
def _sample_unbiased(self, n_samples: int = 1, seed: int | None = None):
@ -2196,6 +2201,310 @@ class Mixture(Univariate):
return new_dist
class DecaySpectrum(Univariate):
"""Energy distribution from decay photon spectra of a mixture of nuclides.
This distribution stores nuclide names, their atom densities, and the volume
of the region. When written to XML and read by the C++ solver, the nuclide
names are resolved against the depletion chain to obtain the decay photon
energy spectra and decay constants. The resulting distribution is a mixture
of per-nuclide photon spectra weighted by absolute activity. The volume is
necessary so that the C++ solver can compute the total photon emission rate
in [photons/s], which is used as the source strength.
.. versionadded:: 0.15.4
Parameters
----------
nuclides : dict
Dictionary mapping nuclide name (str) to atom density (float) in units
of [atom/b-cm].
volume : float
Volume of the source region in [cm³]. Used together with atom densities
to compute the absolute photon emission rate.
Attributes
----------
nuclides : dict
Dictionary mapping nuclide name to atom density in [atom/b-cm].
volume : float
Volume of the source region in [cm³].
"""
def __init__(self, nuclides: dict[str, float], volume: float):
super().__init__(bias=None)
self._dist_cache = None
self._dist_cache_key = None
self.nuclides = nuclides
self.volume = volume
def __len__(self):
return len(self.nuclides)
@property
def nuclides(self):
return self._nuclides
@nuclides.setter
def nuclides(self, nuclides):
cv.check_type('nuclides', nuclides, dict)
for name, density in nuclides.items():
cv.check_type('nuclide name', name, str)
cv.check_type(f'atom density for {name}', density, Real)
cv.check_greater_than(f'atom density for {name}', density, 0.0, True)
self._nuclides = dict(nuclides)
self._dist_cache = None
self._dist_cache_key = None
@property
def volume(self):
return self._volume
@volume.setter
def volume(self, volume):
cv.check_type('volume', volume, Real)
cv.check_greater_than('volume', volume, 0.0)
self._volume = float(volume)
self._dist_cache = None
self._dist_cache_key = None
@staticmethod
def _chain_file_cache_key():
"""Return a hashable key for the active depletion chain."""
chain_file = openmc.config.get('chain_file')
if chain_file is None:
return None
path = Path(chain_file).resolve()
try:
stat = path.stat()
except OSError:
return (path, None, None)
return (path, stat.st_mtime, stat.st_size)
def to_distribution(self):
"""Convert to a concrete distribution using decay chain data.
Builds a combined photon energy distribution by looking up each nuclide
in the depletion chain via :func:`openmc.data.decay_photon_energy` and
weighting by absolute atom count (``density * 1e24 * volume``). The
result is cached on the object; the cache is invalidated automatically
when :attr:`nuclides` or :attr:`volume` are reassigned.
Requires ``openmc.config['chain_file']`` to be set.
Returns
-------
openmc.stats.Univariate or None
Combined photon energy distribution, or ``None`` if no nuclide in
:attr:`nuclides` has a photon source in the chain.
"""
chain_key = self._chain_file_cache_key()
if self._dist_cache is not None and self._dist_cache_key == chain_key:
return self._dist_cache
dists = []
weights = []
for name, density in self.nuclides.items():
dist = openmc.data.decay_photon_energy(name)
if dist is not None:
dists.append(dist)
weights.append(density * 1e24 * self.volume)
if not dists:
return None
self._dist_cache = combine_distributions(dists, weights)
self._dist_cache_key = chain_key
return self._dist_cache
def to_xml_element(self, element_name: str):
"""Return XML representation of the decay photon distribution
Parameters
----------
element_name : str
XML element name
Returns
-------
element : lxml.etree._Element
XML element containing decay photon distribution data
"""
element = ET.Element(element_name)
element.set("type", "decay_spectrum")
element.set("volume", str(self.volume))
nuclides = ET.SubElement(element, "nuclides")
nuclides.text = ' '.join(self.nuclides)
parameters = ET.SubElement(element, "parameters")
parameters.text = ' '.join(str(density) for density in self.nuclides.values())
return element
@classmethod
def from_xml_element(cls, elem: ET.Element):
"""Generate decay photon distribution from an XML element
Parameters
----------
elem : lxml.etree._Element
XML element
Returns
-------
openmc.stats.DecaySpectrum
Decay photon distribution generated from XML element
"""
volume = float(elem.get('volume'))
names = get_elem_list(elem, 'nuclides', str)
densities = get_elem_list(elem, 'parameters', float)
nuclides = dict(zip(names, densities))
return cls(nuclides, volume)
def _sample_unbiased(self, n_samples=1, seed=None):
dist = self.to_distribution()
if dist is None:
raise RuntimeError(
"DecaySpectrum._sample_unbiased requires chain data but none "
"was found. Ensure openmc.config['chain_file'] is set and the "
"chain contains photon sources for the nuclides present."
)
return dist.sample(n_samples, seed)[0]
def integral(self):
"""Return integral of the distribution
Returns the total photon emission rate in [photons/s] by delegating to
:meth:`to_distribution`. Returns ``0.0`` when no chain data is
available (e.g., ``openmc.config['chain_file']`` is not set).
Returns
-------
float
Total photon emission rate in [photons/s], or ``0.0`` if chain
data is unavailable.
"""
try:
dist = self.to_distribution()
except Exception:
return 0.0
if dist is None:
return 0.0
return dist.integral()
@staticmethod
@cache
def _photon_integral(nuclide: str, chain_key) -> float | None:
"""Return the per-atom photon emission integral for a nuclide"""
dist = openmc.data.decay_photon_energy(nuclide)
return dist.integral() if dist is not None else None
def clip(self, tolerance: float = 1e-9, inplace: bool = False):
"""Remove nuclides with negligible contribution to photon emission.
Nuclides that are stable or have no photon source in the depletion
chain are removed unconditionally. The remaining nuclides are ranked
by their photon emission rate (proportional to
``atom_density * decay_constant * photon_yield``) and the least
important are discarded until the cumulative discarded fraction of the
total emission rate exceeds *tolerance*.
Requires ``openmc.config['chain_file']`` to be set.
Parameters
----------
tolerance : float
Maximum fraction of total photon emission rate that may be
discarded.
inplace : bool
Whether to modify the current object in-place or return a new one.
Returns
-------
openmc.stats.DecaySpectrum
Distribution with negligible nuclides removed.
"""
# Compute per-nuclide emission rate; drop non-emitters
emitting_names = []
emitting_densities = []
rates = []
chain_key = self._chain_file_cache_key()
for name, density in self.nuclides.items():
integral = DecaySpectrum._photon_integral(name, chain_key)
if integral is None:
continue
emitting_names.append(name)
emitting_densities.append(density)
rates.append(density * self.volume * integral)
if not emitting_names:
new_nuclides = {}
else:
indices = _intensity_clip(rates, tolerance=tolerance)
new_nuclides = {
emitting_names[i]: emitting_densities[i] for i in indices
}
if inplace:
self._nuclides = new_nuclides
self._dist_cache = None
self._dist_cache_key = None
return self
return type(self)(new_nuclides, self.volume)
@property
def support(self):
return (0.0, np.inf)
def evaluate(self, x):
"""Evaluate the probability density at a given value.
Delegates to the combined distribution built from chain data. Raises
``NotImplementedError`` if the combined distribution is a
:class:`~openmc.stats.Mixture` (which does not support
``evaluate()``).
Parameters
----------
x : float
Value at which to evaluate the PDF.
Returns
-------
float
Probability density at *x*.
"""
dist = self.to_distribution()
if dist is None:
raise RuntimeError(
"DecaySpectrum.evaluate requires chain data. Ensure "
"openmc.config['chain_file'] is set."
)
return dist.evaluate(x)
def mean(self):
"""Return the mean of the distribution.
Delegates to the combined distribution built from chain data.
Returns
-------
float
Mean photon energy in [eV].
"""
dist = self.to_distribution()
if dist is None:
raise RuntimeError(
"DecaySpectrum.mean requires chain data. Ensure "
"openmc.config['chain_file'] is set."
)
return dist.mean()
def combine_distributions(
dists: Sequence[Discrete | Tabular | Mixture],
probs: Sequence[float]

2
pyproject.toml
View file

@ -69,7 +69,7 @@ include = ['openmc*']
exclude = ['tests*']
[tool.setuptools.package-data]
"openmc.data.dose" = ["**/*.txt"]
"openmc.data.dose" = ["**/*.txt", "*.h5"]
"openmc.data" = ["*.txt", "*.DAT", "*.json", "*.h5"]
"openmc.lib" = ["libopenmc.dylib", "libopenmc.so"]

202
src/bank.cpp
View file

@ -7,6 +7,7 @@
#include "openmc/vector.h"
#include <cstdint>
#include <limits>
#include <numeric>
namespace openmc {
@ -43,6 +44,13 @@ vector<vector<double>> ifp_fission_lifetime_bank;
// used to efficiently sort the fission bank after each iteration.
vector<int64_t> progeny_per_particle;
// When shared secondary bank mode is enabled, secondaries produced during
// transport are collected in the write bank. When a secondary generation is
// complete, write is moved to read for transport, and a new empty write bank
// is created. This repeats until no secondaries remain.
SharedArray<SourceSite> shared_secondary_bank_read;
SharedArray<SourceSite> shared_secondary_bank_write;
} // namespace simulation
//==============================================================================
@ -60,6 +68,8 @@ void free_memory_bank()
simulation::ifp_source_lifetime_bank.clear();
simulation::ifp_fission_delayed_group_bank.clear();
simulation::ifp_fission_lifetime_bank.clear();
simulation::shared_secondary_bank_read.clear();
simulation::shared_secondary_bank_write.clear();
}
void init_fission_bank(int64_t max)
@ -68,13 +78,13 @@ void init_fission_bank(int64_t max)
simulation::progeny_per_particle.resize(simulation::work_per_rank);
}
// Performs an O(n) sort on the fission bank, by leveraging
// Performs an O(n) sort on a fission or secondary bank, by leveraging
// the parent_id and progeny_id fields of banked particles. See the following
// paper for more details:
// "Reproducibility and Monte Carlo Eigenvalue Calculations," F.B. Brown and
// T.M. Sutton, 1992 ANS Annual Meeting, Transactions of the American Nuclear
// Society, Volume 65, Page 235.
void sort_fission_bank()
void sort_bank(SharedArray<SourceSite>& bank, bool is_fission_bank)
{
// Ensure we don't read off the end of the array if we ran with 0 particles
if (simulation::progeny_per_particle.size() == 0) {
@ -96,44 +106,200 @@ void sort_fission_bank()
vector<vector<double>> sorted_ifp_lifetime_bank;
// If there is not enough space, allocate a temporary vector and point to it
if (simulation::fission_bank.size() >
simulation::fission_bank.capacity() / 2) {
sorted_bank_holder.resize(simulation::fission_bank.size());
if (bank.size() > bank.capacity() / 2) {
sorted_bank_holder.resize(bank.size());
sorted_bank = sorted_bank_holder.data();
} else { // otherwise, point sorted_bank to unused portion of the fission bank
sorted_bank = &simulation::fission_bank[simulation::fission_bank.size()];
sorted_bank = bank.data() + bank.size();
}
if (settings::ifp_on) {
if (settings::ifp_on && is_fission_bank) {
allocate_temporary_vector_ifp(
sorted_ifp_delayed_group_bank, sorted_ifp_lifetime_bank);
}
// Use parent and progeny indices to sort fission bank
for (int64_t i = 0; i < simulation::fission_bank.size(); i++) {
const auto& site = simulation::fission_bank[i];
int64_t offset = site.parent_id - 1 - simulation::work_index[mpi::rank];
int64_t idx = simulation::progeny_per_particle[offset] + site.progeny_id;
if (idx >= simulation::fission_bank.size()) {
// Use parent and progeny indices to sort bank
for (int64_t i = 0; i < bank.size(); i++) {
const auto& site = bank[i];
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 >= bank.size()) {
fatal_error("Mismatch detected between sum of all particle progeny and "
"shared fission bank size.");
"bank size during sorting.");
}
sorted_bank[idx] = site;
if (settings::ifp_on) {
if (settings::ifp_on && is_fission_bank) {
copy_ifp_data_from_fission_banks(
i, sorted_ifp_delayed_group_bank[idx], sorted_ifp_lifetime_bank[idx]);
}
}
// Copy sorted bank into the fission bank
std::copy(sorted_bank, sorted_bank + simulation::fission_bank.size(),
simulation::fission_bank.data());
if (settings::ifp_on) {
std::copy(sorted_bank, sorted_bank + bank.size(), bank.data());
if (settings::ifp_on && is_fission_bank) {
copy_ifp_data_to_fission_banks(
sorted_ifp_delayed_group_bank.data(), sorted_ifp_lifetime_bank.data());
}
}
// This function redistributes SourceSite particles across MPI ranks to
// achieve load balancing while preserving the global ordering of particles.
//
// GUARANTEES:
// -----------
// 1. Global Order Preservation: After redistribution, each rank holds a
// contiguous slice of the original global ordering. For example, if the
// input across 3 ranks was:
// - Rank 0: IDs 0-4
// - Rank 1: IDs 5-6
// - Rank 2: IDs 7-200
// Then after redistribution (assuming ~67 particles per rank):
// - Rank 0: IDs 0-66 (contiguous)
// - Rank 1: IDs 67-133 (contiguous)
// - Rank 2: IDs 134-200 (contiguous)
// The global ordering is always preserved - no rank will ever hold
// non-contiguous ID ranges like "0-4 and 100-200".
//
// 2. Even Load Balancing: Particles are distributed as evenly as possible.
// If total % n_procs != 0, the first 'remainder' ranks each get one extra
// particle (i.e., floor division with remainder distributed to lower
// ranks). This follows the same logic as calculate_work().
//
// HOW IT WORKS:
// -------------
// The algorithm uses overlap-based redistribution:
// 1. Each rank's current data occupies a range [cumulative_before[rank],
// cumulative_before[rank+1]) in the global index space.
// 2. Each rank's target data should occupy [cumulative_target[rank],
// cumulative_target[rank+1]) in the same global index space.
// 3. For each pair of (source_rank, dest_rank), we calculate the overlap
// between what source_rank currently has and what dest_rank needs.
// 4. MPI_Alltoallv transfers exactly these overlapping regions, with
// displacements ensuring data lands at the correct position in the
// receiving buffer.
//
// EDGE CASES HANDLED:
// -------------------
// - Single rank (n_procs == 1): Returns immediately with local size, no MPI.
// - Empty total (all ranks have 0 particles): Returns 0 immediately.
// - Imbalanced input (e.g., one rank has all particles): Works correctly;
// that rank will send portions to all other ranks based on target ranges.
// - Non-divisible totals: First 'remainder' ranks get one extra particle.
int64_t synchronize_global_secondary_bank(
SharedArray<SourceSite>& shared_secondary_bank)
{
// Get current size of local bank
int64_t local_size = shared_secondary_bank.size();
if (mpi::n_procs == 1) {
return local_size;
}
#ifdef OPENMC_MPI
// Gather all sizes to all ranks
vector<int64_t> all_sizes(mpi::n_procs);
MPI_Allgather(&local_size, 1, MPI_INT64_T, all_sizes.data(), 1, MPI_INT64_T,
mpi::intracomm);
// Calculate total and check for empty case
int64_t total = 0;
for (int64_t size : all_sizes) {
total += size;
}
// If we don't have any items to distribute, return
if (total == 0) {
return total;
}
int64_t base_count = total / mpi::n_procs;
int64_t remainder = total % mpi::n_procs;
// Calculate target size for each rank
// First 'remainder' ranks get base_count + 1, rest get base_count
vector<int64_t> target_sizes(mpi::n_procs);
for (int i = 0; i < mpi::n_procs; ++i) {
target_sizes[i] = base_count + (i < remainder ? 1 : 0);
}
// Calculate send and receive counts in terms of SourceSite objects
// (not bytes)
vector<int> send_counts(mpi::n_procs, 0);
vector<int> recv_counts(mpi::n_procs, 0);
vector<int> send_displs(mpi::n_procs, 0);
vector<int> recv_displs(mpi::n_procs, 0);
// Calculate cumulative positions (starting index for each rank in the
// global array)
vector<int64_t> cumulative_before(mpi::n_procs + 1, 0);
vector<int64_t> cumulative_target(mpi::n_procs + 1, 0);
for (int i = 0; i < mpi::n_procs; ++i) {
cumulative_before[i + 1] = cumulative_before[i] + all_sizes[i];
cumulative_target[i + 1] = cumulative_target[i] + target_sizes[i];
}
// Determine send and receive amounts for each rank
int64_t my_start = cumulative_before[mpi::rank];
int64_t my_end = cumulative_before[mpi::rank + 1];
int64_t my_target_start = cumulative_target[mpi::rank];
int64_t my_target_end = cumulative_target[mpi::rank + 1];
for (int r = 0; r < mpi::n_procs; ++r) {
// Send: overlap between my current range and rank r's target range
int64_t send_overlap_start = std::max(my_start, cumulative_target[r]);
int64_t send_overlap_end = std::min(my_end, cumulative_target[r + 1]);
if (send_overlap_start < send_overlap_end) {
int64_t count = send_overlap_end - send_overlap_start;
int64_t displ = send_overlap_start - my_start;
if (count > std::numeric_limits<int>::max() ||
displ > std::numeric_limits<int>::max()) {
fatal_error("Secondary bank size exceeds MPI_Alltoallv int limit.");
}
send_counts[r] = static_cast<int>(count);
send_displs[r] = static_cast<int>(displ);
}
// Recv: overlap between rank r's current range and my target range
int64_t recv_overlap_start =
std::max(cumulative_before[r], my_target_start);
int64_t recv_overlap_end =
std::min(cumulative_before[r + 1], my_target_end);
if (recv_overlap_start < recv_overlap_end) {
int64_t count = recv_overlap_end - recv_overlap_start;
int64_t displ = recv_overlap_start - my_target_start;
if (count > std::numeric_limits<int>::max() ||
displ > std::numeric_limits<int>::max()) {
fatal_error("Secondary bank size exceeds MPI_Alltoallv int limit.");
}
recv_counts[r] = static_cast<int>(count);
recv_displs[r] = static_cast<int>(displ);
}
}
// Prepare receive buffer with target size
SharedArray<SourceSite> new_bank(target_sizes[mpi::rank]);
// Perform all-to-all redistribution using the custom MPI type
MPI_Alltoallv(shared_secondary_bank.data(), send_counts.data(),
send_displs.data(), mpi::source_site, new_bank.data(), recv_counts.data(),
recv_displs.data(), mpi::source_site, mpi::intracomm);
// Replace old bank with redistributed data
shared_secondary_bank = std::move(new_bank);
return total;
#else
return local_size;
#endif
}
//==============================================================================
// C API
//==============================================================================

94
src/cell.cpp
View file

@ -340,6 +340,63 @@ void Cell::to_hdf5(hid_t cell_group) const
close_group(group);
}
//==============================================================================
// XML parsing helpers for <cell> nodes
//==============================================================================
vector<int32_t> parse_cell_material_xml(pugi::xml_node node, int32_t cell_id)
{
vector<std::string> mats {
get_node_array<std::string>(node, "material", true)};
if (mats.empty()) {
fatal_error(fmt::format(
"An empty material element was specified for cell {}", cell_id));
}
vector<int32_t> material;
material.reserve(mats.size());
for (const auto& mat : mats) {
if (mat == "void") {
material.push_back(MATERIAL_VOID);
} else {
material.push_back(std::stoi(mat));
}
}
return material;
}
vector<double> parse_cell_temperature_xml(pugi::xml_node node, int32_t cell_id)
{
auto temperatures = get_node_array<double>(node, "temperature");
if (temperatures.empty()) {
fatal_error(fmt::format(
"An empty temperature element was specified for cell {}", cell_id));
}
for (auto T : temperatures) {
if (T < 0) {
fatal_error(fmt::format(
"Cell {} was specified with a negative temperature", cell_id));
}
}
return temperatures;
}
vector<double> parse_cell_density_xml(pugi::xml_node node, int32_t cell_id)
{
auto densities = get_node_array<double>(node, "density");
if (densities.empty()) {
fatal_error(fmt::format(
"An empty density element was specified for cell {}", cell_id));
}
for (auto rho : densities) {
if (rho <= 0) {
fatal_error(fmt::format(
"Cell {} was specified with a density less than or equal to zero",
cell_id));
}
}
return densities;
}
//==============================================================================
// CSGCell implementation
//==============================================================================
@ -390,26 +447,12 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
// universe), more than one material (distribmats), and some materials may
// be "void".
if (material_present) {
vector<std::string> mats {
get_node_array<std::string>(cell_node, "material", true)};
if (mats.size() > 0) {
material_.reserve(mats.size());
for (std::string mat : mats) {
if (mat.compare("void") == 0) {
material_.push_back(MATERIAL_VOID);
} else {
material_.push_back(std::stoi(mat));
}
}
} else {
fatal_error(fmt::format(
"An empty material element was specified for cell {}", id_));
}
material_ = parse_cell_material_xml(cell_node, id_);
}
// Read the temperature element which may be distributed like materials.
if (check_for_node(cell_node, "temperature")) {
sqrtkT_ = get_node_array<double>(cell_node, "temperature");
sqrtkT_ = parse_cell_temperature_xml(cell_node, id_);
sqrtkT_.shrink_to_fit();
// Make sure this is a material-filled cell.
@ -420,14 +463,6 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
id_));
}
// Make sure all temperatures are non-negative.
for (auto T : sqrtkT_) {
if (T < 0) {
fatal_error(fmt::format(
"Cell {} was specified with a negative temperature", id_));
}
}
// Convert to sqrt(k*T).
for (auto& T : sqrtkT_) {
T = std::sqrt(K_BOLTZMANN * T);
@ -440,7 +475,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
// Note: calculating the actual density multiplier is deferred until materials
// are finalized. density_mult_ contains the true density in the meantime.
if (check_for_node(cell_node, "density")) {
density_mult_ = get_node_array<double>(cell_node, "density");
density_mult_ = parse_cell_density_xml(cell_node, id_);
density_mult_.shrink_to_fit();
// Make sure this is a material-filled cell.
@ -461,15 +496,6 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
id_));
}
}
// Make sure all densities are non-negative and greater than zero.
for (auto rho : density_mult_) {
if (rho <= 0) {
fatal_error(fmt::format(
"Cell {} was specified with a density less than or equal to zero",
id_));
}
}
}
// Read the region specification.

8
src/chain.cpp
View file

@ -98,6 +98,8 @@ vector<unique_ptr<ChainNuclide>> chain_nuclides;
void read_chain_file_xml()
{
free_memory_chain();
char* chain_file_path = std::getenv("OPENMC_CHAIN_FILE");
if (!chain_file_path) {
return;
@ -120,4 +122,10 @@ void read_chain_file_xml()
}
}
void free_memory_chain()
{
data::chain_nuclides.clear();
data::chain_nuclide_map.clear();
}
} // namespace openmc

12
src/collision_track.cpp
View file

@ -200,8 +200,13 @@ void collision_track_record(Particle& particle)
return;
int cell_id = model::cells[cell_index]->id_;
const auto* nuclide_ptr = data::nuclides[particle.event_nuclide()].get();
std::string nuclide = nuclide_ptr->name_;
std::string nuclide {};
int nuclide_id = 0;
if (particle.event_nuclide() != NUCLIDE_NONE) {
const auto* nuclide_ptr = data::nuclides[particle.event_nuclide()].get();
nuclide = nuclide_ptr->name_;
nuclide_id = nuclide_ptr->particle_type().pdg_number();
}
int universe_id = model::universes[particle.lowest_coord().universe()]->id_;
double delta_E = particle.E_last() - particle.E();
int material_index = particle.material();
@ -224,8 +229,7 @@ void collision_track_record(Particle& particle)
site.event_mt = particle.event_mt();
site.delayed_group = particle.delayed_group();
site.cell_id = cell_id;
site.nuclide_id =
10000 * nuclide_ptr->Z_ + 10 * nuclide_ptr->A_ + nuclide_ptr->metastable_;
site.nuclide_id = nuclide_id;
site.material_id = material_id;
site.universe_id = universe_id;
site.n_collision = particle.n_collision();

227
src/dagmc.cpp
View file

@ -50,6 +50,10 @@ namespace openmc {
DAGUniverse::DAGUniverse(pugi::xml_node node)
{
MaterialOverrides material_overrides;
TemperatureOverrides temperature_overrides;
DensityOverrides density_overrides;
if (check_for_node(node, "id")) {
id_ = std::stoi(get_node_value(node, "id"));
} else {
@ -76,24 +80,77 @@ DAGUniverse::DAGUniverse(pugi::xml_node node)
adjust_material_ids_ = get_node_value_bool(node, "auto_mat_ids");
}
// get material assignment overloading
if (check_for_node(node, "material_overrides")) {
auto mat_node = node.child("material_overrides");
// loop over all subelements (each subelement corresponds to a material)
for (pugi::xml_node cell_node : mat_node.children("cell_override")) {
// Store assignment reference name
int32_t ref_assignment = std::stoi(get_node_value(cell_node, "id"));
// Get material assignment overrides from nested DAGMC cell elements.
if (node.child("cell")) {
for (pugi::xml_node cell_node : node.children("cell")) {
if (!check_for_node(cell_node, "id")) {
fatal_error(
"Must specify id for each DAGMC cell override in <dagmc_universe>.");
}
// Get mat name for each assignement instances
vector<int32_t> instance_mats =
get_node_array<int32_t>(cell_node, "material_ids");
int32_t cell_id = std::stoi(get_node_value(cell_node, "id"));
// Store mat name for each instances
material_overrides_.emplace(ref_assignment, instance_mats);
if (check_for_node(cell_node, "region")) {
fatal_error(fmt::format(
"DAGMC cell {} override cannot specify a region.", cell_id));
}
if (check_for_node(cell_node, "fill")) {
fatal_error(fmt::format(
"DAGMC cell {} override currently only supports material fills.",
cell_id));
}
if (check_for_node(cell_node, "universe")) {
fatal_error(fmt::format(
"DAGMC cell {} override cannot specify a universe.", cell_id));
}
if (check_for_node(cell_node, "translation") ||
check_for_node(cell_node, "rotation")) {
fatal_error(fmt::format(
"DAGMC cell {} override does not support translation or rotation.",
cell_id));
}
if (!check_for_node(cell_node, "material")) {
fatal_error(fmt::format(
"DAGMC cell {} override must specify material.", cell_id));
}
auto inserted = material_overrides.emplace(
cell_id, parse_cell_material_xml(cell_node, cell_id));
if (!inserted.second) {
fatal_error(fmt::format(
"Duplicate DAGMC cell override specified for cell {}", cell_id));
}
if (check_for_node(cell_node, "temperature")) {
temperature_overrides.emplace(
cell_id, parse_cell_temperature_xml(cell_node, cell_id));
}
if (check_for_node(cell_node, "density")) {
density_overrides.emplace(
cell_id, parse_cell_density_xml(cell_node, cell_id));
}
}
} else if (check_for_node(node, "material_overrides")) {
if (node.child("cell")) {
fatal_error("DAGMCUniverse cannot specify both <material_overrides> and "
"<cell> sub-elements. Use <cell> elements only.");
}
warning("DAGMCUniverse <material_overrides> is deprecated. Use nested "
"<cell> elements under <dagmc_universe> instead.");
for (pugi::xml_node co :
node.child("material_overrides").children("cell_override")) {
int32_t cell_id = std::stoi(get_node_value(co, "id"));
std::istringstream iss(co.child("material_ids").text().get());
vector<int32_t> mats;
for (std::string s; iss >> s;) {
mats.push_back(s == "void" ? MATERIAL_VOID : std::stoi(s));
}
material_overrides.emplace(cell_id, mats);
}
}
initialize();
initialize(material_overrides, temperature_overrides, density_overrides);
}
DAGUniverse::DAGUniverse(
@ -110,9 +167,12 @@ DAGUniverse::DAGUniverse(std::shared_ptr<moab::DagMC> dagmc_ptr,
: dagmc_instance_(dagmc_ptr), filename_(filename),
adjust_geometry_ids_(auto_geom_ids), adjust_material_ids_(auto_mat_ids)
{
MaterialOverrides material_overrides;
TemperatureOverrides temperature_overrides;
DensityOverrides density_overrides;
set_id();
init_metadata();
init_geometry();
init_geometry(material_overrides, temperature_overrides, density_overrides);
}
void DAGUniverse::set_id()
@ -130,6 +190,15 @@ void DAGUniverse::set_id()
}
void DAGUniverse::initialize()
{
MaterialOverrides material_overrides;
TemperatureOverrides temperature_overrides;
initialize(material_overrides, temperature_overrides);
}
void DAGUniverse::initialize(const MaterialOverrides& material_overrides,
const TemperatureOverrides& temperature_overrides,
const DensityOverrides& density_overrides)
{
#ifdef OPENMC_UWUW_ENABLED
// read uwuw materials from the .h5m file if present
@ -140,7 +209,7 @@ void DAGUniverse::initialize()
init_metadata();
init_geometry();
init_geometry(material_overrides, temperature_overrides, density_overrides);
}
void DAGUniverse::init_dagmc()
@ -176,7 +245,9 @@ void DAGUniverse::init_metadata()
MB_CHK_ERR_CONT(rval);
}
void DAGUniverse::init_geometry()
void DAGUniverse::init_geometry(const MaterialOverrides& material_overrides,
const TemperatureOverrides& temperature_overrides,
const DensityOverrides& density_overrides)
{
moab::ErrorCode rval;
@ -202,6 +273,9 @@ void DAGUniverse::init_geometry()
: dagmc_instance_->id_by_index(3, c->dag_index());
c->universe_ = this->id_;
c->fill_ = C_NONE; // no fill, single universe
if (dagmc_instance_->is_implicit_complement(vol_handle)) {
c->name_ = "implicit complement";
}
auto in_map = model::cell_map.find(c->id_);
if (in_map == model::cell_map.end()) {
@ -230,16 +304,68 @@ void DAGUniverse::init_geometry()
if (mat_str == "graveyard") {
graveyard = vol_handle;
}
// material void checks
if (mat_str == "void" || mat_str == "vacuum" || mat_str == "graveyard") {
if (material_overrides.count(c->id_)) {
override_assign_material(c, material_overrides);
} else if (mat_str == "void" || mat_str == "vacuum" ||
mat_str == "graveyard") {
c->material_.push_back(MATERIAL_VOID);
} else if (uses_uwuw()) {
uwuw_assign_material(vol_handle, c);
} else {
if (material_overrides_.count(c->id_)) {
override_assign_material(c);
} else if (uses_uwuw()) {
uwuw_assign_material(vol_handle, c);
} else {
legacy_assign_material(mat_str, c);
legacy_assign_material(mat_str, c);
}
if (temperature_overrides.count(c->id_)) {
if (c->material_.empty() || c->material_[0] == MATERIAL_VOID) {
fatal_error(fmt::format("DAGMC cell {} was specified with a "
"temperature but no non-void material.",
c->id_));
}
c->sqrtkT_.clear();
const auto& temp_overrides = temperature_overrides.at(c->id_);
c->sqrtkT_.reserve(temp_overrides.size());
for (auto T : temp_overrides) {
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * T));
}
if (settings::verbosity >= 10) {
std::stringstream override_values;
for (size_t i = 0; i < temp_overrides.size(); ++i) {
if (i > 0) {
override_values << " ";
}
override_values << temp_overrides[i];
}
auto msg = fmt::format("Overriding DAGMC cell {} property "
"'temperature [K]' with value(s): {}",
c->id_, override_values.str());
write_message(msg, 10);
}
}
if (density_overrides.count(c->id_)) {
if (c->material_.empty() || c->material_[0] == MATERIAL_VOID) {
fatal_error(fmt::format("DAGMC cell {} was specified with a density "
"but no non-void material.",
c->id_));
}
// density_mult_ holds the true density until materials are finalized,
// at which point it is converted to a proper multiplier (same as CSG).
c->density_mult_ = density_overrides.at(c->id_);
if (settings::verbosity >= 10) {
const auto& dens = density_overrides.at(c->id_);
std::stringstream override_values;
for (size_t i = 0; i < dens.size(); ++i) {
if (i > 0)
override_values << " ";
override_values << dens[i];
}
write_message(fmt::format("Overriding DAGMC cell {} property "
"'density [g/cm³]' with value(s): {}",
c->id_, override_values.str()),
10);
}
}
@ -252,18 +378,21 @@ void DAGUniverse::init_geometry()
continue;
}
// assign cell temperature
const auto& mat = model::materials[model::material_map.at(c->material_[0])];
if (dagmc_instance_->has_prop(vol_handle, "temp")) {
rval = dagmc_instance_->prop_value(vol_handle, "temp", temp_value);
MB_CHK_ERR_CONT(rval);
double temp = std::stod(temp_value);
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * temp));
} else if (mat->temperature() > 0.0) {
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature()));
} else {
c->sqrtkT_.push_back(
std::sqrt(K_BOLTZMANN * settings::temperature_default));
// assign cell temperature if not explicitly overridden
if (c->sqrtkT_.empty()) {
const auto& mat =
model::materials[model::material_map.at(c->material_[0])];
if (dagmc_instance_->has_prop(vol_handle, "temp")) {
rval = dagmc_instance_->prop_value(vol_handle, "temp", temp_value);
MB_CHK_ERR_CONT(rval);
double temp = std::stod(temp_value);
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * temp));
} else if (mat->temperature() > 0.0) {
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature()));
} else {
c->sqrtkT_.push_back(
std::sqrt(K_BOLTZMANN * settings::temperature_default));
}
}
model::cells.emplace_back(std::move(c));
@ -630,7 +759,8 @@ void DAGUniverse::uwuw_assign_material(
#endif // OPENMC_UWUW_ENABLED
}
void DAGUniverse::override_assign_material(std::unique_ptr<DAGCell>& c) const
void DAGUniverse::override_assign_material(std::unique_ptr<DAGCell>& c,
const MaterialOverrides& material_overrides) const
{
// if Cell ID matches an override key, use it to override the material
// assignment else if UWUW is used, get the material assignment from the DAGMC
@ -638,17 +768,30 @@ void DAGUniverse::override_assign_material(std::unique_ptr<DAGCell>& c) const
// Notify User that an override is being applied on a DAGMCCell
write_message(fmt::format("Applying override for DAGMCCell {}", c->id_), 8);
const auto& mat_overrides = material_overrides.at(c->id_);
if (settings::verbosity >= 10) {
auto msg = fmt::format("Assigning DAGMC cell {} material(s) based on "
"override information (see input XML).",
c->id_);
std::stringstream override_values;
for (size_t i = 0; i < mat_overrides.size(); ++i) {
if (i > 0) {
override_values << " ";
}
if (mat_overrides[i] == MATERIAL_VOID) {
override_values << "void";
} else {
override_values << mat_overrides[i];
}
}
auto msg = fmt::format("Overriding DAGMC cell {} property 'material' "
"with value(s): {}",
c->id_, override_values.str());
write_message(msg, 10);
}
// Override the material assignment for each cell instance using the legacy
// assignement
for (auto mat_id : material_overrides_.at(c->id_)) {
if (model::material_map.find(mat_id) == model::material_map.end()) {
for (auto mat_id : mat_overrides) {
if (mat_id != MATERIAL_VOID &&
model::material_map.find(mat_id) == model::material_map.end()) {
fatal_error(fmt::format(
"Material with ID '{}' not found for DAGMC cell {}", mat_id, c->id_));
}

124
src/distribution.cpp
View file

@ -7,7 +7,9 @@
#include <numeric> // for accumulate
#include <stdexcept> // for runtime_error
#include <string> // for string, stod
#include <unordered_set>
#include "openmc/chain.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
@ -15,6 +17,10 @@
#include "openmc/random_lcg.h"
#include "openmc/xml_interface.h"
namespace {
std::unordered_set<std::string> decay_spectrum_missing_chain_nuclides;
}
namespace openmc {
//==============================================================================
@ -758,6 +764,8 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
dist = UPtrDist {new Tabular(node)};
} else if (type == "mixture") {
dist = UPtrDist {new Mixture(node)};
} else if (type == "decay_spectrum") {
dist = UPtrDist {new DecaySpectrum(node)};
} else if (type == "muir") {
openmc::fatal_error(
"'muir' distributions are now specified using the openmc.stats.muir() "
@ -768,4 +776,120 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
return dist;
}
//==============================================================================
// DecaySpectrum implementation
//==============================================================================
DecaySpectrum::DecaySpectrum(pugi::xml_node node)
{
// Read the region volume [cm^3] needed for absolute emission rate
if (!check_for_node(node, "volume"))
fatal_error("DecaySpectrum: 'volume' attribute is required.");
double volume = std::stod(get_node_value(node, "volume"));
// Read nuclide names and atom densities from XML
vector<int> nuclide_indices;
vector<double> atoms;
auto names = get_node_array<std::string>(node, "nuclides");
auto densities = get_node_array<double>(node, "parameters");
if (names.size() != densities.size()) {
fatal_error("DecaySpectrum nuclides and parameters must have the same "
"length.");
}
for (size_t i = 0; i < names.size(); ++i) {
const auto& name = names[i];
double density = densities[i];
// Look up nuclide in the depletion chain
auto it = data::chain_nuclide_map.find(name);
if (it == data::chain_nuclide_map.end()) {
if (decay_spectrum_missing_chain_nuclides.insert(name).second) {
warning("Nuclide '" + name +
"' appears in a DecaySpectrum source but is not present in "
"the depletion chain; it will be ignored.");
}
continue;
}
int nuclide_index = it->second;
const auto& chain_nuc = data::chain_nuclides[nuclide_index];
const Distribution* photon_dist = chain_nuc->photon_energy();
if (!photon_dist)
continue;
// Skip non-positive densities and warn if negative
if (density <= 0.0) {
if (density < 0.0) {
warning("Nuclide '" + name +
"' has a negative density in a DecaySpectrum source; it will "
"be ignored.");
}
continue;
}
// atoms = density [atom/b-cm] * 1e24 [b/cm^2] * volume [cm^3]
double atoms_i = density * 1.0e24 * volume;
nuclide_indices.push_back(nuclide_index);
atoms.push_back(atoms_i);
}
init(std::move(nuclide_indices), atoms);
}
void DecaySpectrum::init(
vector<int> nuclide_indices, const vector<double>& atoms)
{
if (nuclide_indices.size() != atoms.size()) {
fatal_error("DecaySpectrum nuclide index and atoms arrays must have "
"the same length.");
}
vector<double> probs;
probs.reserve(nuclide_indices.size());
for (size_t i = 0; i < nuclide_indices.size(); ++i) {
// Distribution integral is in [photons/s/atom]; multiplying by atoms gives
// the total emission rate [photons/s] for this nuclide.
const auto* dist =
data::chain_nuclides[nuclide_indices[i]]->photon_energy();
probs.push_back(atoms[i] * dist->integral());
}
nuclide_indices_ = std::move(nuclide_indices);
integral_ = std::accumulate(probs.begin(), probs.end(), 0.0);
if (nuclide_indices_.empty() || integral_ <= 0.0) {
fatal_error("DecaySpectrum source did not resolve any nuclides with decay "
"photon spectra and positive atom densities. Ensure "
"OPENMC_CHAIN_FILE is set and matches the nuclides in the "
"source definition.");
}
di_.assign(probs);
}
DecaySpectrum::Sample DecaySpectrum::sample_with_parent(uint64_t* seed) const
{
size_t idx = di_.sample(seed);
int parent_nuclide = nuclide_indices_[idx];
const auto* dist = data::chain_nuclides[parent_nuclide]->photon_energy();
auto [energy, weight] = dist->sample(seed);
return {energy, weight, parent_nuclide};
}
std::pair<double, double> DecaySpectrum::sample(uint64_t* seed) const
{
auto sample = sample_with_parent(seed);
return {sample.energy, sample.weight};
}
double DecaySpectrum::integral() const
{
return integral_;
}
double DecaySpectrum::sample_unbiased(uint64_t* seed) const
{
return sample_with_parent(seed).energy;
}
} // namespace openmc

51
src/event.cpp
View file

@ -1,6 +1,9 @@
#include "openmc/event.h"
#include "openmc/bank.h"
#include "openmc/error.h"
#include "openmc/material.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/timer.h"
@ -64,7 +67,8 @@ void process_init_events(int64_t n_particles, int64_t source_offset)
simulation::time_event_init.start();
#pragma omp parallel for schedule(runtime)
for (int64_t i = 0; i < n_particles; i++) {
initialize_history(simulation::particles[i], source_offset + i + 1);
initialize_particle_track(
simulation::particles[i], source_offset + i + 1, false);
dispatch_xs_event(i);
}
simulation::time_event_init.stop();
@ -136,7 +140,7 @@ void process_surface_crossing_events()
int64_t buffer_idx = simulation::surface_crossing_queue[i].idx;
Particle& p = simulation::particles[buffer_idx];
p.event_cross_surface();
p.event_revive_from_secondary();
p.event_check_limit_and_revive();
if (p.alive())
dispatch_xs_event(buffer_idx);
}
@ -155,7 +159,7 @@ void process_collision_events()
int64_t buffer_idx = simulation::collision_queue[i].idx;
Particle& p = simulation::particles[buffer_idx];
p.event_collide();
p.event_revive_from_secondary();
p.event_check_limit_and_revive();
if (p.alive())
dispatch_xs_event(buffer_idx);
}
@ -176,4 +180,45 @@ void process_death_events(int64_t n_particles)
simulation::time_event_death.stop();
}
void process_transport_events()
{
while (true) {
int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
simulation::calculate_nonfuel_xs_queue.size(),
simulation::advance_particle_queue.size(),
simulation::surface_crossing_queue.size(),
simulation::collision_queue.size()});
if (max == 0) {
break;
} else if (max == simulation::calculate_fuel_xs_queue.size()) {
process_calculate_xs_events(simulation::calculate_fuel_xs_queue);
} else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
} else if (max == simulation::advance_particle_queue.size()) {
process_advance_particle_events();
} else if (max == simulation::surface_crossing_queue.size()) {
process_surface_crossing_events();
} else if (max == simulation::collision_queue.size()) {
process_collision_events();
}
}
}
void process_init_secondary_events(int64_t n_particles, int64_t offset,
const SharedArray<SourceSite>& shared_secondary_bank)
{
simulation::time_event_init.start();
#pragma omp parallel for schedule(runtime)
for (int64_t i = 0; i < n_particles; i++) {
initialize_particle_track(simulation::particles[i], offset + i + 1, true);
const SourceSite& site = shared_secondary_bank[offset + i];
simulation::particles[i].event_revive_from_secondary(site);
if (simulation::particles[i].alive()) {
dispatch_xs_event(i);
}
}
simulation::time_event_init.stop();
}
} // namespace openmc

5
src/finalize.cpp
View file

@ -2,6 +2,7 @@
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/chain.h"
#include "openmc/cmfd_solver.h"
#include "openmc/collision_track.h"
#include "openmc/constants.h"
@ -44,6 +45,7 @@ void free_memory()
free_memory_photon();
free_memory_settings();
free_memory_thermal();
free_memory_chain();
library_clear();
nuclides_clear();
free_memory_source();
@ -147,6 +149,8 @@ int openmc_finalize()
settings::uniform_source_sampling = false;
settings::ufs_on = false;
settings::urr_ptables_on = true;
settings::use_decay_photons = false;
settings::use_shared_secondary_bank = false;
settings::verbosity = -1;
settings::weight_cutoff = 0.25;
settings::weight_survive = 1.0;
@ -217,6 +221,7 @@ int openmc_reset()
settings::cmfd_run = false;
simulation::n_lost_particles = 0;
simulation::simulation_tracks_completed = 0;
return 0;
}

4
src/ifp.cpp
View file

@ -32,13 +32,13 @@ void ifp(const Particle& p, int64_t idx)
{
if (is_beta_effective_or_both()) {
const auto& delayed_groups =
simulation::ifp_source_delayed_group_bank[p.current_work() - 1];
simulation::ifp_source_delayed_group_bank[p.current_work()];
simulation::ifp_fission_delayed_group_bank[idx] =
_ifp(p.delayed_group(), delayed_groups);
}
if (is_generation_time_or_both()) {
const auto& lifetimes =
simulation::ifp_source_lifetime_bank[p.current_work() - 1];
simulation::ifp_source_lifetime_bank[p.current_work()];
simulation::ifp_fission_lifetime_bank[idx] = _ifp(p.lifetime(), lifetimes);
}
}

78
src/initialize.cpp
View file

@ -161,7 +161,7 @@ void initialize_mpi(MPI_Comm intracomm)
// Create bank datatype
SourceSite b;
MPI_Aint disp[11];
MPI_Aint disp[14];
MPI_Get_address(&b.r, &disp[0]);
MPI_Get_address(&b.u, &disp[1]);
MPI_Get_address(&b.E, &disp[2]);
@ -173,14 +173,35 @@ void initialize_mpi(MPI_Comm intracomm)
MPI_Get_address(&b.parent_nuclide, &disp[8]);
MPI_Get_address(&b.parent_id, &disp[9]);
MPI_Get_address(&b.progeny_id, &disp[10]);
for (int i = 10; i >= 0; --i) {
MPI_Get_address(&b.wgt_born, &disp[11]);
MPI_Get_address(&b.wgt_ww_born, &disp[12]);
MPI_Get_address(&b.n_split, &disp[13]);
for (int i = 13; i >= 0; --i) {
disp[i] -= disp[0];
}
int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site);
// Block counts for each field
int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
// Types for each field
MPI_Datatype types[] = {
MPI_DOUBLE, // r (3 doubles)
MPI_DOUBLE, // u (3 doubles)
MPI_DOUBLE, // E
MPI_DOUBLE, // time
MPI_DOUBLE, // wgt
MPI_INT, // delayed_group
MPI_INT, // surf_id
MPI_INT, // particle (enum)
MPI_INT, // parent_nuclide
MPI_INT64_T, // parent_id
MPI_INT64_T, // progeny_id
MPI_DOUBLE, // wgt_born
MPI_DOUBLE, // wgt_ww_born
MPI_INT64_T // n_split
};
MPI_Type_create_struct(14, blocks, disp, types, &mpi::source_site);
MPI_Type_commit(&mpi::source_site);
CollisionTrackSite bc;
@ -300,6 +321,11 @@ int parse_command_line(int argc, char* argv[])
settings::run_mode = RunMode::VOLUME;
} else if (arg == "-s" || arg == "--threads") {
// Read number of threads
if (i + 1 >= argc) {
std::string msg {"Number of threads not specified."};
strcpy(openmc_err_msg, msg.c_str());
return OPENMC_E_INVALID_ARGUMENT;
}
i += 1;
#ifdef _OPENMP
@ -361,6 +387,28 @@ int parse_command_line(int argc, char* argv[])
return 0;
}
// TODO: Pulse-height tallies require per-history scoring across the full
// particle tree (parent + all descendants). The shared secondary bank
// transports each secondary as an independent Particle, breaking this
// assumption. A proper fix would defer pulse-height scoring: save
// (root_source_id, cell, pht_storage) per particle, then aggregate by
// root_source_id after all secondary generations complete before scoring
// into the histogram. For now, disable shared secondary when pulse-height
// tallies are present.
static void check_pulse_height_compatibility()
{
if (settings::use_shared_secondary_bank) {
for (const auto& t : model::tallies) {
if (t->type_ == TallyType::PULSE_HEIGHT) {
settings::use_shared_secondary_bank = false;
warning("Pulse-height tallies are not yet compatible with the shared "
"secondary bank. Disabling shared secondary bank.");
break;
}
}
}
}
bool read_model_xml()
{
std::string model_filename = settings::path_input;
@ -405,6 +453,10 @@ bool read_model_xml()
write_message(
fmt::format("Reading model XML file '{}' ...", model_filename), 5);
// Read chain data before settings so DecaySpectrum source distributions can
// resolve nuclides while sources are constructed.
read_chain_file_xml();
read_settings_xml(settings_root);
// If other XML files are present, display warning
@ -420,9 +472,6 @@ bool read_model_xml()
}
}
// Read data from chain file
read_chain_file_xml();
// Read materials and cross sections
if (!check_for_node(root, "materials")) {
fatal_error(fmt::format(
@ -454,6 +503,8 @@ bool read_model_xml()
if (check_for_node(root, "tallies"))
read_tallies_xml(root.child("tallies"));
check_pulse_height_compatibility();
// Initialize distribcell_filters
prepare_distribcell();
@ -475,14 +526,15 @@ bool read_model_xml()
void read_separate_xml_files()
{
// Read chain data before settings so DecaySpectrum source distributions can
// resolve nuclides while sources are constructed.
read_chain_file_xml();
read_settings_xml();
if (settings::run_mode != RunMode::PLOTTING) {
read_cross_sections_xml();
}
// Read data from chain file
read_chain_file_xml();
read_materials_xml();
read_geometry_xml();
@ -498,6 +550,8 @@ void read_separate_xml_files()
read_tallies_xml();
check_pulse_height_compatibility();
// Initialize distribcell_filters
prepare_distribcell();

4
src/mesh.cpp
View file

@ -3670,7 +3670,7 @@ Position LibMesh::sample_element(int32_t bin, uint64_t* seed) const
// Get tet vertex coordinates from LibMesh
std::array<Position, 4> tet_verts;
for (int i = 0; i < elem.n_nodes(); i++) {
auto node_ref = elem.node_ref(i);
const auto& node_ref = elem.node_ref(i);
tet_verts[i] = {node_ref(0), node_ref(1), node_ref(2)};
}
// Samples position within tet using Barycentric coordinates
@ -3700,7 +3700,7 @@ int LibMesh::n_vertices() const
Position LibMesh::vertex(int vertex_id) const
{
const auto node_ref = m_->node_ref(vertex_id);
const auto& node_ref = m_->node_ref(vertex_id);
if (length_multiplier_ > 0.0) {
return length_multiplier_ * Position(node_ref(0), node_ref(1), node_ref(2));
} else {

8
src/output.cpp
View file

@ -501,6 +501,14 @@ void print_runtime()
show_rate("Calculation Rate (inactive)", speed_inactive);
}
show_rate("Calculation Rate (active)", speed_active);
// Display track rate when weight windows are enabled
if (settings::weight_windows_on) {
double speed_tracks =
simulation::simulation_tracks_completed / time_active.elapsed();
fmt::print(
" {:<33} = {:.6} tracks/second\n", "Track Rate (active)", speed_tracks);
}
}
//==============================================================================

178
src/particle.cpp
View file

@ -94,7 +94,7 @@ bool Particle::create_secondary(
// Increment number of secondaries created (for ParticleProductionFilter)
n_secondaries()++;
auto& bank = secondary_bank().emplace_back();
SourceSite bank;
bank.particle = type;
bank.wgt = wgt;
bank.r = r();
@ -102,12 +102,21 @@ bool Particle::create_secondary(
bank.E = settings::run_CE ? E : g();
bank.time = time();
bank_second_E() += bank.E;
bank.parent_id = current_work();
if (settings::use_shared_secondary_bank) {
bank.progeny_id = n_progeny()++;
}
bank.wgt_born = wgt_born();
bank.wgt_ww_born = wgt_ww_born();
bank.n_split = n_split();
local_secondary_bank().emplace_back(bank);
return true;
}
void Particle::split(double wgt)
{
auto& bank = secondary_bank().emplace_back();
SourceSite bank;
bank.particle = type();
bank.wgt = wgt;
bank.r = r();
@ -122,6 +131,16 @@ void Particle::split(double wgt)
int surf_id = model::surfaces[surface_index()]->id_;
bank.surf_id = (surface() > 0) ? surf_id : -surf_id;
}
bank.wgt_born = wgt_born();
bank.wgt_ww_born = wgt_ww_born();
bank.n_split = n_split();
bank.parent_id = current_work();
if (settings::use_shared_secondary_bank) {
bank.progeny_id = n_progeny()++;
}
local_secondary_bank().emplace_back(bank);
}
void Particle::from_source(const SourceSite* src)
@ -168,6 +187,10 @@ void Particle::from_source(const SourceSite* src)
int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1;
surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
}
wgt_born() = src->wgt_born;
wgt_ww_born() = src->wgt_ww_born;
n_split() = src->n_split;
}
void Particle::event_calculate_xs()
@ -449,61 +472,72 @@ void Particle::event_collide()
#endif
}
void Particle::event_revive_from_secondary()
void Particle::event_revive_from_secondary(const SourceSite& site)
{
// Write final position for the previous track (skip if this is a freshly
// constructed particle with no prior track, e.g., Phase 2 of shared
// secondary transport)
if (write_track() && n_event() > 0) {
write_particle_track(*this);
}
from_source(&site);
n_event() = 0;
if (!settings::use_shared_secondary_bank) {
n_tracks()++;
}
bank_second_E() = 0.0;
// Subtract secondary particle energy from interim pulse-height results.
// In shared secondary mode, this subtraction was already done on the parent
// particle during create_secondary(), so skip it here.
if (!settings::use_shared_secondary_bank &&
!model::active_pulse_height_tallies.empty() && this->type().is_photon()) {
// Since the birth cell of the particle has not been set we
// have to determine it before the energy of the secondary particle can be
// removed from the pulse-height of this cell.
if (lowest_coord().cell() == C_NONE) {
bool verbose = settings::verbosity >= 10 || trace();
if (!exhaustive_find_cell(*this, verbose)) {
mark_as_lost("Could not find the cell containing particle " +
std::to_string(id()));
return;
}
// Set birth cell attribute
if (cell_born() == C_NONE)
cell_born() = lowest_coord().cell();
// Initialize last cells from current cell
for (int j = 0; j < n_coord(); ++j) {
cell_last(j) = coord(j).cell();
}
n_coord_last() = n_coord();
}
pht_secondary_particles();
}
// Enter new particle in particle track file
if (write_track())
add_particle_track(*this);
}
void Particle::event_check_limit_and_revive()
{
// If particle has too many events, display warning and kill it
++n_event();
n_event()++;
if (n_event() == settings::max_particle_events) {
warning("Particle " + std::to_string(id()) +
" underwent maximum number of events.");
wgt() = 0.0;
}
// Check for secondary particles if this particle is dead
if (!alive()) {
// Write final position for this particle
if (write_track()) {
write_particle_track(*this);
}
// If no secondary particles, break out of event loop
if (secondary_bank().empty())
return;
from_source(&secondary_bank().back());
secondary_bank().pop_back();
n_event() = 0;
bank_second_E() = 0.0;
// Subtract secondary particle energy from interim pulse-height results
if (!model::active_pulse_height_tallies.empty() &&
this->type().is_photon()) {
// Since the birth cell of the particle has not been set we
// have to determine it before the energy of the secondary particle can be
// removed from the pulse-height of this cell.
if (lowest_coord().cell() == C_NONE) {
bool verbose = settings::verbosity >= 10 || trace();
if (!exhaustive_find_cell(*this, verbose)) {
mark_as_lost("Could not find the cell containing particle " +
std::to_string(id()));
return;
}
// Set birth cell attribute
if (cell_born() == C_NONE)
cell_born() = lowest_coord().cell();
// Initialize last cells from current cell
for (int j = 0; j < n_coord(); ++j) {
cell_last(j) = coord(j).cell();
}
n_coord_last() = n_coord();
}
pht_secondary_particles();
}
// Enter new particle in particle track file
if (write_track())
add_particle_track(*this);
// In non-shared-secondary mode, revive from local secondary bank
if (!alive() && !settings::use_shared_secondary_bank &&
!local_secondary_bank().empty()) {
SourceSite& site = local_secondary_bank().back();
event_revive_from_secondary(site);
local_secondary_bank().pop_back();
}
}
@ -515,6 +549,7 @@ void Particle::event_death()
// Finish particle track output.
if (write_track()) {
write_particle_track(*this);
finalize_particle_track(*this);
}
@ -538,11 +573,17 @@ void Particle::event_death()
score_pulse_height_tally(*this, model::active_pulse_height_tallies);
}
// Accumulate track count for this particle history
if (!settings::use_shared_secondary_bank) {
#pragma omp atomic
simulation::simulation_tracks_completed += n_tracks();
}
// Record the number of progeny created by this particle.
// This data will be used to efficiently sort the fission bank.
if (settings::run_mode == RunMode::EIGENVALUE) {
int64_t offset = id() - 1 - simulation::work_index[mpi::rank];
simulation::progeny_per_particle[offset] = n_progeny();
if (settings::run_mode == RunMode::EIGENVALUE ||
settings::use_shared_secondary_bank) {
simulation::progeny_per_particle[current_work()] = n_progeny();
}
}
@ -863,28 +904,27 @@ void Particle::write_restart() const
write_dataset(file_id, "id", id());
write_dataset(file_id, "type", type().pdg_number());
// Get source site data for the particle that got lost
int64_t i = current_work();
SourceSite site;
if (settings::run_mode == RunMode::EIGENVALUE) {
// take source data from primary bank for eigenvalue simulation
write_dataset(file_id, "weight", simulation::source_bank[i - 1].wgt);
write_dataset(file_id, "energy", simulation::source_bank[i - 1].E);
write_dataset(file_id, "xyz", simulation::source_bank[i - 1].r);
write_dataset(file_id, "uvw", simulation::source_bank[i - 1].u);
write_dataset(file_id, "time", simulation::source_bank[i - 1].time);
site = simulation::source_bank[i];
} else if (settings::run_mode == RunMode::FIXED_SOURCE &&
settings::use_shared_secondary_bank &&
i < simulation::shared_secondary_bank_read.size()) {
site = simulation::shared_secondary_bank_read[i];
} else if (settings::run_mode == RunMode::FIXED_SOURCE) {
// re-sample using rng random number seed used to generate source particle
int64_t id = (simulation::total_gen + overall_generation() - 1) *
settings::n_particles +
simulation::work_index[mpi::rank] + i;
// Re-sample using the same seed used to generate the source particle.
// current_work() is 0-indexed, compute_particle_id expects 1-indexed.
int64_t id = compute_transport_seed(compute_particle_id(i + 1));
uint64_t seed = init_seed(id, STREAM_SOURCE);
// re-sample source site
auto site = sample_external_source(&seed);
write_dataset(file_id, "weight", site.wgt);
write_dataset(file_id, "energy", site.E);
write_dataset(file_id, "xyz", site.r);
write_dataset(file_id, "uvw", site.u);
write_dataset(file_id, "time", site.time);
site = sample_external_source(&seed);
}
write_dataset(file_id, "weight", site.wgt);
write_dataset(file_id, "energy", site.E);
write_dataset(file_id, "xyz", site.r);
write_dataset(file_id, "uvw", site.u);
write_dataset(file_id, "time", site.time);
// Close file
file_close(file_id);

23
src/particle_restart.cpp
View file

@ -1,6 +1,7 @@
#include "openmc/particle_restart.h"
#include "openmc/array.h"
#include "openmc/bank.h"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/mgxs_interface.h"
@ -106,20 +107,16 @@ void run_particle_restart()
// Set all tallies to 0 for now (just tracking errors)
model::tallies.clear();
// Compute random number seed
int64_t particle_seed;
switch (previous_run_mode) {
case RunMode::EIGENVALUE:
case RunMode::FIXED_SOURCE:
particle_seed = (simulation::total_gen + overall_generation() - 1) *
settings::n_particles +
p.id();
break;
default:
throw std::runtime_error {
"Unexpected run mode: " +
std::to_string(static_cast<int>(previous_run_mode))};
// Allocate progeny_per_particle if needed for shared secondary mode
// (event_death() writes to this array). Set current_work to 0 since we
// only have one particle being restarted.
if (settings::use_shared_secondary_bank) {
p.current_work() = 0;
simulation::progeny_per_particle.resize(1, 0);
}
// Compute random number seed
int64_t particle_seed = compute_transport_seed(p.id());
init_particle_seeds(particle_seed, p.seeds());
// Force calculation of cross-sections by setting last energy to zero

25
src/physics.cpp
View file

@ -44,7 +44,7 @@ void collision(Particle& p)
{
// Add to collision counter for particle
++(p.n_collision());
p.secondary_bank_index() = p.secondary_bank().size();
p.secondary_bank_index() = p.local_secondary_bank().size();
// Sample reaction for the material the particle is in
switch (p.type().pdg_number()) {
@ -127,7 +127,8 @@ void sample_neutron_reaction(Particle& p)
// Make sure particle population doesn't grow out of control for
// subcritical multiplication problems.
if (p.secondary_bank().size() >= settings::max_secondaries) {
if (p.local_secondary_bank().size() >= settings::max_secondaries &&
!settings::use_shared_secondary_bank) {
fatal_error(
"The secondary particle bank appears to be growing without "
"bound. You are likely running a subcritical multiplication problem "
@ -228,7 +229,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
}
// Set parent and progeny IDs
site.parent_id = p.id();
site.parent_id = p.current_work();
site.progeny_id = p.n_progeny()++;
// Store fission site in bank
@ -253,7 +254,10 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
ifp(p, idx);
}
} else {
p.secondary_bank().push_back(site);
site.wgt_born = p.wgt_born();
site.wgt_ww_born = p.wgt_ww_born();
site.n_split = p.n_split();
p.local_secondary_bank().push_back(site);
p.n_secondaries()++;
}
@ -1223,9 +1227,20 @@ void sample_secondary_photons(Particle& p, int i_nuclide)
// Create the secondary photon
bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon());
// Pre-add photon energy to pht_storage so pht_secondary_particles()
// subtraction results in net zero
if (created_photon && !model::active_pulse_height_tallies.empty()) {
auto it = std::find(model::pulse_height_cells.begin(),
model::pulse_height_cells.end(), p.lowest_coord().cell());
if (it != model::pulse_height_cells.end()) {
int index = std::distance(model::pulse_height_cells.begin(), it);
p.pht_storage()[index] += E;
}
}
// Tag secondary particle with parent nuclide
if (created_photon && settings::use_decay_photons) {
p.secondary_bank().back().parent_nuclide =
p.local_secondary_bank().back().parent_nuclide =
rx->products_[i_product].parent_nuclide_;
}
}

9
src/physics_mg.cpp
View file

@ -27,7 +27,7 @@ void collision_mg(Particle& p)
{
// Add to the collision counter for the particle
p.n_collision()++;
p.secondary_bank_index() = p.secondary_bank().size();
p.secondary_bank_index() = p.local_secondary_bank().size();
// Sample the reaction type
sample_reaction(p);
@ -179,7 +179,7 @@ void create_fission_sites(Particle& p)
}
// Set parent and progeny ID
site.parent_id = p.id();
site.parent_id = p.current_work();
site.progeny_id = p.n_progeny()++;
// Store fission site in bank
@ -200,7 +200,10 @@ void create_fission_sites(Particle& p)
break;
}
} else {
p.secondary_bank().push_back(site);
site.wgt_born = p.wgt_born();
site.wgt_ww_born = p.wgt_ww_born();
site.n_split = p.n_split();
p.local_secondary_bank().push_back(site);
p.n_secondaries()++;
}

202
src/plot.cpp
View file

@ -33,6 +33,7 @@
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
#include "openmc/tallies/filter.h"
namespace openmc {
@ -44,10 +45,12 @@ constexpr int PLOT_LEVEL_LOWEST {-1}; //!< lower bound on plot universe level
constexpr int32_t NOT_FOUND {-2};
constexpr int32_t OVERLAP {-3};
IdData::IdData(size_t h_res, size_t v_res) : data_({v_res, h_res, 3}, NOT_FOUND)
IdData::IdData(size_t h_res, size_t v_res, bool /*include_filter*/)
: data_({v_res, h_res, 3}, NOT_FOUND)
{}
void IdData::set_value(size_t y, size_t x, const GeometryState& p, int level)
void IdData::set_value(size_t y, size_t x, const Particle& p, int level,
Filter* /*filter*/, FilterMatch* /*match*/)
{
// set cell data
if (p.n_coord() <= level) {
@ -64,7 +67,6 @@ void IdData::set_value(size_t y, size_t x, const GeometryState& p, int level)
Cell* c = model::cells.at(p.lowest_coord().cell()).get();
if (p.material() == MATERIAL_VOID) {
data_(y, x, 2) = MATERIAL_VOID;
return;
} else if (c->type_ == Fill::MATERIAL) {
Material* m = model::materials.at(p.material()).get();
data_(y, x, 2) = m->id_;
@ -77,12 +79,12 @@ void IdData::set_overlap(size_t y, size_t x)
data_(y, x, k) = OVERLAP;
}
PropertyData::PropertyData(size_t h_res, size_t v_res)
PropertyData::PropertyData(size_t h_res, size_t v_res, bool /*include_filter*/)
: data_({v_res, h_res, 2}, NOT_FOUND)
{}
void PropertyData::set_value(
size_t y, size_t x, const GeometryState& p, int level)
void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level,
Filter* /*filter*/, FilterMatch* /*match*/)
{
Cell* c = model::cells.at(p.lowest_coord().cell()).get();
data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN;
@ -97,6 +99,74 @@ void PropertyData::set_overlap(size_t y, size_t x)
data_(y, x) = OVERLAP;
}
//==============================================================================
// RasterData implementation
//==============================================================================
RasterData::RasterData(size_t h_res, size_t v_res, bool include_filter)
: id_data_({v_res, h_res, include_filter ? 4u : 3u}, NOT_FOUND),
property_data_({v_res, h_res, 2}, static_cast<double>(NOT_FOUND)),
include_filter_(include_filter)
{}
void RasterData::set_value(size_t y, size_t x, const Particle& p, int level,
Filter* filter, FilterMatch* match)
{
// set cell data
if (p.n_coord() <= level) {
id_data_(y, x, 0) = NOT_FOUND;
id_data_(y, x, 1) = NOT_FOUND;
} else {
id_data_(y, x, 0) = model::cells.at(p.coord(level).cell())->id_;
id_data_(y, x, 1) = level == p.n_coord() - 1
? p.cell_instance()
: cell_instance_at_level(p, level);
}
// set material data
Cell* c = model::cells.at(p.lowest_coord().cell()).get();
if (p.material() == MATERIAL_VOID) {
id_data_(y, x, 2) = MATERIAL_VOID;
} else if (c->type_ == Fill::MATERIAL) {
Material* m = model::materials.at(p.material()).get();
id_data_(y, x, 2) = m->id_;
}
// set filter index (only if filter is being used)
if (include_filter_ && filter) {
filter->get_all_bins(p, TallyEstimator::COLLISION, *match);
if (match->bins_.empty()) {
id_data_(y, x, 3) = -1;
} else {
id_data_(y, x, 3) = match->bins_[0];
}
match->bins_.clear();
match->weights_.clear();
}
// set temperature (in K)
property_data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN;
// 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_;
}
}
void RasterData::set_overlap(size_t y, size_t x)
{
// Set cell, instance, and material to OVERLAP, but preserve filter bin
id_data_(y, x, 0) = OVERLAP;
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;
}
//==============================================================================
// Global variables
//==============================================================================
@ -450,6 +520,22 @@ void Plot::set_width(pugi::xml_node plot_node)
if (pl_width.size() == 2) {
width_.x = pl_width[0];
width_.y = pl_width[1];
switch (basis_) {
case PlotBasis::xy:
u_span_ = {width_.x, 0.0, 0.0};
v_span_ = {0.0, width_.y, 0.0};
break;
case PlotBasis::xz:
u_span_ = {width_.x, 0.0, 0.0};
v_span_ = {0.0, 0.0, width_.y};
break;
case PlotBasis::yz:
u_span_ = {0.0, width_.x, 0.0};
v_span_ = {0.0, 0.0, width_.y};
break;
default:
UNREACHABLE();
}
} else {
fatal_error(
fmt::format("<width> must be length 2 in slice plot {}", id()));
@ -765,7 +851,7 @@ Plot::Plot(pugi::xml_node plot_node, PlotType type)
set_width(plot_node);
set_meshlines(plot_node);
slice_level_ = level_; // Copy level employed in SlicePlotBase::get_map
slice_color_overlaps_ = color_overlaps_;
show_overlaps_ = color_overlaps_;
}
//==============================================================================
@ -862,23 +948,39 @@ void Plot::draw_mesh_lines(ImageData& data) const
rgb = meshlines_color_;
int ax1, ax2;
Position expected_u {};
Position expected_v {};
switch (basis_) {
case PlotBasis::xy:
ax1 = 0;
ax2 = 1;
expected_u = {width_[0], 0.0, 0.0};
expected_v = {0.0, width_[1], 0.0};
break;
case PlotBasis::xz:
ax1 = 0;
ax2 = 2;
expected_u = {width_[0], 0.0, 0.0};
expected_v = {0.0, 0.0, width_[1]};
break;
case PlotBasis::yz:
ax1 = 1;
ax2 = 2;
expected_u = {0.0, width_[0], 0.0};
expected_v = {0.0, 0.0, width_[1]};
break;
default:
UNREACHABLE();
}
// Meshlines rely on axis-aligned indexing in global coordinates.
constexpr double rel_tol {1e-12};
double span_tol = rel_tol * (1.0 + u_span_.norm() + v_span_.norm());
if ((u_span_ - expected_u).norm() > span_tol ||
(v_span_ - expected_v).norm() > span_tol) {
fatal_error("Meshlines are only supported for axis-aligned slice plots.");
}
Position ll_plot {origin_};
Position ur_plot {origin_};
@ -1008,11 +1110,11 @@ void Plot::create_voxel() const
voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace);
SlicePlotBase pltbase;
pltbase.width_ = width_;
pltbase.origin_ = origin_;
pltbase.basis_ = PlotBasis::xy;
pltbase.u_span_ = {width_.x, 0.0, 0.0};
pltbase.v_span_ = {0.0, width_.y, 0.0};
pltbase.pixels() = pixels();
pltbase.slice_color_overlaps_ = color_overlaps_;
pltbase.show_overlaps_ = color_overlaps_;
ProgressBar pb;
for (int z = 0; z < pixels()[2]; z++) {
@ -1794,6 +1896,12 @@ void PhongRay::on_intersection()
extern "C" int openmc_id_map(const void* plot, int32_t* data_out)
{
static bool warned {false};
if (!warned) {
warning("openmc_id_map is deprecated and will be removed in a future "
"release. Use openmc_slice_data.");
warned = true;
}
auto plt = reinterpret_cast<const SlicePlotBase*>(plot);
if (!plt) {
@ -1801,7 +1909,7 @@ extern "C" int openmc_id_map(const void* plot, int32_t* data_out)
return OPENMC_E_INVALID_ARGUMENT;
}
if (plt->slice_color_overlaps_ && model::overlap_check_count.size() == 0) {
if (plt->show_overlaps_ && model::overlap_check_count.size() == 0) {
model::overlap_check_count.resize(model::cells.size());
}
@ -1815,14 +1923,20 @@ extern "C" int openmc_id_map(const void* plot, int32_t* data_out)
extern "C" int openmc_property_map(const void* plot, double* data_out)
{
static bool warned {false};
if (!warned) {
warning("openmc_property_map is deprecated and will be removed in a future "
"release. Use openmc_slice_data.");
warned = true;
}
auto plt = reinterpret_cast<const SlicePlotBase*>(plot);
if (!plt) {
set_errmsg("Invalid slice pointer passed to openmc_id_map");
set_errmsg("Invalid slice pointer passed to openmc_property_map");
return OPENMC_E_INVALID_ARGUMENT;
}
if (plt->slice_color_overlaps_ && model::overlap_check_count.size() == 0) {
if (plt->show_overlaps_ && model::overlap_check_count.size() == 0) {
model::overlap_check_count.resize(model::cells.size());
}
@ -1834,6 +1948,68 @@ extern "C" int openmc_property_map(const void* plot, double* data_out)
return 0;
}
extern "C" int openmc_slice_data(const double origin[3], const double u_span[3],
const double v_span[3], const size_t pixels[2], bool color_overlaps,
int level, int32_t filter_index, int32_t* geom_data, double* property_data)
{
// Validate span vectors
Direction u_span_pos {u_span[0], u_span[1], u_span[2]};
Direction v_span_pos {v_span[0], v_span[1], v_span[2]};
double u_norm = u_span_pos.norm();
double v_norm = v_span_pos.norm();
if (u_norm == 0.0 || v_norm == 0.0) {
set_errmsg("Slice span vectors must be non-zero.");
return OPENMC_E_INVALID_ARGUMENT;
}
constexpr double ORTHO_REL_TOL = 1e-10;
double dot = u_span_pos.dot(v_span_pos);
if (std::abs(dot) > ORTHO_REL_TOL * u_norm * v_norm) {
set_errmsg("Slice span vectors must be orthogonal.");
return OPENMC_E_INVALID_ARGUMENT;
}
// Validate filter index if provided
if (filter_index >= 0) {
if (int err = verify_filter(filter_index))
return err;
}
// Initialize overlap check vector if needed
if (color_overlaps && model::overlap_check_count.size() == 0) {
model::overlap_check_count.resize(model::cells.size());
}
try {
// Create a temporary SlicePlotBase object to reuse get_map logic
SlicePlotBase plot_params;
plot_params.origin_ = Position {origin[0], origin[1], origin[2]};
plot_params.u_span_ = u_span_pos;
plot_params.v_span_ = v_span_pos;
plot_params.pixels_[0] = pixels[0];
plot_params.pixels_[1] = pixels[1];
plot_params.show_overlaps_ = color_overlaps;
plot_params.slice_level_ = level;
// 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
if (property_data != nullptr) {
std::copy(
data.property_data_.begin(), data.property_data_.end(), property_data);
}
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
return 0;
}
extern "C" int openmc_get_plot_index(int32_t id, int32_t* index)
{
auto it = model::plot_map.find(id);

22
src/settings.cpp
View file

@ -83,6 +83,7 @@ bool uniform_source_sampling {false};
bool ufs_on {false};
bool urr_ptables_on {true};
bool use_decay_photons {false};
bool use_shared_secondary_bank {false};
bool weight_windows_on {false};
bool weight_window_checkpoint_surface {false};
bool weight_window_checkpoint_collision {true};
@ -1320,6 +1321,27 @@ void read_settings_xml(pugi::xml_node root)
settings::use_decay_photons =
get_node_value_bool(root, "use_decay_photons");
}
// If weight windows are on, also enable shared secondary bank (unless
// explicitly disabled by user).
if (check_for_node(root, "shared_secondary_bank")) {
bool val = get_node_value_bool(root, "shared_secondary_bank");
if (val && run_mode == RunMode::EIGENVALUE) {
warning(
"Shared secondary bank is not supported in eigenvalue calculations. "
"Setting will be ignored.");
} else {
settings::use_shared_secondary_bank = val;
}
} else if (settings::weight_windows_on) {
if (run_mode == RunMode::EIGENVALUE) {
warning(
"Shared secondary bank is not supported in eigenvalue calculations. "
"Particle local secondary banks will be used instead.");
} else if (run_mode == RunMode::FIXED_SOURCE) {
settings::use_shared_secondary_bank = true;
}
}
}
void free_memory_settings()

379
src/simulation.cpp
View file

@ -86,7 +86,7 @@ int openmc_simulation_init()
}
// Determine how much work each process should do
calculate_work();
calculate_work(settings::n_particles);
// Allocate source, fission and surface source banks.
allocate_banks();
@ -214,6 +214,20 @@ int openmc_simulation_finalize()
// Stop timers and show timing statistics
simulation::time_finalize.stop();
simulation::time_total.stop();
#ifdef OPENMC_MPI
// Reduce track count across ranks for correct reporting. In shared secondary
// bank mode, all ranks already have the global count; in non-shared mode,
// each rank only has its own count.
if (settings::weight_windows_on && !settings::use_shared_secondary_bank) {
int64_t total_tracks;
MPI_Reduce(&simulation::simulation_tracks_completed, &total_tracks, 1,
MPI_INT64_T, MPI_SUM, 0, mpi::intracomm);
if (mpi::master)
simulation::simulation_tracks_completed = total_tracks;
}
#endif
if (mpi::master) {
if (settings::solver_type != SolverType::RANDOM_RAY) {
if (settings::verbosity >= 6)
@ -254,9 +268,17 @@ int openmc_next_batch(int* status)
// Transport loop
if (settings::event_based) {
transport_event_based();
if (settings::use_shared_secondary_bank) {
transport_event_based_shared_secondary();
} else {
transport_event_based();
}
} else {
transport_history_based();
if (settings::use_shared_secondary_bank) {
transport_history_based_shared_secondary();
} else {
transport_history_based();
}
}
// Accumulate time for transport
@ -324,6 +346,8 @@ const RegularMesh* ufs_mesh {nullptr};
vector<double> k_generation;
vector<int64_t> work_index;
int64_t simulation_tracks_completed {0};
} // namespace simulation
//==============================================================================
@ -547,7 +571,7 @@ void finalize_generation()
// If using shared memory, stable sort the fission bank (by parent IDs)
// so as to allow for reproducibility regardless of which order particles
// are run in.
sort_fission_bank();
sort_bank(simulation::fission_bank, true);
// Distribute fission bank across processors evenly
synchronize_bank();
@ -571,26 +595,35 @@ void finalize_generation()
}
}
void initialize_history(Particle& p, int64_t index_source)
void sample_source_particle(Particle& p, int64_t index_source)
{
// set defaults
// Sample a particle from the source bank
if (settings::run_mode == RunMode::EIGENVALUE) {
// set defaults for eigenvalue simulations from primary bank
p.from_source(&simulation::source_bank[index_source - 1]);
} else if (settings::run_mode == RunMode::FIXED_SOURCE) {
// initialize random number seed
int64_t id = (simulation::total_gen + overall_generation() - 1) *
settings::n_particles +
simulation::work_index[mpi::rank] + index_source;
int64_t id = compute_transport_seed(compute_particle_id(index_source));
uint64_t seed = init_seed(id, STREAM_SOURCE);
// sample from external source distribution or custom library then set
auto site = sample_external_source(&seed);
p.from_source(&site);
}
p.current_work() = index_source;
}
void initialize_particle_track(
Particle& p, int64_t index_source, bool is_secondary)
{
// Note: index_source is 1-based (first particle = 1), but current_work() is
// stored as 0-based for direct use as an array index into
// progeny_per_particle, source_bank, ifp banks, etc.
if (!is_secondary) {
sample_source_particle(p, index_source);
}
p.current_work() = index_source - 1;
// set identifier for particle
p.id() = simulation::work_index[mpi::rank] + index_source;
p.id() = compute_particle_id(index_source);
// set progeny count to zero
p.n_progeny() = 0;
@ -598,6 +631,9 @@ void initialize_history(Particle& p, int64_t index_source)
// Reset particle event counter
p.n_event() = 0;
// Initialize track counter (1 for this primary/secondary track)
p.n_tracks() = 1;
// Reset split counter
p.n_split() = 0;
@ -611,9 +647,7 @@ void initialize_history(Particle& p, int64_t index_source)
std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0);
// set random number seed
int64_t particle_seed =
(simulation::total_gen + overall_generation() - 1) * settings::n_particles +
p.id();
int64_t particle_seed = compute_transport_seed(p.id());
init_particle_seeds(particle_seed, p.seeds());
// set particle trace
@ -627,17 +661,21 @@ void initialize_history(Particle& p, int64_t index_source)
p.write_track() = check_track_criteria(p);
// Set the particle's initial weight window value.
p.wgt_ww_born() = -1.0;
apply_weight_windows(p);
if (!is_secondary) {
p.wgt_ww_born() = -1.0;
apply_weight_windows(p);
}
// Display message if high verbosity or trace is on
if (settings::verbosity >= 9 || p.trace()) {
write_message("Simulating Particle {}", p.id());
}
// Add particle's starting weight to count for normalizing tallies later
// Add particle's starting weight to count for normalizing tallies later
if (!is_secondary) {
#pragma omp atomic
simulation::total_weight += p.wgt();
simulation::total_weight += p.wgt();
}
// Force calculation of cross-sections by setting last energy to zero
if (settings::run_CE) {
@ -655,13 +693,34 @@ int overall_generation()
return settings::gen_per_batch * (current_batch - 1) + current_gen;
}
void calculate_work()
int64_t compute_particle_id(int64_t index_source)
{
if (settings::use_shared_secondary_bank) {
return simulation::work_index[mpi::rank] + index_source +
simulation::simulation_tracks_completed;
} else {
return simulation::work_index[mpi::rank] + index_source;
}
}
int64_t compute_transport_seed(int64_t particle_id)
{
if (settings::use_shared_secondary_bank) {
return particle_id;
} else {
return (simulation::total_gen + overall_generation() - 1) *
settings::n_particles +
particle_id;
}
}
void calculate_work(int64_t n_particles)
{
// Determine minimum amount of particles to simulate on each processor
int64_t min_work = settings::n_particles / mpi::n_procs;
int64_t min_work = n_particles / mpi::n_procs;
// Determine number of processors that have one extra particle
int64_t remainder = settings::n_particles % mpi::n_procs;
int64_t remainder = n_particles % mpi::n_procs;
int64_t i_bank = 0;
simulation::work_index.resize(mpi::n_procs + 1);
@ -816,7 +875,7 @@ void transport_history_based_single_particle(Particle& p)
p.event_collide();
}
}
p.event_revive_from_secondary();
p.event_check_limit_and_revive();
}
p.event_death();
}
@ -826,11 +885,137 @@ void transport_history_based()
#pragma omp parallel for schedule(runtime)
for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
Particle p;
initialize_history(p, i_work);
initialize_particle_track(p, i_work, false);
transport_history_based_single_particle(p);
}
}
// The shared secondary bank transport algorithm works in two phases. In the
// first phase, all primary particles are sampled then transported, and their
// secondary particles are deposited into a shared secondary bank. The second
// phase occurs in a loop, where all secondary tracks in the shared secondary
// bank are transported. Any secondary particles generated during this phase are
// deposited back into the shared secondary bank. The shared secondary bank is
// sorted for consistent ordering and load balanced across MPI ranks. This loop
// continues until there are no more secondary tracks left to transport.
void transport_history_based_shared_secondary()
{
// Clear shared secondary banks from any prior use
simulation::shared_secondary_bank_read.clear();
simulation::shared_secondary_bank_write.clear();
if (mpi::master) {
write_message(fmt::format(" Primary source particles: {}",
settings::n_particles),
6);
}
simulation::progeny_per_particle.resize(simulation::work_per_rank);
std::fill(simulation::progeny_per_particle.begin(),
simulation::progeny_per_particle.end(), 0);
// Phase 1: Transport primary particles and deposit first generation of
// secondaries in the shared secondary bank
#pragma omp parallel
{
vector<SourceSite> thread_bank;
#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);
}
}
}
simulation::simulation_tracks_completed += settings::n_particles;
// Phase 2: Now that the secondary bank has been populated, enter loop over
// all secondary generations
int n_generation_depth = 1;
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
// all MPI ranks.
alive_secondary = synchronize_global_secondary_bank(
simulation::shared_secondary_bank_write);
// Recalculate work for each MPI rank based on number of alive secondary
// tracks
calculate_work(alive_secondary);
// Display the number of secondary tracks in this generation. This
// is useful for user monitoring so as to see if the secondary population is
// exploding and to determine how many generations of secondaries are being
// transported.
if (mpi::master) {
write_message(fmt::format(" Secondary generation {:<2} tracks: {}",
n_generation_depth, alive_secondary),
6);
}
simulation::shared_secondary_bank_read =
std::move(simulation::shared_secondary_bank_write);
simulation::shared_secondary_bank_write = SharedArray<SourceSite>();
simulation::progeny_per_particle.resize(
simulation::shared_secondary_bank_read.size());
std::fill(simulation::progeny_per_particle.begin(),
simulation::progeny_per_particle.end(), 0);
// Transport all secondary tracks from the shared secondary bank
#pragma omp parallel
{
vector<SourceSite> thread_bank;
#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);
transport_history_based_single_particle(p);
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);
}
}
} // End of transport loop over tracks in shared secondary bank
n_generation_depth++;
simulation::simulation_tracks_completed += alive_secondary;
} // End of loop over secondary generations
// Reset work so that fission bank etc works correctly
calculate_work(settings::n_particles);
}
void transport_event_based()
{
int64_t remaining_work = simulation::work_per_rank;
@ -848,33 +1033,7 @@ void transport_event_based()
// Initialize all particle histories for this subiteration
process_init_events(n_particles, source_offset);
// Event-based transport loop
while (true) {
// Determine which event kernel has the longest queue
int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
simulation::calculate_nonfuel_xs_queue.size(),
simulation::advance_particle_queue.size(),
simulation::surface_crossing_queue.size(),
simulation::collision_queue.size()});
// Execute event with the longest queue
if (max == 0) {
break;
} else if (max == simulation::calculate_fuel_xs_queue.size()) {
process_calculate_xs_events(simulation::calculate_fuel_xs_queue);
} else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
} else if (max == simulation::advance_particle_queue.size()) {
process_advance_particle_events();
} else if (max == simulation::surface_crossing_queue.size()) {
process_surface_crossing_events();
} else if (max == simulation::collision_queue.size()) {
process_collision_events();
}
}
// Execute death event for all particles
process_transport_events();
process_death_events(n_particles);
// Adjust remaining work and source offset variables
@ -883,4 +1042,122 @@ void transport_event_based()
}
}
void transport_event_based_shared_secondary()
{
// Clear shared secondary banks from any prior use
simulation::shared_secondary_bank_read.clear();
simulation::shared_secondary_bank_write.clear();
if (mpi::master) {
write_message(fmt::format(" Primary source particles: {}",
settings::n_particles),
6);
}
simulation::progeny_per_particle.resize(simulation::work_per_rank);
std::fill(simulation::progeny_per_particle.begin(),
simulation::progeny_per_particle.end(), 0);
// Phase 1: Transport primary particles using event-based processing and
// deposit first generation of secondaries in the shared secondary bank
int64_t remaining_work = simulation::work_per_rank;
int64_t source_offset = 0;
while (remaining_work > 0) {
int64_t n_particles =
std::min(remaining_work, settings::max_particles_in_flight);
process_init_events(n_particles, source_offset);
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();
}
remaining_work -= n_particles;
source_offset += n_particles;
}
simulation::simulation_tracks_completed += settings::n_particles;
// Phase 2: Now that the secondary bank has been populated, enter loop over
// all secondary generations
int n_generation_depth = 1;
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.
alive_secondary = synchronize_global_secondary_bank(
simulation::shared_secondary_bank_write);
// Recalculate work for each MPI rank based on number of alive secondary
// tracks
calculate_work(alive_secondary);
if (mpi::master) {
write_message(fmt::format(" Secondary generation {:<2} tracks: {}",
n_generation_depth, alive_secondary),
6);
}
simulation::shared_secondary_bank_read =
std::move(simulation::shared_secondary_bank_write);
simulation::shared_secondary_bank_write = SharedArray<SourceSite>();
simulation::progeny_per_particle.resize(
simulation::shared_secondary_bank_read.size());
std::fill(simulation::progeny_per_particle.begin(),
simulation::progeny_per_particle.end(), 0);
// Ensure particle buffer is large enough for this secondary generation
int64_t sec_buffer_length = std::min(
static_cast<int64_t>(simulation::shared_secondary_bank_read.size()),
settings::max_particles_in_flight);
if (sec_buffer_length >
static_cast<int64_t>(simulation::particles.size())) {
init_event_queues(sec_buffer_length);
}
// Transport secondary tracks using event-based processing
int64_t sec_remaining = simulation::shared_secondary_bank_read.size();
int64_t sec_offset = 0;
while (sec_remaining > 0) {
int64_t n_particles =
std::min(sec_remaining, settings::max_particles_in_flight);
process_init_secondary_events(
n_particles, sec_offset, simulation::shared_secondary_bank_read);
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();
}
sec_remaining -= n_particles;
sec_offset += n_particles;
} // End of subiteration loop over secondary tracks
n_generation_depth++;
simulation::simulation_tracks_completed += alive_secondary;
} // End of loop over secondary generations
// Reset work so that fission bank etc works correctly
calculate_work(settings::n_particles);
}
} // namespace openmc

30
src/source.cpp
View file

@ -354,6 +354,19 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node)
if (check_for_node(node, "energy")) {
pugi::xml_node node_dist = node.child("energy");
energy_ = distribution_from_xml(node_dist);
// For decay photon sources, use the absolute photon emission rate in
// [photons/s] as the source strength
if (dynamic_cast<DecaySpectrum*>(energy_.get())) {
if (strength_ != 1.0) {
warning(fmt::format(
"Source strength of {} is ignored because the source uses a "
"DecaySpectrum energy distribution. The source strength will be "
"set from the DecaySpectrum emission rate.",
strength_));
}
strength_ = energy_->integral();
}
} else {
// Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
energy_ = UPtrDist {new Watt(0.988e6, 2.249e-6)};
@ -414,6 +427,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
// Check for monoenergetic source above maximum particle energy
auto p = particle_.transport_index();
auto energy_ptr = dynamic_cast<Discrete*>(energy_.get());
auto decay_spectrum = dynamic_cast<DecaySpectrum*>(energy_.get());
if (energy_ptr) {
auto energies =
tensor::Tensor<double>(energy_ptr->x().data(), energy_ptr->x().size());
@ -424,10 +438,18 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
}
while (true) {
// Sample energy spectrum
auto [E, E_wgt_temp] = energy_->sample(seed);
site.E = E;
E_wgt = E_wgt_temp;
// Sample energy spectrum. For decay photon sources, also get the parent
// nuclide index to store in the source site for tallying purposes.
if (decay_spectrum) {
auto sample = decay_spectrum->sample_with_parent(seed);
site.E = sample.energy;
E_wgt = sample.weight;
site.parent_nuclide = sample.parent_nuclide;
} else {
auto [E, E_wgt_temp] = energy_->sample(seed);
site.E = E;
E_wgt = E_wgt_temp;
}
// Resample if energy falls above maximum particle energy
if (site.E < data::energy_max[p] &&

2
src/tallies/filter_particle_production.cpp
View file

@ -19,7 +19,7 @@ void ParticleProductionFilter::get_all_bins(
// Loop over secondary bank entries
for (int bank_idx = start_idx; bank_idx < end_idx; bank_idx++) {
const auto& site = p.secondary_bank(bank_idx);
const auto& site = p.local_secondary_bank(bank_idx);
// Find which particle-type slot this secondary belongs to
auto it = type_to_index_.find(site.particle.pdg_number());

20
src/tallies/tally.cpp
View file

@ -644,8 +644,24 @@ void Tally::set_scores(const vector<std::string>& scores)
break;
case HEATING:
if (settings::photon_transport)
estimator_ = TallyEstimator::COLLISION;
if (settings::photon_transport) {
// Photon heating requires a collision estimator (analog energy
// balance). However, if the tally only scores neutrons, we can keep the
// tracklength estimator since neutron heating uses kerma coefficients
// that support tracklength scoring.
bool neutron_only = false;
for (auto i_filt : filters_) {
auto pf =
dynamic_cast<ParticleFilter*>(model::tally_filters[i_filt].get());
if (pf && pf->particles().size() == 1 &&
pf->particles()[0].is_neutron()) {
neutron_only = true;
break;
}
}
if (!neutron_only)
estimator_ = TallyEstimator::COLLISION;
}
break;
case SCORE_PULSE_HEIGHT: {

9
src/tallies/tally_scoring.cpp
View file

@ -945,7 +945,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
if (p.type().is_neutron() && p.fission()) {
if (is_generation_time_or_both()) {
const auto& lifetimes =
simulation::ifp_source_lifetime_bank[p.current_work() - 1];
simulation::ifp_source_lifetime_bank[p.current_work()];
if (lifetimes.size() == settings::ifp_n_generation) {
score = lifetimes[0] * p.wgt_last();
}
@ -959,7 +959,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
if (p.type().is_neutron() && p.fission()) {
if (is_beta_effective_or_both()) {
const auto& delayed_groups =
simulation::ifp_source_delayed_group_bank[p.current_work() - 1];
simulation::ifp_source_delayed_group_bank[p.current_work()];
if (delayed_groups.size() == settings::ifp_n_generation) {
if (delayed_groups[0] > 0) {
score = p.wgt_last();
@ -985,12 +985,11 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
int ifp_data_size;
if (is_beta_effective_or_both()) {
ifp_data_size = static_cast<int>(
simulation::ifp_source_delayed_group_bank[p.current_work() - 1]
simulation::ifp_source_delayed_group_bank[p.current_work()]
.size());
} else {
ifp_data_size = static_cast<int>(
simulation::ifp_source_lifetime_bank[p.current_work() - 1]
.size());
simulation::ifp_source_lifetime_bank[p.current_work()].size());
}
if (ifp_data_size == settings::ifp_n_generation) {
score = p.wgt_last();

BIN
tests/regression_tests/collision_track/case_1_Reactions/collision_track_true.h5 Normal file
View file

Binary file not shown.

6
tests/regression_tests/collision_track/case_1_Reactions/inputs_true.dat
View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -51,7 +51,7 @@
</source>
<collision_track>
<reactions>(n,fission) 101</reactions>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

2
tests/regression_tests/collision_track/case_1_Reactions/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

BIN
tests/regression_tests/collision_track/case_2_Cell_ID/collision_track_true.h5 Normal file
View file

Binary file not shown.

6
tests/regression_tests/collision_track/case_2_Cell_ID/inputs_true.dat
View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -51,7 +51,7 @@
</source>
<collision_track>
<cell_ids>22</cell_ids>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

2
tests/regression_tests/collision_track/case_2_Cell_ID/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

BIN
tests/regression_tests/collision_track/case_3_Material_ID/collision_track_true.h5 Normal file
View file

Binary file not shown.

6
tests/regression_tests/collision_track/case_3_Material_ID/inputs_true.dat
View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -51,7 +51,7 @@
</source>
<collision_track>
<material_ids>1</material_ids>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

2
tests/regression_tests/collision_track/case_3_Material_ID/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

BIN
tests/regression_tests/collision_track/case_4_Nuclide_ID/collision_track_true.h5 Normal file
View file

Binary file not shown.

6
tests/regression_tests/collision_track/case_4_Nuclide_ID/inputs_true.dat
View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -51,7 +51,7 @@
</source>
<collision_track>
<nuclides>O16 U235</nuclides>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

2
tests/regression_tests/collision_track/case_4_Nuclide_ID/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

BIN
tests/regression_tests/collision_track/case_5_Universe_ID/collision_track_true.h5 Normal file
View file

Binary file not shown.

6
tests/regression_tests/collision_track/case_5_Universe_ID/inputs_true.dat
View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -52,7 +52,7 @@
<collision_track>
<cell_ids>22</cell_ids>
<universe_ids>77</universe_ids>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

2
tests/regression_tests/collision_track/case_5_Universe_ID/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

BIN
tests/regression_tests/collision_track/case_6_deposited_energy_threshold/collision_track_true.h5 Normal file
View file

Binary file not shown.

6
tests/regression_tests/collision_track/case_6_deposited_energy_threshold/inputs_true.dat
View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -51,7 +51,7 @@
</source>
<collision_track>
<deposited_E_threshold>550000.0</deposited_E_threshold>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

2
tests/regression_tests/collision_track/case_6_deposited_energy_threshold/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

BIN
tests/regression_tests/collision_track/case_7_all_parameters_used_together/collision_track_true.h5 Normal file
View file

Binary file not shown.

6
tests/regression_tests/collision_track/case_7_all_parameters_used_together/inputs_true.dat
View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -56,7 +56,7 @@
<material_ids>1 11</material_ids>
<nuclides>U238 U235 H1 U234</nuclides>
<deposited_E_threshold>100000.0</deposited_E_threshold>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

2
tests/regression_tests/collision_track/case_7_all_parameters_used_together/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

57
tests/regression_tests/collision_track/case_8_2threads/inputs_true.dat
View file

@ -1,57 +0,0 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="11.0" units="g/cm3"/>
<nuclide name="U234" ao="0.0004524"/>
<nuclide name="U235" ao="0.0506068"/>
<nuclide name="U238" ao="0.948709"/>
<nuclide name="U236" ao="0.0002318"/>
<nuclide name="O16" ao="2.0"/>
</material>
<material id="11">
<density value="1.0" units="g/cm3"/>
<nuclide name="H1" ao="2.0"/>
<nuclide name="O16" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="5" fill="77" region="-5 4 -7 6 -9 8" universe="1"/>
<cell id="8" material="11" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="1"/>
<cell id="22" material="1" region="-1 2 -3" universe="77"/>
<cell id="33" material="11" region="1 | -2 | 3" universe="77"/>
<surface id="1" type="z-cylinder" coeffs="0.0 0.0 2.0"/>
<surface id="2" type="z-plane" coeffs="-2.0"/>
<surface id="3" type="z-plane" coeffs="2.0"/>
<surface id="4" type="x-plane" coeffs="-3.0"/>
<surface id="5" type="x-plane" coeffs="3.0"/>
<surface id="6" type="y-plane" coeffs="-3.0"/>
<surface id="7" type="y-plane" coeffs="3.0"/>
<surface id="8" type="z-plane" coeffs="-3.0"/>
<surface id="9" type="z-plane" coeffs="3.0"/>
<surface id="10" type="x-plane" boundary="vacuum" coeffs="-4.0"/>
<surface id="11" type="x-plane" boundary="vacuum" coeffs="4.0"/>
<surface id="12" type="y-plane" boundary="vacuum" coeffs="-4.0"/>
<surface id="13" type="y-plane" boundary="vacuum" coeffs="4.0"/>
<surface id="14" type="z-plane" boundary="vacuum" coeffs="-4.0"/>
<surface id="15" type="z-plane" boundary="vacuum" coeffs="4.0"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>5</batches>
<inactive>1</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
</space>
<constraints>
<fissionable>true</fissionable>
</constraints>
</source>
<collision_track>
<max_collisions>200</max_collisions>
</collision_track>
<seed>1</seed>
</settings>
</model>

2
tests/regression_tests/collision_track/case_8_2threads/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

51
tests/regression_tests/collision_track/test.py
View file

@ -59,28 +59,13 @@ TODO:
"""
import os
import openmc
import openmc.lib
import pytest
from tests.testing_harness import CollisionTrackTestHarness
from tests.regression_tests import config
@pytest.fixture(scope="function")
def two_threads(monkeypatch):
"""Set the number of OMP threads to 2 for the test."""
monkeypatch.setenv("OMP_NUM_THREADS", "2")
@pytest.fixture(scope="function")
def single_process(monkeypatch):
"""Set the number of MPI process to 1 for the test."""
monkeypatch.setitem(config, "mpi_np", "1")
@pytest.fixture(scope="module")
def model_1():
"""Cylindrical core contained in a first box which is contained in a larger box.
@ -181,9 +166,9 @@ def model_1():
# =============================================================================
model.settings = openmc.Settings()
model.settings.particles = 100
model.settings.particles = 80
model.settings.batches = 5
model.settings.inactive = 1
model.settings.inactive = 4
model.settings.seed = 1
bounds = [
@ -203,19 +188,19 @@ def model_1():
@pytest.mark.parametrize(
"folder, model_name, parameter",
[("case_1_Reactions", "model_1", {"max_collisions": 300, "reactions": ["(n,fission)", 101]}),
[("case_1_Reactions", "model_1", {"max_collisions": 100, "reactions": ["(n,fission)", 101]}),
("case_2_Cell_ID", "model_1", {
"max_collisions": 300, "cell_ids": [22]}),
"max_collisions": 100, "cell_ids": [22]}),
("case_3_Material_ID", "model_1", {
"max_collisions": 300, "material_ids": [1]}),
"max_collisions": 100, "material_ids": [1]}),
("case_4_Nuclide_ID", "model_1", {
"max_collisions": 300, "nuclides": ["O16", "U235"]}),
"max_collisions": 100, "nuclides": ["O16", "U235"]}),
("case_5_Universe_ID", "model_1", {
"max_collisions": 300, "cell_ids": [22], "universe_ids": [77]}),
"max_collisions": 100, "cell_ids": [22], "universe_ids": [77]}),
("case_6_deposited_energy_threshold", "model_1", {
"max_collisions": 300, "deposited_E_threshold": 5.5e5}),
"max_collisions": 100, "deposited_E_threshold": 5.5e5}),
("case_7_all_parameters_used_together", "model_1", {
"max_collisions": 300,
"max_collisions": 100,
"reactions": ["elastic", 18, "(n,disappear)"],
"material_ids": [1, 11],
"universe_ids": [77],
@ -235,21 +220,3 @@ def test_collision_track_several_cases(
"statepoint.5.h5", model=model, workdir=folder
)
harness.main()
@pytest.mark.skipif(config["event"], reason="Results from history-based mode.")
def test_collision_track_2threads(model_1, two_threads, single_process):
# This test checks that the `max_collisions` setting is honored:
# no collisions beyond the specified limit should be recorded.
#
# For the result to be reproducible, the number of threads and
# the transport mode (history vs. event) must remain fixed.
assert os.environ["OMP_NUM_THREADS"] == "2"
assert config["mpi_np"] == "1"
model_1.settings.collision_track = {
"max_collisions": 200
}
harness = CollisionTrackTestHarness(
"statepoint.5.h5", model=model_1, workdir="case_8_2threads"
)
harness.main()

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

Some files were not shown because too many files have changed in this diff Show more