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

Merge branch 'openmc-dev:develop' into depletion-fy-energy-interp

Browse source
This commit is contained in:
YimengChan0 2026-05-25 19:59:25 +02:00 • committed by GitHub
commit 11cfdc9abd
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
93 changed files with 4231 additions and 793 deletions
Download patch file Download diff file Expand all files Collapse all files

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

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

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

@ -741,14 +741,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 +1060,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 +1090,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 +1123,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 +1158,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 +1207,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 +1332,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

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

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`,

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

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

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

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

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)

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

5
openmc/lib/core.py
View file

@ -31,7 +31,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

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,

116
openmc/mesh.py
View file

@ -268,10 +268,11 @@ class MeshBase(IDManagerMixin, ABC):
return string
def _volume_dim_check(self):
if self.n_dimension != 3 or \
any([d == 0 for d in self.dimension]):
raise RuntimeError(f'Mesh {self.id} is not 3D. '
'Volumes cannot be provided.')
if any(d == 0 for d in self.dimension):
raise RuntimeError(
f'Mesh {self.id} has a zero-size dimension. '
'Volumes cannot be provided.'
)
@classmethod
def from_hdf5(cls, group: h5py.Group):
@ -2737,7 +2738,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,
@ -3048,7 +3050,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:
@ -3105,35 +3109,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:
@ -3155,42 +3166,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):

2
openmc/model/model.py
View file

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

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

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

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

14
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()++;
}
@ -1225,7 +1229,7 @@ void sample_secondary_photons(Particle& p, int i_nuclide)
// 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()++;
}

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

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

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

65
tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/inputs_true.dat Normal file
View file

@ -0,0 +1,65 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>2g.h5</cross_sections>
<material id="1" name="mat_1">
<density value="1.0" units="macro"/>
<macroscopic name="mat_1"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="1 -2" universe="0"/>
<surface id="1" type="x-plane" boundary="reflective" coeffs="0.0"/>
<surface id="2" type="x-plane" boundary="vacuum" coeffs="929.45"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>2</batches>
<inactive>0</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>0.0 -1000.0 -1000.0 929.45 1000.0 1000.0</parameters>
</space>
</source>
<output>
<summary>false</summary>
</output>
<energy_mode>multi-group</energy_mode>
<tabular_legendre>
<enable>false</enable>
</tabular_legendre>
<create_fission_neutrons>true</create_fission_neutrons>
<weight_windows id="1">
<mesh>1</mesh>
<particle_type>neutron</particle_type>
<energy_bounds>0.0 0.625 20000000.0</energy_bounds>
<lower_ww_bounds>0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5</lower_ww_bounds>
<upper_ww_bounds>2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5</upper_ww_bounds>
<survival_ratio>3.0</survival_ratio>
<max_split>10</max_split>
<weight_cutoff>1e-38</weight_cutoff>
</weight_windows>
<mesh id="1">
<dimension>5 1 1</dimension>
<lower_left>0.0 -1000.0 -1000.0</lower_left>
<upper_right>929.45 1000.0 1000.0</upper_right>
</mesh>
<shared_secondary_bank>true</shared_secondary_bank>
<max_history_splits>100</max_history_splits>
</settings>
<tallies>
<mesh id="2">
<dimension>5 1 1</dimension>
<lower_left>0.0 -1000.0 -1000.0</lower_left>
<upper_right>929.45 1000.0 1000.0</upper_right>
</mesh>
<filter id="1" type="mesh">
<bins>2</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>flux</scores>
</tally>
</tallies>
</model>

11
tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/results_true.dat Normal file
View file

@ -0,0 +1,11 @@
tally 1:
1.765369E+04
3.087787E+08
1.708316E+04
2.842002E+08
9.444106E+03
7.488341E+07
2.066528E+03
2.142445E+06
8.689619E+02
5.652099E+05

91
tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/test.py Normal file
View file

@ -0,0 +1,91 @@
import os
import numpy as np
import openmc
from openmc.examples import slab_mg
from tests.testing_harness import PyAPITestHarness
def create_library():
# Instantiate the energy group data and file object
groups = openmc.mgxs.EnergyGroups([0.0, 0.625, 20.0e6])
mg_cross_sections_file = openmc.MGXSLibrary(groups)
# Make the base, isotropic data
nu = [2.50, 2.50]
fiss = np.array([0.002817, 0.097])
capture = [0.008708, 0.02518]
absorption = np.add(capture, fiss)
scatter = np.array(
[[[0.31980, 0.06694], [0.004555, -0.0003972]],
[[0.00000, 0.00000], [0.424100, 0.05439000]]])
total = [0.33588, 0.54628]
chi = [1., 0.]
mat_1 = openmc.XSdata('mat_1', groups)
mat_1.order = 1
mat_1.set_nu_fission(np.multiply(nu, fiss))
mat_1.set_absorption(absorption)
mat_1.set_scatter_matrix(scatter)
mat_1.set_total(total)
mat_1.set_chi(chi)
mg_cross_sections_file.add_xsdata(mat_1)
# Write the file
mg_cross_sections_file.export_to_hdf5('2g.h5')
class MGXSTestHarness(PyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = '2g.h5'
if os.path.exists(f):
os.remove(f)
def test_mg_fixed_source_ww_fission_shared_secondary():
create_library()
model = slab_mg()
# Override settings for fixed-source mode with shared secondary bank
model.settings.run_mode = 'fixed source'
model.settings.inactive = 0
model.settings.batches = 2
model.settings.particles = 100
model.settings.create_fission_neutrons = True
model.settings.shared_secondary_bank = True
model.settings.max_history_splits = 100
# Add weight windows on a simple 1D mesh
ww_mesh = openmc.RegularMesh()
ww_mesh.lower_left = (0.0, -1000.0, -1000.0)
ww_mesh.upper_right = (929.45, 1000.0, 1000.0)
ww_mesh.dimension = (5, 1, 1)
# Uniform lower bounds for 2 energy groups, 5 spatial bins
lower_bounds = np.full((2, 5, 1, 1), 0.5)
ww = openmc.WeightWindows(
ww_mesh,
lower_bounds.flatten(),
None,
5.0,
[0.0, 0.625, 20.0e6],
'neutron'
)
model.settings.weight_windows = [ww]
# Add a flux tally
mesh = openmc.RegularMesh()
mesh.lower_left = (0.0, -1000.0, -1000.0)
mesh.upper_right = (929.45, 1000.0, 1000.0)
mesh.dimension = (5, 1, 1)
tally = openmc.Tally()
tally.filters = [openmc.MeshFilter(mesh)]
tally.scores = ['flux']
model.tallies = [tally]
harness = MGXSTestHarness('statepoint.2.h5', model)
harness.main()

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

38
tests/regression_tests/particle_production_fission/local/inputs_true.dat Normal file
View file

@ -0,0 +1,38 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="18.0" units="g/cm3"/>
<nuclide name="U235" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-1" universe="1"/>
<surface id="1" type="sphere" boundary="vacuum" coeffs="0.0 0.0 0.0 5.0"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>2</batches>
<source type="independent" strength="1.0" particle="neutron">
<energy type="discrete">
<parameters>1000000.0 1.0</parameters>
</energy>
</source>
<create_fission_neutrons>true</create_fission_neutrons>
<shared_secondary_bank>false</shared_secondary_bank>
</settings>
<tallies>
<filter id="2" type="particle">
<bins>neutron</bins>
</filter>
<filter id="1" type="particleproduction">
<particles>neutron</particles>
</filter>
<tally id="1">
<filters>2 1</filters>
<scores>events</scores>
<estimator>analog</estimator>
</tally>
</tallies>
</model>

3
tests/regression_tests/particle_production_fission/local/results_true.dat Normal file
View file

@ -0,0 +1,3 @@
tally 1:
4.570000E+00
1.047890E+01

38
tests/regression_tests/particle_production_fission/shared/inputs_true.dat Normal file
View file

@ -0,0 +1,38 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="18.0" units="g/cm3"/>
<nuclide name="U235" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-1" universe="1"/>
<surface id="1" type="sphere" boundary="vacuum" coeffs="0.0 0.0 0.0 5.0"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>2</batches>
<source type="independent" strength="1.0" particle="neutron">
<energy type="discrete">
<parameters>1000000.0 1.0</parameters>
</energy>
</source>
<create_fission_neutrons>true</create_fission_neutrons>
<shared_secondary_bank>true</shared_secondary_bank>
</settings>
<tallies>
<filter id="2" type="particle">
<bins>neutron</bins>
</filter>
<filter id="1" type="particleproduction">
<particles>neutron</particles>
</filter>
<tally id="1">
<filters>2 1</filters>
<scores>events</scores>
<estimator>analog</estimator>
</tally>
</tallies>
</model>

3
tests/regression_tests/particle_production_fission/shared/results_true.dat Normal file
View file

@ -0,0 +1,3 @@
tally 1:
5.180000E+00
1.355140E+01

49
tests/regression_tests/particle_production_fission/test.py Normal file
View file

@ -0,0 +1,49 @@
import openmc
import pytest
from openmc.utility_funcs import change_directory
from tests.testing_harness import PyAPITestHarness
@pytest.mark.parametrize("shared_secondary,subdir", [
(False, "local"),
(True, "shared"),
])
def test_particle_production_fission(shared_secondary, subdir):
"""Fixed-source model with fissionable material to test that
ParticleProductionFilter correctly counts fission-born neutrons,
with both local and shared secondary bank modes."""
with change_directory(subdir):
openmc.reset_auto_ids()
model = openmc.Model()
mat = openmc.Material()
mat.set_density('g/cm3', 18.0)
mat.add_nuclide('U235', 1.0)
model.materials.append(mat)
sph = openmc.Sphere(r=5.0, boundary_type='vacuum')
cell = openmc.Cell(fill=mat, region=-sph)
model.geometry = openmc.Geometry([cell])
source = openmc.IndependentSource()
source.energy = openmc.stats.delta_function(1.0e6)
model.settings.particles = 100
model.settings.run_mode = 'fixed source'
model.settings.batches = 2
model.settings.source = source
model.settings.create_fission_neutrons = True
model.settings.shared_secondary_bank = shared_secondary
# ParticleProductionFilter tracking fission neutron production
ppf = openmc.ParticleProductionFilter(['neutron'])
neutron_filter = openmc.ParticleFilter(['neutron'])
tally = openmc.Tally()
tally.filters = [neutron_filter, ppf]
tally.scores = ['events']
tally.estimator = 'analog'
model.tallies = [tally]
harness = PyAPITestHarness('statepoint.2.h5', model)
harness.main()

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

9
tests/regression_tests/particle_restart_fixed_shared_secondary/geometry.xml Normal file
View file

@ -0,0 +1,9 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with no boundary condition -->
<surface id="1" type="sphere" coeffs="0.0 0.0 0.0 80.0" />
<surface id="2" type="sphere" coeffs="0.0 0.0 0.0 10000.0" boundary="vacuum" />
<cell id="1" material="1" region="-1" />
</geometry>

9
tests/regression_tests/particle_restart_fixed_shared_secondary/materials.xml Normal file
View file

@ -0,0 +1,9 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="1.4" units="g/cc" />
<nuclide name="B10" ao="1.0" />
</material>
</materials>

16
tests/regression_tests/particle_restart_fixed_shared_secondary/results_true.dat Normal file
View file

@ -0,0 +1,16 @@
current batch:
4.000000E+00
current generation:
1.000000E+00
particle id:
3.241000E+03
run mode:
fixed source
particle weight:
1.000000E+00
particle energy:
3.896365E+06
particle xyz:
8.710681E-01 3.698823E+00 -2.286229E+00
particle uvw:
-5.882735E-01 4.665422E-01 -6.605093E-01

15
tests/regression_tests/particle_restart_fixed_shared_secondary/settings.xml Normal file
View file

@ -0,0 +1,15 @@
<?xml version="1.0"?>
<settings>
<run_mode>fixed source</run_mode>
<batches>12</batches>
<particles>1000</particles>
<shared_secondary_bank>true</shared_secondary_bank>
<source>
<space type="box">
<parameters>-10 -10 -5 10 10 5</parameters>
</space>
</source>
</settings>

6
tests/regression_tests/particle_restart_fixed_shared_secondary/test.py Normal file
View file

@ -0,0 +1,6 @@
from tests.testing_harness import ParticleRestartTestHarness
def test_particle_restart_fixed_shared_secondary():
harness = ParticleRestartTestHarness('particle_4_3241.h5')
harness.main()

4
tests/regression_tests/pulse_height/inputs_true.dat → tests/regression_tests/pulse_height/local/inputs_true.dat
View file

@ -18,14 +18,12 @@
<particles>100</particles>
<batches>5</batches>
<source type="independent" strength="1.0" particle="photon">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
<energy type="discrete">
<parameters>1000000.0 1.0</parameters>
</energy>
</source>
<photon_transport>true</photon_transport>
<shared_secondary_bank>false</shared_secondary_bank>
</settings>
<tallies>
<filter id="1" type="cell">

0
tests/regression_tests/pulse_height/results_true.dat → tests/regression_tests/pulse_height/local/results_true.dat
View file

40
tests/regression_tests/pulse_height/shared/inputs_true.dat Normal file
View file

@ -0,0 +1,40 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1">
<density value="3.7" units="g/cc"/>
<nuclide name="Na23" ao="1.0"/>
<nuclide name="I127" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" name="inner sphere" material="1" region="-1" universe="1"/>
<cell id="2" name="outer sphere" material="void" region="1 -2" universe="1"/>
<surface id="1" type="sphere" coeffs="0.0 0.0 0.0 1"/>
<surface id="2" type="sphere" boundary="vacuum" coeffs="0.0 0.0 0.0 2"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>5</batches>
<source type="independent" strength="1.0" particle="photon">
<energy type="discrete">
<parameters>1000000.0 1.0</parameters>
</energy>
</source>
<photon_transport>true</photon_transport>
<shared_secondary_bank>true</shared_secondary_bank>
</settings>
<tallies>
<filter id="1" type="cell">
<bins>1</bins>
</filter>
<filter id="2" type="energy">
<bins>0.0 10000.0 20000.0 30000.0 40000.0 50000.0 60000.0 70000.0 80000.0 90000.0 100000.0 110000.0 120000.0 130000.0 140000.0 150000.0 160000.0 170000.0 180000.0 190000.0 200000.0 210000.0 220000.0 230000.0 240000.0 250000.0 260000.0 270000.0 280000.0 290000.0 300000.0 310000.0 320000.0 330000.0 340000.0 350000.0 360000.0 370000.0 380000.0 390000.0 400000.0 410000.0 420000.0 430000.0 440000.0 450000.0 460000.0 470000.0 480000.0 490000.0 500000.0 510000.0 520000.0 530000.0 540000.0 550000.0 560000.0 570000.0 580000.0 590000.0 600000.0 610000.0 620000.0 630000.0 640000.0 650000.0 660000.0 670000.0 680000.0 690000.0 700000.0 710000.0 720000.0 730000.0 740000.0 750000.0 760000.0 770000.0 780000.0 790000.0 800000.0 810000.0 820000.0 830000.0 840000.0 850000.0 860000.0 870000.0 880000.0 890000.0 900000.0 910000.0 920000.0 930000.0 940000.0 950000.0 960000.0 970000.0 980000.0 990000.0 1000000.0</bins>
</filter>
<tally id="1" name="pht tally">
<filters>1 2</filters>
<scores>pulse-height</scores>
</tally>
</tallies>
</model>

201
tests/regression_tests/pulse_height/shared/results_true.dat Normal file
View file

@ -0,0 +1,201 @@
tally 1:
4.140000E+00
3.443000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
3.000000E-02
5.000000E-04
2.000000E-02
4.000000E-04
2.000000E-02
4.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.000000E-02
3.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.000000E-02
6.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
2.000000E-02
4.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
3.000000E-02
5.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
2.000000E-02
2.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
2.000000E-02
2.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
3.000000E-02
5.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-01
8.600000E-03

80
tests/regression_tests/pulse_height/test.py
View file

@ -1,53 +1,53 @@
import numpy as np
import openmc
import pytest
from openmc.utility_funcs import change_directory
from tests.testing_harness import PyAPITestHarness
@pytest.fixture
def sphere_model():
model = openmc.model.Model()
@pytest.mark.parametrize("shared_secondary,subdir", [
(False, "local"),
(True, "shared"),
])
def test_pulse_height(shared_secondary, subdir):
with change_directory(subdir):
openmc.reset_auto_ids()
model = openmc.Model()
# Define materials
NaI = openmc.Material()
NaI.set_density('g/cc', 3.7)
NaI.add_element('Na', 1.0)
NaI.add_element('I', 1.0)
# Define materials
NaI = openmc.Material()
NaI.set_density('g/cc', 3.7)
NaI.add_element('Na', 1.0)
NaI.add_element('I', 1.0)
model.materials = openmc.Materials([NaI])
model.materials = openmc.Materials([NaI])
# Define geometry: two spheres in each other
s1 = openmc.Sphere(r=1)
s2 = openmc.Sphere(r=2, boundary_type='vacuum')
inner_sphere = openmc.Cell(name='inner sphere', fill=NaI, region=-s1)
outer_sphere = openmc.Cell(name='outer sphere', region=+s1 & -s2)
model.geometry = openmc.Geometry([inner_sphere, outer_sphere])
# Define geometry: two spheres in each other
s1 = openmc.Sphere(r=1)
s2 = openmc.Sphere(r=2, boundary_type='vacuum')
inner_sphere = openmc.Cell(name='inner sphere', fill=NaI, region=-s1)
outer_sphere = openmc.Cell(name='outer sphere', region=+s1 & -s2)
model.geometry = openmc.Geometry([inner_sphere, outer_sphere])
# Define settings
model.settings.run_mode = 'fixed source'
model.settings.batches = 5
model.settings.particles = 100
model.settings.photon_transport = True
model.settings.source = openmc.IndependentSource(
space=openmc.stats.Point(),
energy=openmc.stats.Discrete([1e6], [1]),
particle='photon'
)
# Define settings
model.settings.run_mode = 'fixed source'
model.settings.batches = 5
model.settings.particles = 100
model.settings.photon_transport = True
model.settings.shared_secondary_bank = shared_secondary
model.settings.source = openmc.IndependentSource(
energy=openmc.stats.delta_function(1e6),
particle='photon'
)
# Define tallies
tally = openmc.Tally(name="pht tally")
tally.scores = ['pulse-height']
cell_filter = openmc.CellFilter(inner_sphere)
energy_filter = openmc.EnergyFilter(np.linspace(0, 1_000_000, 101))
tally.filters = [cell_filter, energy_filter]
model.tallies = [tally]
# Define tallies
tally = openmc.Tally(name="pht tally")
tally.scores = ['pulse-height']
cell_filter = openmc.CellFilter(inner_sphere)
energy_filter = openmc.EnergyFilter(np.linspace(0, 1_000_000, 101))
tally.filters = [cell_filter, energy_filter]
model.tallies = [tally]
return model
def test_pulse_height(sphere_model):
harness = PyAPITestHarness('statepoint.5.h5', sphere_model)
harness.main()
harness = PyAPITestHarness('statepoint.5.h5', model)
harness.main()

94
tests/regression_tests/weightwindows/local/inputs_true.dat Normal file
View file

@ -0,0 +1,94 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1">
<density value="2.3" units="g/cc"/>
<nuclide name="H1" ao="0.168018676"/>
<nuclide name="O16" ao="0.561814465"/>
<nuclide name="O17" ao="0.00021401"/>
<nuclide name="Na23" ao="0.021365"/>
<nuclide name="Al27" ao="0.021343"/>
<nuclide name="Si28" ao="0.187439342"/>
<nuclide name="Si29" ao="0.009517714"/>
<nuclide name="Si30" ao="0.006273944"/>
<nuclide name="Ca40" ao="0.018026179"/>
<nuclide name="Ca42" ao="0.00012031"/>
<nuclide name="Ca44" ao="0.000387892"/>
<nuclide name="Fe54" ao="0.000248179"/>
<nuclide name="Fe56" ao="0.003895875"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-1" universe="1"/>
<cell id="2" material="void" region="1 -2" universe="1"/>
<surface id="1" type="sphere" coeffs="0.0 0.0 0.0 240"/>
<surface id="2" type="sphere" boundary="vacuum" coeffs="0.0 0.0 0.0 250"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>500</particles>
<batches>2</batches>
<source type="independent" strength="1.0" particle="neutron">
<space type="point">
<parameters>0.001 0.001 0.001</parameters>
</space>
<energy type="discrete">
<parameters>14000000.0 1.0</parameters>
</energy>
</source>
<photon_transport>true</photon_transport>
<weight_windows id="1">
<mesh>2</mesh>
<particle_type>neutron</particle_type>
<energy_bounds>0.0 0.5 20000000.0</energy_bounds>
<lower_ww_bounds>-1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.0008758251780046591 -1.0 -1.0 -1.0 -1.0 0.00044494017319042853 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 1.6620271125067025e-05 0.0006278777700923334 3.3816651814344154e-05 -1.0 -1.0 0.0024363004681119066 0.04404124277352227 0.002389900091734746 -1.0 -1.0 0.001096572213298872 0.04862206884590391 0.0011530054432113332 -1.0 -1.0 -1.0 0.00029204950777272335 2.2332845991385424e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.0001298973206651331 0.0028826281529980655 0.00018941326535850932 -1.0 5.4657521927731274e-05 0.014670839729146354 0.49999999999999994 0.011271060592765664 -1.0 -1.0 0.014846491082523524 0.4897211700090108 0.013284874451810723 -1.0 -1.0 5.14657989105535e-05 0.0010115278438204965 0.0008429411802845685 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.0009344129479768359 -1.0 -1.0 -1.0 0.0024828273390431962 0.04299740304329489 0.0020539079252555113 -1.0 -1.0 0.003034165905819096 0.04870927636605937 0.00313101668086297 -1.0 -1.0 -1.0 6.339910999436737e-05 7.442176086066386e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.00011276372995589871 0.0002236100157024887 9.56304298913265e-08 -1.0 -1.0 0.0001596955110781239 9.960576598335084e-06 5.3833030623153676e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 4.0453018186074215e-05 4.6013290110121894e-05 2.709738801641897e-05 -1.0 -1.0 -1.0 9.538410811938703e-05 1.992978141244964e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 1.0723231851298364e-05 -1.0 -1.0 -1.0 -1.0 0.0008030100296698905 0.017386220476187222 0.0005027758935317528 -1.0 -1.0 0.00103568411326969 0.015609071124009234 0.0007473731339616588 -1.0 -1.0 -1.0 6.979537549963374e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.0014257756927429828 5.9098070539559466e-05 -1.0 5.0168071459366625e-06 0.005588863190166347 0.5 0.003930588100414563 -1.0 -1.0 0.005163345217476071 0.48250750993887326 0.003510432129522241 -1.0 -1.0 3.371977841720992e-05 0.0006640103276501794 9.554899988419713e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.000187872910697387 -1.0 -1.0 -1.0 0.0009134587866163143 0.017081858078684467 0.0002972747357542354 -1.0 -1.0 0.0007973472495507653 0.019300903254809747 0.000902203032280694 -1.0 -1.0 6.170015233787292e-05 1.898984300674969e-05 9.85411583595722e-07 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 1.0173913765490095e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0</lower_ww_bounds>
<upper_ww_bounds>-10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.008758251780046591 -10.0 -10.0 -10.0 -10.0 0.004449401731904285 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00016620271125067024 0.006278777700923334 0.00033816651814344155 -10.0 -10.0 0.024363004681119065 0.4404124277352227 0.02389900091734746 -10.0 -10.0 0.01096572213298872 0.4862206884590391 0.011530054432113333 -10.0 -10.0 -10.0 0.0029204950777272335 0.00022332845991385425 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.0012989732066513312 0.028826281529980655 0.0018941326535850931 -10.0 0.0005465752192773128 0.14670839729146354 4.999999999999999 0.11271060592765664 -10.0 -10.0 0.14846491082523525 4.897211700090108 0.13284874451810724 -10.0 -10.0 0.000514657989105535 0.010115278438204964 0.008429411802845685 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00934412947976836 -10.0 -10.0 -10.0 0.024828273390431962 0.4299740304329489 0.020539079252555114 -10.0 -10.0 0.03034165905819096 0.48709276366059373 0.0313101668086297 -10.0 -10.0 -10.0 0.0006339910999436737 0.0007442176086066385 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.0011276372995589871 0.002236100157024887 9.56304298913265e-07 -10.0 -10.0 0.001596955110781239 9.960576598335084e-05 0.0005383303062315367 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00040453018186074213 0.00046013290110121896 0.0002709738801641897 -10.0 -10.0 -10.0 0.0009538410811938703 0.00019929781412449641 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00010723231851298364 -10.0 -10.0 -10.0 -10.0 0.008030100296698905 0.17386220476187222 0.0050277589353175285 -10.0 -10.0 0.010356841132696899 0.15609071124009233 0.007473731339616587 -10.0 -10.0 -10.0 0.0006979537549963374 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.014257756927429827 0.0005909807053955946 -10.0 5.016807145936663e-05 0.055888631901663474 5.0 0.039305881004145636 -10.0 -10.0 0.05163345217476071 4.825075099388733 0.03510432129522241 -10.0 -10.0 0.00033719778417209923 0.006640103276501793 0.0009554899988419713 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00187872910697387 -10.0 -10.0 -10.0 0.009134587866163143 0.17081858078684467 0.002972747357542354 -10.0 -10.0 0.007973472495507653 0.19300903254809748 0.009022030322806941 -10.0 -10.0 0.0006170015233787292 0.0001898984300674969 9.85411583595722e-06 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00010173913765490096 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0</upper_ww_bounds>
<survival_ratio>3.0</survival_ratio>
<max_lower_bound_ratio>1.5</max_lower_bound_ratio>
<max_split>10</max_split>
<weight_cutoff>1e-38</weight_cutoff>
</weight_windows>
<mesh id="2">
<dimension>5 6 7</dimension>
<lower_left>-240 -240 -240</lower_left>
<upper_right>240 240 240</upper_right>
</mesh>
<weight_windows id="2">
<mesh>2</mesh>
<particle_type>photon</particle_type>
<energy_bounds>0.0 0.5 20000000.0</energy_bounds>
<lower_ww_bounds>-1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 2.0445691096074744e-05 -1.0 -1.0 -1.0 0.00013281328509288383 0.0006267128082339883 -1.0 -1.0 -1.0 -1.0 0.0004209808495056742 5.340221214300681e-05 -1.0 -1.0 -1.0 1.553693492042344e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 3.7096847855980484e-05 -1.0 -1.0 5.41969614042385e-05 0.0007918090215756442 9.763572147337743e-05 -1.0 -1.0 0.0026398319229115233 0.04039871468258457 0.002806654735003014 -1.0 -1.0 0.0027608366482897244 0.040190978087723005 0.0025084416103979975 -1.0 -1.0 9.774868433298412e-05 0.0006126404916879651 0.00014841730774317252 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 9.75342408717107e-06 1.7610849405640072e-05 -1.0 -1.0 0.0001160370359598608 0.002822592083222468 0.0003897177301239376 -1.0 -1.0 0.014965179733172532 0.5 0.011723491704290401 -1.0 2.260493623875525e-05 0.015157364795884922 0.49376993953402387 0.013111419473204967 -1.0 -1.0 0.00014919915366377564 0.002197964829133137 0.0003209711829219673 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 1.4806422867046436e-05 -1.0 -1.0 -1.0 9.207268175745281e-05 0.0007066033561638983 -1.0 -1.0 -1.0 0.003246064903858183 0.03905493028065908 0.0025380574501073605 -1.0 -1.0 0.0032499260211917933 0.04335567738779479 0.0031577214557017056 4.521352809838806e-05 -1.0 0.00018268090756001903 0.0004353654810741932 0.00011273259573092372 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 9.650194311662424e-06 0.00029403373970835075 0.00013506203419326438 -1.0 -1.0 4.5065302725431816e-06 0.00027725323638744237 3.547290864255937e-05 -1.0 -1.0 -1.0 8.786172732576295e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0</lower_ww_bounds>
<upper_ww_bounds>-10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00020445691096074745 -10.0 -10.0 -10.0 0.0013281328509288383 0.006267128082339883 -10.0 -10.0 -10.0 -10.0 0.004209808495056742 0.0005340221214300681 -10.0 -10.0 -10.0 0.00015536934920423439 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00037096847855980485 -10.0 -10.0 0.000541969614042385 0.007918090215756443 0.0009763572147337743 -10.0 -10.0 0.026398319229115234 0.4039871468258457 0.02806654735003014 -10.0 -10.0 0.027608366482897245 0.40190978087723006 0.025084416103979976 -10.0 -10.0 0.0009774868433298411 0.0061264049168796506 0.0014841730774317252 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 9.753424087171071e-05 0.00017610849405640073 -10.0 -10.0 0.001160370359598608 0.02822592083222468 0.003897177301239376 -10.0 -10.0 0.14965179733172532 5.0 0.11723491704290401 -10.0 0.0002260493623875525 0.15157364795884923 4.937699395340239 0.13111419473204966 -10.0 -10.0 0.0014919915366377564 0.02197964829133137 0.003209711829219673 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00014806422867046437 -10.0 -10.0 -10.0 0.0009207268175745281 0.007066033561638983 -10.0 -10.0 -10.0 0.03246064903858183 0.39054930280659084 0.025380574501073606 -10.0 -10.0 0.03249926021191793 0.4335567738779479 0.03157721455701706 0.0004521352809838806 -10.0 0.0018268090756001904 0.004353654810741932 0.001127325957309237 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 9.650194311662424e-05 0.0029403373970835075 0.0013506203419326437 -10.0 -10.0 4.5065302725431816e-05 0.002772532363874424 0.0003547290864255937 -10.0 -10.0 -10.0 0.0008786172732576295 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0</upper_ww_bounds>
<survival_ratio>3.0</survival_ratio>
<max_lower_bound_ratio>1.5</max_lower_bound_ratio>
<max_split>10</max_split>
<weight_cutoff>1e-38</weight_cutoff>
</weight_windows>
<shared_secondary_bank>false</shared_secondary_bank>
<weight_window_checkpoints>
<collision>true</collision>
<surface>true</surface>
</weight_window_checkpoints>
<max_history_splits>200</max_history_splits>
</settings>
<tallies>
<mesh id="1">
<dimension>5 10 15</dimension>
<lower_left>-240 -240 -240</lower_left>
<upper_right>240 240 240</upper_right>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<filter id="2" type="energy">
<bins>0.0 0.5 20000000.0</bins>
</filter>
<filter id="3" type="particle">
<bins>neutron photon</bins>
</filter>
<tally id="1">
<filters>1 2 3</filters>
<scores>flux</scores>
</tally>
</tallies>
</model>

1
tests/regression_tests/weightwindows/local/results_true.dat Normal file
View file

@ -0,0 +1 @@
a78972fe3c0dadfc256cfb139c5ca8c7b634738e1895ad2d6286ed5e2567c6dc518e499363908e9058432684148a706a179a39110f703d5322491184a1d0c3e4

1
tests/regression_tests/weightwindows/results_true.dat
View file

@ -1 +0,0 @@
a4a3ccca43666e2ca1e71800201b152cca20c387b93d67522c5339807348dcee5cada9acbed3238f37e2e86e76b374b06988742f07d4ea1b413e4e75d0c180b1

1
tests/regression_tests/weightwindows/inputs_true.dat → tests/regression_tests/weightwindows/shared/inputs_true.dat
View file

@ -64,6 +64,7 @@
<max_split>10</max_split>
<weight_cutoff>1e-38</weight_cutoff>
</weight_windows>
<shared_secondary_bank>true</shared_secondary_bank>
<weight_window_checkpoints>
<collision>true</collision>
<surface>true</surface>

1
tests/regression_tests/weightwindows/shared/results_true.dat Normal file
View file

@ -0,0 +1 @@
27c2d218ecf46161088003fbb39ccd63495f80d6e0b27dbc3b3e845b49b282dabf0eaba4fb619be15c6116fc963ec5ce1fe4197efbe4084972ec761e02e70eb5

69
tests/regression_tests/weightwindows/survival_biasing/local/inputs_true.dat Normal file
View file

@ -0,0 +1,69 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" name="Tungsten">
<density value="19.25" units="g/cm3"/>
<nuclide name="W180" ao="0.0012"/>
<nuclide name="W182" ao="0.265"/>
<nuclide name="W183" ao="0.1431"/>
<nuclide name="W184" ao="0.3064"/>
<nuclide name="W186" ao="0.2843"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="1 -2 3 -4 5 -6" universe="1"/>
<surface id="1" type="x-plane" boundary="reflective" coeffs="0.0"/>
<surface id="2" type="x-plane" boundary="vacuum" coeffs="160.0"/>
<surface id="3" type="y-plane" boundary="reflective" coeffs="0.0"/>
<surface id="4" type="y-plane" boundary="reflective" coeffs="160.0"/>
<surface id="5" type="z-plane" boundary="reflective" coeffs="0.0"/>
<surface id="6" type="z-plane" boundary="reflective" coeffs="160.0"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>50</particles>
<batches>5</batches>
<source type="independent" strength="1.0" particle="neutron">
<space type="cartesian">
<x type="discrete">
<parameters>0.01 1.0</parameters>
</x>
<y type="uniform" parameters="0.0 160.0"/>
<z type="uniform" parameters="0.0 160.0"/>
</space>
<angle type="monodirectional" reference_uvw="1.0 0.0 0.0"/>
<energy type="discrete">
<parameters>14100000.0 1.0</parameters>
</energy>
</source>
<survival_biasing>true</survival_biasing>
<weight_windows id="1">
<mesh>1</mesh>
<particle_type>neutron</particle_type>
<lower_ww_bounds>0.46135961568957096 0.2874485390202603 0.13256013950494605 0.052765209609807295 0.020958440832685638 0.006317250035749876 0.002627037175682506 0.00030032343592437284 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4611178493390162 0.28624862560842024 0.13999870622621136 0.05508479408959756 0.018281241958254993 0.0055577764872361555 0.0015265432226293397 0.00024298772352636966 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.46294957913063317 0.2857734367546051 0.13365623190336945 0.05034742166227103 0.017525851812913266 0.005096725451851077 0.0026297640951531403 0.00033732424392026144 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.45772598750472554 0.281494921471495 0.13604397962762246 0.054792881472295364 0.019802376898755525 0.0073351745579128495 0.002067392946066426 9.529474085926143e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4739762983490318 0.28419061537046764 0.12757509901507888 0.0516920293019733 0.01919167874787235 0.007593035964707848 0.0017488176314107277 9.678267909681988e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.47560199098558276 0.27728389146956994 0.12760802572535623 0.05251965811713552 0.022498960644951632 0.01063372459325151 0.004242071054914633 0.00045603112202665183 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.47659199471934693 0.28886666944252215 0.1258637573398161 0.05818862375955657 0.020476313720744762 0.007143193244018263 0.0022364083022515273 0.0003953123781408474 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4693487369606823 0.27180862139490125 0.12833455570423483 0.055146743559938344 0.020667403529470333 0.007891508723137146 0.001643551705407358 0.0006501416781353278 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4538724931023337 0.2824491421084296 0.1275341706351566 0.055027501141893705 0.02305114640508087 0.008599303158607198 0.0025082611194161804 0.000518530311231696 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.45258444753220123 0.26839102466484255 0.12991633918797083 0.05442263709947767 0.019695895223471486 0.005353204961887518 0.0015074698293840157 0.0001598680898180058 8.872047880811405e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4482467985982804 0.2813276470183359 0.13449860790768647 0.056103278190533436 0.026274320334196962 0.006884071908429303 0.0033099390738735458 0.0005337454397674922 0.0001386380399465575 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4875748950139487 0.2908586705316248 0.14212060658392278 0.06068108009637272 0.019467970468789477 0.00637102427946399 0.0028510425266191687 0.0012184421778115488 0.00029404541567146187 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.5 0.2852101551297824 0.13622307172321643 0.058302519884339946 0.023124734915152625 0.007401527839638796 0.002527635652743992 0.0008509612315162482 0.00018213334574554768 7.815977984795461e-06 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.47706397688984176 0.28313436859495966 0.13170624744221976 0.0525453512766549 0.021111753158319105 0.0067451713955609645 0.003204270539718037 0.0010308107242263192 0.0002759210144794204 0.00013727805365387318 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.48641988659284513 0.2816173231471242 0.131541850061199 0.054793200248327956 0.016517701691156576 0.006757591781856257 0.0025425350750908644 0.0006383278085839443 2.279421641064226e-05 0.0003597432472224371 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4307755159245372 0.2652713784080849 0.12753857957535653 0.05470396852881577 0.02049817309420563 0.00566559741247352 0.0007213846765465147 6.848024591311009e-05 -1.0 8.609308814924014e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.43631324024956025 0.2639077825530567 0.13368737962742414 0.05213531603720763 0.020628860469743833 0.008268384844678654 0.0028599221013934375 0.00013059757129982714 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.45926287238219 0.2635824148883888 0.13104280996799478 0.055379812186041905 0.019255042561941074 0.007252095690246872 0.0018026593488140996 0.0003215836458542421 9.14840176000331e-06 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4574378998400604 0.28714971179927723 0.1363507911051337 0.05010315954654347 0.017963278989571074 0.0058998176297971605 0.0021393135666168 0.00015677575033083482 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4821442782978812 0.27704132233083506 0.14149074035662307 0.05546042815834048 0.016800649382269998 0.004096804603054233 0.0019161108398578026 0.0005289600253155307 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0</lower_ww_bounds>
<upper_ww_bounds>2.306798078447855 1.4372426951013015 0.6628006975247303 0.2638260480490365 0.10479220416342819 0.03158625017874938 0.013135185878412529 0.0015016171796218643 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.3055892466950807 1.4312431280421012 0.6999935311310568 0.2754239704479878 0.09140620979127496 0.027788882436180776 0.007632716113146698 0.0012149386176318483 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.314747895653166 1.4288671837730256 0.6682811595168472 0.25173710831135515 0.08762925906456634 0.025483627259255386 0.013148820475765701 0.0016866212196013073 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.288629937523628 1.407474607357475 0.6802198981381123 0.2739644073614768 0.09901188449377762 0.036675872789564246 0.01033696473033213 0.00047647370429630714 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.369881491745159 1.4209530768523382 0.6378754950753944 0.2584601465098665 0.09595839373936176 0.03796517982353924 0.008744088157053638 0.00048391339548409945 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.3780099549279137 1.3864194573478497 0.6380401286267812 0.26259829058567763 0.11249480322475816 0.053168622966257545 0.021210355274573163 0.0022801556101332593 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.382959973596735 1.4443333472126108 0.6293187866990806 0.29094311879778284 0.10238156860372381 0.035715966220091315 0.011182041511257637 0.001976561890704237 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.3467436848034113 1.3590431069745064 0.6416727785211741 0.2757337177996917 0.10333701764735166 0.03945754361568573 0.00821775852703679 0.0032507083906766388 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.2693624655116684 1.412245710542148 0.637670853175783 0.27513750570946854 0.11525573202540434 0.04299651579303599 0.012541305597080901 0.00259265155615848 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.262922237661006 1.3419551233242126 0.6495816959398542 0.2721131854973884 0.09847947611735743 0.02676602480943759 0.007537349146920078 0.000799340449090029 0.0004436023940405703 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.2412339929914022 1.4066382350916795 0.6724930395384323 0.28051639095266717 0.1313716016709848 0.03442035954214652 0.016549695369367727 0.0026687271988374613 0.0006931901997327875 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.4378744750697434 1.4542933526581239 0.7106030329196139 0.3034054004818636 0.09733985234394739 0.03185512139731995 0.014255212633095843 0.006092210889057744 0.0014702270783573093 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.5 1.4260507756489118 0.6811153586160822 0.2915125994216997 0.11562367457576313 0.03700763919819398 0.012638178263719959 0.004254806157581241 0.0009106667287277384 3.9079889923977307e-05 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.385319884449209 1.4156718429747983 0.6585312372110987 0.2627267563832745 0.10555876579159552 0.03372585697780482 0.016021352698590185 0.005154053621131596 0.001379605072397102 0.000686390268269366 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.432099432964226 1.4080866157356209 0.657709250305995 0.2739660012416398 0.08258850845578287 0.03378795890928128 0.012712675375454322 0.0031916390429197216 0.0001139710820532113 0.0017987162361121855 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.1538775796226863 1.3263568920404245 0.6376928978767826 0.27351984264407886 0.10249086547102815 0.028327987062367603 0.0036069233827325737 0.0003424012295655504 -5.0 0.0004304654407462007 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.181566201247801 1.3195389127652835 0.6684368981371207 0.2606765801860381 0.10314430234871916 0.041341924223393264 0.014299610506967188 0.0006529878564991358 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.29631436191095 1.317912074441944 0.6552140498399739 0.27689906093020955 0.09627521280970537 0.03626047845123436 0.009013296744070498 0.0016079182292712106 4.574200880001655e-05 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.287189499200302 1.435748558996386 0.6817539555256685 0.25051579773271737 0.08981639494785537 0.029499088148985803 0.010696567833083998 0.0007838787516541741 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.410721391489406 1.3852066116541752 0.7074537017831153 0.2773021407917024 0.08400324691134999 0.020484023015271167 0.009580554199289014 0.0026448001265776534 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0</upper_ww_bounds>
<survival_ratio>3.0</survival_ratio>
<max_split>10</max_split>
<weight_cutoff>1e-38</weight_cutoff>
</weight_windows>
<mesh id="1">
<dimension>20 20 1</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>160.0 160.0 160.0</upper_right>
</mesh>
<shared_secondary_bank>false</shared_secondary_bank>
<weight_window_checkpoints>
<collision>true</collision>
<surface>true</surface>
</weight_window_checkpoints>
</settings>
<tallies>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<tally id="1" name="flux">
<filters>1</filters>
<scores>flux</scores>
</tally>
</tallies>
</model>

0
tests/regression_tests/weightwindows/survival_biasing/results_true.dat → tests/regression_tests/weightwindows/survival_biasing/local/results_true.dat
View file

1
tests/regression_tests/weightwindows/survival_biasing/inputs_true.dat → tests/regression_tests/weightwindows/survival_biasing/shared/inputs_true.dat
View file

@ -51,6 +51,7 @@
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>160.0 160.0 160.0</upper_right>
</mesh>
<shared_secondary_bank>true</shared_secondary_bank>
<weight_window_checkpoints>
<collision>true</collision>
<surface>true</surface>

1
tests/regression_tests/weightwindows/survival_biasing/shared/results_true.dat Normal file
View file

@ -0,0 +1 @@
d17a437262d3316985fba4b48e21a7fcd5f32ac2d96ef0e66849f567deaeadc1e0f18d5d0bf5e9cab4246dbf823e7f43f249a1f67e928b27ae8c70b89f3cefbf

28
tests/regression_tests/weightwindows/survival_biasing/test.py
View file

@ -1,14 +1,16 @@
from pathlib import Path
import pytest
import numpy as np
import openmc
from openmc.stats import Discrete, Point
from openmc.utility_funcs import change_directory
from tests.testing_harness import HashedPyAPITestHarness
@pytest.fixture
def model():
def build_model(shared_secondary):
openmc.reset_auto_ids()
# Material
w = openmc.Material(name='Tungsten')
w.add_element('W', 1.0)
@ -38,13 +40,14 @@ def model():
angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0))
energy = openmc.stats.Discrete([14.1e6], [1.0])
source = openmc.Source(space=space, angle=angle, energy=energy)
source = openmc.IndependentSource(space=space, angle=angle, energy=energy)
settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.batches = 5
settings.particles = 50
settings.source = source
settings.shared_secondary_bank = shared_secondary
model = openmc.Model(geometry=geometry, materials=materials, settings=settings)
@ -61,7 +64,8 @@ def model():
tallies = openmc.Tallies([flux_tally])
model.tallies = tallies
lower_ww_bounds = np.loadtxt('ww_n.txt')
parent_dir = Path(__file__).parent
lower_ww_bounds = np.loadtxt(parent_dir / 'ww_n.txt')
weight_windows = openmc.WeightWindows(mesh,
lower_ww_bounds,
@ -76,6 +80,12 @@ def model():
return model
def test_weight_windows_with_survival_biasing(model):
harness = HashedPyAPITestHarness('statepoint.5.h5', model)
harness.main()
@pytest.mark.parametrize("shared_secondary,subdir", [
(False, "local"),
(True, "shared"),
])
def test_weight_windows_with_survival_biasing(shared_secondary, subdir):
with change_directory(subdir):
model = build_model(shared_secondary)
harness = HashedPyAPITestHarness('statepoint.5.h5', model)
harness.main()

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

@ -1,14 +1,17 @@
from pathlib import Path
import pytest
import numpy as np
import openmc
from openmc.stats import Discrete, Point
from openmc.utility_funcs import change_directory
from tests.testing_harness import HashedPyAPITestHarness
@pytest.fixture
def model():
def build_model(shared_secondary):
openmc.reset_auto_ids()
model = openmc.Model()
# materials (M4 steel alloy)
@ -43,6 +46,7 @@ def model():
settings.batches = 2
settings.max_history_splits = 200
settings.photon_transport = True
settings.shared_secondary_bank = shared_secondary
settings.weight_window_checkpoints = {'surface': True,
'collision': True}
space = Point((0.001, 0.001, 0.001))
@ -71,10 +75,10 @@ def model():
# weight windows
# load pre-generated weight windows
# (created using the same tally as above)
ww_n_lower_bnds = np.loadtxt('ww_n.txt')
ww_p_lower_bnds = np.loadtxt('ww_p.txt')
# load pre-generated weight windows from parent directory
parent_dir = Path(__file__).parent
ww_n_lower_bnds = np.loadtxt(parent_dir / 'ww_n.txt')
ww_p_lower_bnds = np.loadtxt(parent_dir / 'ww_p.txt')
# create a mesh matching the one used
# to generate the weight windows
@ -104,9 +108,15 @@ def model():
return model
def test_weightwindows(model):
test = HashedPyAPITestHarness('statepoint.2.h5', model)
test.main()
@pytest.mark.parametrize("shared_secondary,subdir", [
(False, "local"),
(True, "shared"),
])
def test_weightwindows(shared_secondary, subdir):
with change_directory(subdir):
model = build_model(shared_secondary)
test = HashedPyAPITestHarness('statepoint.2.h5', model)
test.main()
def test_wwinp_cylindrical():

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

60
tests/regression_tests/weightwindows_pulse_height/local/inputs_true.dat Normal file
View file

@ -0,0 +1,60 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1">
<density value="3.7" units="g/cc"/>
<nuclide name="Na23" ao="1.0"/>
<nuclide name="I127" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" name="inner sphere" material="1" region="-1" universe="1"/>
<cell id="2" name="outer sphere" material="void" region="1 -2" universe="1"/>
<surface id="1" type="sphere" coeffs="0.0 0.0 0.0 1"/>
<surface id="2" type="sphere" boundary="vacuum" coeffs="0.0 0.0 0.0 2"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>5</batches>
<source type="independent" strength="1.0" particle="photon">
<energy type="discrete">
<parameters>1000000.0 1.0</parameters>
</energy>
</source>
<photon_transport>true</photon_transport>
<weight_windows id="1">
<mesh>1</mesh>
<particle_type>photon</particle_type>
<energy_bounds>0.0 2000000.0</energy_bounds>
<lower_ww_bounds>0.01</lower_ww_bounds>
<upper_ww_bounds>0.05</upper_ww_bounds>
<survival_ratio>3.0</survival_ratio>
<max_split>10</max_split>
<weight_cutoff>1e-38</weight_cutoff>
</weight_windows>
<mesh id="1">
<dimension>1 1 1</dimension>
<lower_left>-2 -2 -2</lower_left>
<upper_right>2 2 2</upper_right>
</mesh>
<shared_secondary_bank>false</shared_secondary_bank>
<weight_window_checkpoints>
<collision>true</collision>
<surface>true</surface>
</weight_window_checkpoints>
<max_history_splits>50</max_history_splits>
</settings>
<tallies>
<filter id="1" type="cell">
<bins>1</bins>
</filter>
<filter id="2" type="energy">
<bins>0.0 10000.0 20000.0 30000.0 40000.0 50000.0 60000.0 70000.0 80000.0 90000.0 100000.0 110000.0 120000.0 130000.0 140000.0 150000.0 160000.0 170000.0 180000.0 190000.0 200000.0 210000.0 220000.0 230000.0 240000.0 250000.0 260000.0 270000.0 280000.0 290000.0 300000.0 310000.0 320000.0 330000.0 340000.0 350000.0 360000.0 370000.0 380000.0 390000.0 400000.0 410000.0 420000.0 430000.0 440000.0 450000.0 460000.0 470000.0 480000.0 490000.0 500000.0 510000.0 520000.0 530000.0 540000.0 550000.0 560000.0 570000.0 580000.0 590000.0 600000.0 610000.0 620000.0 630000.0 640000.0 650000.0 660000.0 670000.0 680000.0 690000.0 700000.0 710000.0 720000.0 730000.0 740000.0 750000.0 760000.0 770000.0 780000.0 790000.0 800000.0 810000.0 820000.0 830000.0 840000.0 850000.0 860000.0 870000.0 880000.0 890000.0 900000.0 910000.0 920000.0 930000.0 940000.0 950000.0 960000.0 970000.0 980000.0 990000.0 1000000.0</bins>
</filter>
<tally id="1" name="pht tally">
<filters>1 2</filters>
<scores>pulse-height</scores>
</tally>
</tallies>
</model>

201
tests/regression_tests/weightwindows_pulse_height/local/results_true.dat Normal file
View file

@ -0,0 +1,201 @@
tally 1:
4.140000E+00
3.443000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
3.000000E-02
5.000000E-04
2.000000E-02
4.000000E-04
2.000000E-02
4.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.000000E-02
3.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.000000E-02
6.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
2.000000E-02
4.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
3.000000E-02
5.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
2.000000E-02
2.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
2.000000E-02
2.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
3.000000E-02
5.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-01
8.600000E-03

60
tests/regression_tests/weightwindows_pulse_height/shared/inputs_true.dat Normal file
View file

@ -0,0 +1,60 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1">
<density value="3.7" units="g/cc"/>
<nuclide name="Na23" ao="1.0"/>
<nuclide name="I127" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" name="inner sphere" material="1" region="-1" universe="1"/>
<cell id="2" name="outer sphere" material="void" region="1 -2" universe="1"/>
<surface id="1" type="sphere" coeffs="0.0 0.0 0.0 1"/>
<surface id="2" type="sphere" boundary="vacuum" coeffs="0.0 0.0 0.0 2"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>5</batches>
<source type="independent" strength="1.0" particle="photon">
<energy type="discrete">
<parameters>1000000.0 1.0</parameters>
</energy>
</source>
<photon_transport>true</photon_transport>
<weight_windows id="1">
<mesh>1</mesh>
<particle_type>photon</particle_type>
<energy_bounds>0.0 2000000.0</energy_bounds>
<lower_ww_bounds>0.01</lower_ww_bounds>
<upper_ww_bounds>0.05</upper_ww_bounds>
<survival_ratio>3.0</survival_ratio>
<max_split>10</max_split>
<weight_cutoff>1e-38</weight_cutoff>
</weight_windows>
<mesh id="1">
<dimension>1 1 1</dimension>
<lower_left>-2 -2 -2</lower_left>
<upper_right>2 2 2</upper_right>
</mesh>
<shared_secondary_bank>true</shared_secondary_bank>
<weight_window_checkpoints>
<collision>true</collision>
<surface>true</surface>
</weight_window_checkpoints>
<max_history_splits>50</max_history_splits>
</settings>
<tallies>
<filter id="1" type="cell">
<bins>1</bins>
</filter>
<filter id="2" type="energy">
<bins>0.0 10000.0 20000.0 30000.0 40000.0 50000.0 60000.0 70000.0 80000.0 90000.0 100000.0 110000.0 120000.0 130000.0 140000.0 150000.0 160000.0 170000.0 180000.0 190000.0 200000.0 210000.0 220000.0 230000.0 240000.0 250000.0 260000.0 270000.0 280000.0 290000.0 300000.0 310000.0 320000.0 330000.0 340000.0 350000.0 360000.0 370000.0 380000.0 390000.0 400000.0 410000.0 420000.0 430000.0 440000.0 450000.0 460000.0 470000.0 480000.0 490000.0 500000.0 510000.0 520000.0 530000.0 540000.0 550000.0 560000.0 570000.0 580000.0 590000.0 600000.0 610000.0 620000.0 630000.0 640000.0 650000.0 660000.0 670000.0 680000.0 690000.0 700000.0 710000.0 720000.0 730000.0 740000.0 750000.0 760000.0 770000.0 780000.0 790000.0 800000.0 810000.0 820000.0 830000.0 840000.0 850000.0 860000.0 870000.0 880000.0 890000.0 900000.0 910000.0 920000.0 930000.0 940000.0 950000.0 960000.0 970000.0 980000.0 990000.0 1000000.0</bins>
</filter>
<tally id="1" name="pht tally">
<filters>1 2</filters>
<scores>pulse-height</scores>
</tally>
</tallies>
</model>

201
tests/regression_tests/weightwindows_pulse_height/shared/results_true.dat Normal file
View file

@ -0,0 +1,201 @@
tally 1:
4.140000E+00
3.443000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
3.000000E-02
5.000000E-04
2.000000E-02
4.000000E-04
2.000000E-02
4.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.000000E-02
3.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.000000E-02
6.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
2.000000E-02
4.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
3.000000E-02
5.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
2.000000E-02
2.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
2.000000E-02
2.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
3.000000E-02
5.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
1.000000E-02
1.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-02
2.000000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.000000E-01
8.600000E-03

73
tests/regression_tests/weightwindows_pulse_height/test.py Normal file
View file

@ -0,0 +1,73 @@
import numpy as np
import openmc
import pytest
from openmc.utility_funcs import change_directory
from tests.testing_harness import PyAPITestHarness
@pytest.mark.parametrize("shared_secondary,subdir", [
(False, "local"),
(True, "shared"),
])
def test_weightwindows_pulse_height(shared_secondary, subdir):
with change_directory(subdir):
openmc.reset_auto_ids()
model = openmc.Model()
# Define materials (NaI scintillator)
NaI = openmc.Material()
NaI.set_density('g/cc', 3.7)
NaI.add_element('Na', 1.0)
NaI.add_element('I', 1.0)
model.materials = openmc.Materials([NaI])
# Define geometry: NaI sphere inside vacuum sphere
s1 = openmc.Sphere(r=1)
s2 = openmc.Sphere(r=2, boundary_type='vacuum')
inner_sphere = openmc.Cell(name='inner sphere', fill=NaI, region=-s1)
outer_sphere = openmc.Cell(name='outer sphere', region=+s1 & -s2)
model.geometry = openmc.Geometry([inner_sphere, outer_sphere])
# Define settings
model.settings.run_mode = 'fixed source'
model.settings.batches = 5
model.settings.particles = 100
model.settings.photon_transport = True
model.settings.shared_secondary_bank = shared_secondary
model.settings.max_history_splits = 50
model.settings.weight_window_checkpoints = {
'surface': True,
'collision': True,
}
model.settings.source = openmc.IndependentSource(
energy=openmc.stats.delta_function(1e6),
particle='photon'
)
# Define pulse-height tally
tally = openmc.Tally(name="pht tally")
tally.scores = ['pulse-height']
cell_filter = openmc.CellFilter(inner_sphere)
energy_filter = openmc.EnergyFilter(np.linspace(0, 1_000_000, 101))
tally.filters = [cell_filter, energy_filter]
model.tallies = [tally]
# Define weight windows on a simple mesh covering the geometry
ww_mesh = openmc.RegularMesh()
ww_mesh.lower_left = (-2, -2, -2)
ww_mesh.upper_right = (2, 2, 2)
ww_mesh.dimension = (1, 1, 1)
# Single energy bin for photons
e_bnds = [0.0, 2e6]
# Uniform weight window bounds (low enough to trigger some splitting)
lower_bounds = np.array([0.01])
ww = openmc.WeightWindows(ww_mesh, lower_bounds, None, 5.0, e_bnds, 'photon')
model.settings.weight_windows = [ww]
harness = PyAPITestHarness('statepoint.5.h5', model)
harness.main()

253
tests/unit_tests/dagmc/test_model.py
View file

@ -70,28 +70,20 @@ def model(request):
openmc.reset_auto_ids()
def test_dagmc_replace_material_assignment(model):
mats = {}
mats["foo"] = openmc.Material(name="foo")
mats["foo"].add_nuclide("H1", 2.0)
mats["foo"].add_element("O", 1.0)
mats["foo"].set_density("g/cm3", 1.0)
mats["foo"].add_s_alpha_beta("c_H_in_H2O")
def test_dagmc_sync_cell_names(model):
dag_univ = None
for univ in model.geometry.get_all_universes().values():
if not isinstance(univ, openmc.DAGMCUniverse):
if isinstance(univ, openmc.DAGMCUniverse):
dag_univ = univ
break
cells_with_41 = []
for cell in univ.cells.values():
if cell.fill is None:
continue
if cell.fill.name == "41":
cells_with_41.append(cell.id)
univ.replace_material_assignment("41", mats["foo"])
for cell_id in cells_with_41:
assert univ.cells[cell_id] == mats["foo"]
assert dag_univ is not None
for cell_id, cell in dag_univ.cells.items():
assert cell.name == openmc.lib.cells[cell_id].name
assert any(cell.name == "implicit complement"
for cell in dag_univ.cells.values())
def test_dagmc_add_material_override_with_id(model):
@ -114,7 +106,7 @@ def test_dagmc_add_material_override_with_id(model):
cells_with_41.append(cell.id)
univ.add_material_override(cell.id, mats["foo"])
for cell_id in cells_with_41:
assert univ.cells[cell_id] == mats["foo"]
assert univ.cells[cell_id].fill == mats["foo"]
def test_dagmc_add_material_override_with_cell(model):
@ -137,7 +129,7 @@ def test_dagmc_add_material_override_with_cell(model):
cells_with_41.append(cell.id)
univ.add_material_override(cell, mats["foo"])
for cell_id in cells_with_41:
assert univ.cells[cell_id] == mats["foo"]
assert univ.cells[cell_id].fill == mats["foo"]
def test_model_differentiate_depletable_with_dagmc(model, run_in_tmpdir):
@ -174,57 +166,20 @@ def test_model_differentiate_with_dagmc(model):
assert len(model.materials) == 4*2 + 4
def test_bad_override_cell_id(model):
for univ in model.geometry.get_all_universes().values():
if isinstance(univ, openmc.DAGMCUniverse):
break
with pytest.raises(ValueError, match="Cell ID '1' not found in DAGMC universe"):
univ.material_overrides = {1: model.materials[0]}
def test_bad_override_type(model):
not_a_dag_cell = openmc.Cell()
for univ in model.geometry.get_all_universes().values():
if isinstance(univ, openmc.DAGMCUniverse):
break
with pytest.raises(ValueError, match="Unrecognized key type. Must be an integer or openmc.DAGMCCell object"):
univ.material_overrides = {not_a_dag_cell: model.materials[0]}
def test_bad_replacement_mat_name(model):
for univ in model.geometry.get_all_universes().values():
if isinstance(univ, openmc.DAGMCUniverse):
break
with pytest.raises(ValueError, match="No material with name 'not_a_mat' found in the DAGMC universe"):
univ.replace_material_assignment("not_a_mat", model.materials[0])
def test_dagmc_xml(model):
# Set the environment
mats = {}
mats["no-void fuel"] = openmc.Material(1, name="no-void fuel")
mats["no-void fuel"].add_nuclide("U235", 0.03)
mats["no-void fuel"].add_nuclide("U238", 0.97)
mats["no-void fuel"].add_nuclide("O16", 2.0)
mats["no-void fuel"].set_density("g/cm3", 10.0)
mats[5] = openmc.Material(name="41")
mats[5].add_nuclide("H1", 2.0)
mats[5].add_element("O", 1.0)
mats[5].set_density("g/cm3", 1.0)
mats[5].add_s_alpha_beta("c_H_in_H2O")
override_mat = openmc.Material(name="41")
override_mat.add_nuclide("H1", 2.0)
override_mat.add_element("O", 1.0)
override_mat.set_density("g/cm3", 1.0)
override_mat.add_s_alpha_beta("c_H_in_H2O")
model.materials.append(override_mat)
for univ in model.geometry.get_all_universes().values():
if isinstance(univ, openmc.DAGMCUniverse):
dag_univ = univ
break
for k, v in mats.items():
if isinstance(k, int):
dag_univ.add_material_override(k, v)
model.materials.append(v)
elif isinstance(k, str):
dag_univ.replace_material_assignment(k, v)
dag_univ.add_material_override(5, override_mat)
# Tesing the XML subelement generation
root = ET.Element('dagmc_universe')
@ -236,12 +191,24 @@ def test_dagmc_xml(model):
assert dagmc_ele.get('filename') == str(dag_univ.filename)
assert dagmc_ele.get('auto_geom_ids') == str(dag_univ.auto_geom_ids).lower()
override_eles = dagmc_ele.find('material_overrides').findall('cell_override')
assert len(override_eles) == 4
assert dagmc_ele.find('material_overrides') is None
for i, override_ele in enumerate(override_eles):
cell_id = override_ele.get('id')
assert dag_univ.material_overrides[int(cell_id)][0].id == int(override_ele.find('material_ids').text)
override_elements = dagmc_ele.findall('cell')
assert len(override_elements) == len(dag_univ.cells)
xml_cells = {int(elem.get('id')): elem for elem in override_elements}
for cell_id, cell in dag_univ.cells.items():
assert cell_id in xml_cells
xml_cell = xml_cells[cell_id]
if cell.fill_type == 'void':
assert xml_cell.get('material') == 'void'
elif cell.fill_type == 'material':
assert xml_cell.get('material') == str(cell.fill.id)
elif cell.fill_type == 'distribmat':
mat_list = xml_cell.find('material').text.split()
expected = ["void" if m is None else str(m.id) for m in cell.fill]
assert mat_list == expected
else:
pytest.fail(f"Unexpected DAGMC cell fill type: {cell.fill_type}")
model.export_to_model_xml()
@ -252,7 +219,147 @@ def test_dagmc_xml(model):
xml_dagmc_univ = univ
break
assert xml_dagmc_univ._material_overrides.keys() == dag_univ._material_overrides.keys()
assert xml_dagmc_univ.cells.keys() == dag_univ.cells.keys()
for xml_mats, model_mats in zip(xml_dagmc_univ._material_overrides.values(), dag_univ._material_overrides.values()):
assert all([xml_mat.id == orig_mat.id for xml_mat, orig_mat in zip(xml_mats, model_mats)])
for cell_id, cell in dag_univ.cells.items():
xml_cell = xml_dagmc_univ.cells[cell_id]
assert xml_cell.fill_type == cell.fill_type
if cell.fill_type == 'void':
assert xml_cell.fill is None
elif cell.fill_type == 'material':
assert xml_cell.fill.id == cell.fill.id
elif cell.fill_type == 'distribmat':
xml_ids = [m.id if m is not None else None for m in xml_cell.fill]
model_ids = [m.id if m is not None else None for m in cell.fill]
assert xml_ids == model_ids
else:
pytest.fail(f"Unexpected DAGMC cell fill type: {cell.fill_type}")
def test_dagmc_xml_reject_fill_override():
mats = {'1': openmc.Material(1), 'void': None}
elem = ET.fromstring(
'<dagmc_universe id="1" filename="dagmc.h5m">'
'<cell id="1" fill="2"/>'
'</dagmc_universe>'
)
with pytest.raises(ValueError, match="cannot specify 'fill'"):
openmc.DAGMCUniverse.from_xml_element(elem, mats)
def test_dagmc_xml_reject_region_override():
mats = {'1': openmc.Material(1), 'void': None}
elem = ET.fromstring(
'<dagmc_universe id="1" filename="dagmc.h5m">'
'<cell id="1" material="1" region="-1"/>'
'</dagmc_universe>'
)
with pytest.raises(ValueError, match="cannot specify 'region'"):
openmc.DAGMCUniverse.from_xml_element(elem, mats)
def _legacy_xml(cell_overrides):
"""Helper to build a <dagmc_universe> with old-format <material_overrides>."""
inner = ''.join(
f'<cell_override id="{cid}"><material_ids>{mids}</material_ids></cell_override>'
for cid, mids in cell_overrides.items()
)
return ET.fromstring(
f'<dagmc_universe id="1" filename="dagmc.h5m">'
f'<material_overrides>{inner}</material_overrides>'
f'</dagmc_universe>'
)
def test_dagmc_xml_legacy_single_material_compat():
mat = openmc.Material(1)
mats = {'1': mat, 'void': None}
elem = _legacy_xml({3: '1'})
with pytest.warns(DeprecationWarning, match="deprecated"):
univ = openmc.DAGMCUniverse.from_xml_element(elem, mats)
assert 3 in univ.cells
assert univ.cells[3].fill is mat
def test_dagmc_xml_legacy_distribmat_compat():
mat1, mat2 = openmc.Material(2), openmc.Material(3)
mats = {'2': mat1, '3': mat2, 'void': None}
elem = _legacy_xml({5: '2 3'})
with pytest.warns(DeprecationWarning):
univ = openmc.DAGMCUniverse.from_xml_element(elem, mats)
assert univ.cells[5].fill_type == 'distribmat'
assert list(univ.cells[5].fill) == [mat1, mat2]
def test_dagmc_xml_legacy_void_compat():
mats = {'void': None}
elem = _legacy_xml({7: 'void'})
with pytest.warns(DeprecationWarning):
univ = openmc.DAGMCUniverse.from_xml_element(elem, mats)
assert univ.cells[7].fill_type == 'void'
def test_dagmc_xml_legacy_both_raises():
mat = openmc.Material(1)
mats = {'1': mat, 'void': None}
elem = ET.fromstring(
'<dagmc_universe id="1" filename="dagmc.h5m">'
'<material_overrides>'
'<cell_override id="3"><material_ids>1</material_ids></cell_override>'
'</material_overrides>'
'<cell id="5" material="1"/>'
'</dagmc_universe>'
)
with pytest.raises(ValueError, match="both"):
openmc.DAGMCUniverse.from_xml_element(elem, mats)
def test_dagmc_xml_legacy_deprecation_warning():
mats = {'1': openmc.Material(1), 'void': None}
elem = _legacy_xml({3: '1'})
with pytest.warns(DeprecationWarning):
openmc.DAGMCUniverse.from_xml_element(elem, mats)
def test_dagmc_xml_legacy_roundtrip():
"""Old-format XML loads correctly and re-exports using the new <cell> format."""
mat = openmc.Material(1)
mats = {'1': mat, 'void': None}
elem = _legacy_xml({3: '1'})
with pytest.warns(DeprecationWarning):
univ = openmc.DAGMCUniverse.from_xml_element(elem, mats)
root = ET.Element('geometry')
univ.create_xml_subelement(root)
dagmc_elem = root.find('dagmc_universe')
assert dagmc_elem.find('material_overrides') is None
cell_elems = dagmc_elem.findall('cell')
assert len(cell_elems) == 1
assert int(cell_elems[0].get('id')) == 3
assert cell_elems[0].get('material') == '1'
def test_dagmc_xml_temperature_roundtrip():
mat = openmc.Material(1)
mats = {'1': mat, 'void': None}
elem = ET.fromstring(
'<dagmc_universe id="10" filename="dagmc.h5m">'
'<cell id="7" material="1" temperature="825.0"/>'
'</dagmc_universe>'
)
dag_univ = openmc.DAGMCUniverse.from_xml_element(elem, mats)
assert dag_univ.cells[7].fill.id == 1
assert dag_univ.cells[7].temperature == pytest.approx(825.0)
root = ET.Element('geometry')
dag_univ.create_xml_subelement(root)
dagmc_elem = root.find('dagmc_universe')
xml_cell = dagmc_elem.find('cell')
assert xml_cell.get('temperature') == '825.0'
dag_univ_roundtrip = openmc.DAGMCUniverse.from_xml_element(dagmc_elem, mats)
assert dag_univ_roundtrip.cells[7].fill.id == 1
assert dag_univ_roundtrip.cells[7].temperature == pytest.approx(825.0)

43
tests/unit_tests/test_cell.py
View file

@ -374,6 +374,49 @@ def test_rotation_from_xml(rotation):
np.testing.assert_allclose(new_cell.rotation, cell.rotation)
def test_dagmccell_from_xml_element():
"""DAGMCCell.from_xml_element parses material, temperature, density,
and volume; rejects unsupported attributes."""
mat = openmc.Material(1)
mat.add_nuclide('U235', 1.0)
mats = {'1': mat}
# In practice, from_xml_element is always called from DAGMCUniverse
# during XML parsing. A placeholder universe is used here so the test
# can exercise the method directly without a real DAGMC model file.
placeholder_univ = openmc.DAGMCUniverse('model.h5m')
# material + temperature + density round-trip
xml = '<cell id="5" name="fuel" material="1" temperature="900.0" density="10.5"/>'
cell = openmc.DAGMCCell.from_xml_element(ET.fromstring(xml), mats, placeholder_univ)
assert cell.id == 5
assert cell.name == 'fuel'
assert cell.fill is mat
assert cell.density == 10.5
assert cell.temperature == 900.0
# volume round-trip
xml = '<cell id="6" material="1" volume="42.0"/>'
cell = openmc.DAGMCCell.from_xml_element(ET.fromstring(xml), mats, placeholder_univ)
assert cell.volume == 42.0
# forbidden: region, fill, universe
for tag, val in [('region', '-1'), ('fill', '2'), ('universe', '0')]:
xml = f'<cell id="7" material="1" {tag}="{val}"/>'
with pytest.raises(ValueError, match=tag):
openmc.DAGMCCell.from_xml_element(ET.fromstring(xml), mats, placeholder_univ)
# forbidden: translation, rotation
for tag, val in [('translation', '1 0 0'), ('rotation', '0 0 90')]:
xml = f'<cell id="7" material="1" {tag}="{val}"/>'
with pytest.raises(ValueError, match=tag):
openmc.DAGMCCell.from_xml_element(ET.fromstring(xml), mats, placeholder_univ)
# missing material raises
xml = '<cell id="8"/>'
with pytest.raises(ValueError, match='material'):
openmc.DAGMCCell.from_xml_element(ET.fromstring(xml), mats, placeholder_univ)
def test_plot(run_in_tmpdir):
zcyl = openmc.ZCylinder()
c = openmc.Cell(region=-zcyl)

6
tests/unit_tests/test_material.py
View file

@ -561,6 +561,10 @@ def test_get_activity():
m1.add_element("Fe", 0.7)
m1.add_element("Li", 0.3)
m1.set_density('g/cm3', 1.5)
with pytest.raises(ValueError, match="Volume must be set"):
m1.get_activity(units='Bq')
with pytest.raises(ValueError, match="Volume must be set"):
m1.get_activity(units='Ci')
# activity in Bq/cc and Bq/g should not require volume setting
assert m1.get_activity(units='Bq/cm3') == 0
assert m1.get_activity(units='Bq/g') == 0
@ -619,6 +623,8 @@ def test_get_decay_heat():
m1.add_nuclide("U235", 0.2)
m1.add_nuclide("U238", 0.8)
m1.set_density('g/cm3', 10.5)
with pytest.raises(ValueError, match="Volume must be set"):
m1.get_decay_heat(units='W')
# decay heat in W/cc and W/g should not require volume setting
assert m1.get_decay_heat(units='W/cm3') == 0
assert m1.get_decay_heat(units='W/g') == 0

49
tests/unit_tests/test_materials.py
View file

@ -1,5 +1,7 @@
from pathlib import Path
import pytest
import openmc
from openmc.deplete import Chain
@ -78,3 +80,50 @@ def test_export_duplicate_materials_to_xml(run_in_tmpdir):
materials_in = openmc.Materials.from_xml("materials.xml")
assert len(materials_in) == 2
def test_materials_deplete_length_mismatch():
mats = openmc.Materials([openmc.Material()])
with pytest.raises(ValueError, match="multigroup_fluxes length"):
mats.deplete(
multigroup_fluxes=[],
energy_group_structures=["VITAMIN-J-42"],
timesteps=[1.0],
source_rates=1.0,
)
with pytest.raises(ValueError, match="energy_group_structures length"):
mats.deplete(
multigroup_fluxes=[[1.0]],
energy_group_structures=[],
timesteps=[1.0],
source_rates=1.0,
)
def test_materials_deplete_missing_volume(monkeypatch):
mat = openmc.Material()
mat.add_nuclide("Ni58", 1.0)
mat.set_density("g/cm3", 7.87)
mats = openmc.Materials([mat])
class DummySession:
def __enter__(self):
return self
def __exit__(self, exc_type, exc, tb):
return False
monkeypatch.setattr(openmc.lib, "TemporarySession", DummySession)
chain = Path(__file__).parents[1] / "chain_ni.xml"
with pytest.raises(ValueError, match="has no volume"):
mats.deplete(
multigroup_fluxes=[[1.0]],
energy_group_structures=["VITAMIN-J-42"],
timesteps=[1.0],
source_rates=1.0,
chain_file=chain,
)

45
tests/unit_tests/test_mesh.py
View file

@ -490,13 +490,28 @@ def test_umesh(run_in_tmpdir, simple_umesh, export_type):
simple_umesh.write_data_to_vtk(datasets={'mean': ref_data[:-2]}, filename=filename)
@pytest.mark.skipif(not openmc.lib._dagmc_enabled(), reason="DAGMC not enabled.")
def test_write_vtkhdf(request, run_in_tmpdir):
vtkhdf_tests = [
(
Path("test_mesh_dagmc_tets.vtk"),
"moab"
),
(
Path("test_mesh_hexes.exo"),
"libmesh"
)
]
@pytest.mark.parametrize('mesh_file, mesh_library', vtkhdf_tests)
def test_write_vtkhdf(mesh_file, mesh_library, request, run_in_tmpdir):
"""Performs a minimal UnstructuredMesh simulation, reads in the resulting
statepoint file and writes the mesh data to vtk and vtkhdf files. It is
necessary to read in the unstructured mesh from a statepoint file to ensure
it has all the required attributes
"""
if mesh_library == 'moab' and not openmc.lib._dagmc_enabled():
pytest.skip("DAGMC not enabled.")
if mesh_library == 'libmesh' and not openmc.lib._libmesh_enabled():
pytest.skip("LibMesh not enabled.")
model = openmc.Model()
surf1 = openmc.Sphere(r=1000.0, boundary_type="vacuum")
@ -504,8 +519,8 @@ def test_write_vtkhdf(request, run_in_tmpdir):
model.geometry = openmc.Geometry([cell1])
umesh = openmc.UnstructuredMesh(
request.path.parent / "test_mesh_dagmc_tets.vtk",
"moab",
request.path.parent / mesh_file,
mesh_library,
mesh_id = 1
)
mesh_filter = openmc.MeshFilter(umesh)
@ -514,6 +529,8 @@ def test_write_vtkhdf(request, run_in_tmpdir):
mesh_tally = openmc.Tally(name="test_tally")
mesh_tally.filters = [mesh_filter]
mesh_tally.scores = ["flux"]
if mesh_library == "libmesh":
mesh_tally.estimator = "collision"
model.tallies = [mesh_tally]
@ -556,6 +573,26 @@ def test_write_vtkhdf(request, run_in_tmpdir):
with h5py.File("test_mesh.vtkhdf", "r"):
...
import vtk
reader = vtk.vtkHDFReader()
reader.SetFileName("test_mesh.vtkhdf")
reader.Update()
# Get mean from file and make sure it matches original data
num_elements = reader.GetOutput().GetNumberOfCells()
assert num_elements == umesh_from_sp.n_elements
num_vertices = reader.GetOutput().GetNumberOfPoints()
assert num_vertices == umesh_from_sp.vertices.shape[0]
arr = reader.GetOutput().GetCellData().GetArray("mean")
mean = np.array([arr.GetTuple1(i) for i in range(my_tally.mean.size)])
np.testing.assert_almost_equal(mean, my_tally.mean.flatten()/umesh_from_sp.volumes)
arr = reader.GetOutput().GetCellData().GetArray("std_dev")
std_dev = np.array([arr.GetTuple1(i) for i in range(my_tally.std_dev.size)])
np.testing.assert_almost_equal(std_dev, my_tally.std_dev.flatten()/umesh_from_sp.volumes)
def test_mesh_get_homogenized_materials():
"""Test the get_homogenized_materials method"""

1
tests/unit_tests/test_mesh_hexes.exo Symbolic link
View file

@ -0,0 +1 @@
../regression_tests/unstructured_mesh/test_mesh_hexes.exo

69
tests/unit_tests/test_source.py
View file

@ -1,5 +1,6 @@
from collections import Counter
from math import pi
from pathlib import Path
import openmc
import openmc.lib
@ -91,6 +92,74 @@ def test_point_cloud_strengths(run_in_tmpdir, sphere_box_model):
assert sampled_strength == expected_strength, f'Strength incorrect for {positions[i]}'
def test_decay_spectrum_parent_nuclide(run_in_tmpdir):
chain_file = Path('chain_decay_spectrum_parent.xml')
chain_file.write_text("""<?xml version="1.0"?>
<depletion_chain>
<nuclide name="ParentA" decay_modes="0" reactions="0" half_life="1.0">
<source type="discrete" particle="photon">
<parameters>1000000.0 1.0</parameters>
</source>
</nuclide>
<nuclide name="ParentB" decay_modes="0" reactions="0" half_life="1.0">
<source type="discrete" particle="photon">
<parameters>2000000.0 1.0</parameters>
</source>
</nuclide>
</depletion_chain>
""")
inner_sphere = openmc.Sphere(r=10.0)
outer_sphere = openmc.Sphere(r=20.0, boundary_type='vacuum')
shell_mat = openmc.Material()
shell_mat.add_nuclide('H1', 1.0)
shell_mat.set_density('atom/b-cm', 1.0e-12)
void_cell = openmc.Cell(region=-inner_sphere)
shell_cell = openmc.Cell(fill=shell_mat, region=+inner_sphere & -outer_sphere)
model = openmc.Model()
model.geometry = openmc.Geometry([void_cell, shell_cell])
model.materials = [shell_mat]
model.settings.run_mode = 'fixed source'
model.settings.photon_transport = True
model.settings.particles = 1000
model.settings.batches = 5
model.settings.source = openmc.IndependentSource(
particle='photon',
space=openmc.stats.Point((0.0, 0.0, 0.0)),
energy=openmc.stats.DecaySpectrum(
{'ParentA': 1.0, 'ParentB': 1.0},
volume=1.0
)
)
tally = openmc.Tally()
tally.filters = [
openmc.CellFilter([void_cell]),
openmc.ParticleFilter(['photon']),
openmc.EnergyFilter([0.0, 1.5e6, 2.5e6]),
openmc.ParentNuclideFilter(['ParentA', 'ParentB'])
]
tally.scores = ['flux']
model.tallies = [tally]
with openmc.config.patch('chain_file', chain_file):
sp_filename = model.run()
with openmc.StatePoint(sp_filename) as sp:
tally_out = sp.tallies[tally.id]
mean = tally_out.get_reshaped_data('mean').squeeze()
assert mean.shape == (2, 2)
assert mean[0, 0] > 0.0
assert mean[1, 1] > 0.0
assert mean[0, 1] == 0.0
assert mean[1, 0] == 0.0
assert np.count_nonzero(mean) == 2
def test_source_file():
filename = 'source.h5'
src = openmc.FileSource(path=filename)

135
tests/unit_tests/test_stats.py
View file

@ -1,10 +1,11 @@
from math import pi
from pathlib import Path
import numpy as np
import pytest
import openmc
import openmc.stats
from openmc.stats.univariate import _INTERPOLATION_SCHEMES
from openmc.stats.univariate import _INTERPOLATION_SCHEMES, DecaySpectrum
from scipy.integrate import trapezoid
from tests.unit_tests import assert_sample_mean
@ -1013,3 +1014,135 @@ def test_fusion_spectrum_invalid():
# Temperature above 100 keV should raise an error
with pytest.raises(ValueError):
openmc.stats.fusion_neutron_spectrum(101e3, 'DT')
@pytest.fixture(autouse=False)
def decay_spectrum_chain():
"""Set chain_file for the duration of a test and clear the _photon_integral
cache so results from a different chain don't bleed across tests."""
CHAIN_FILE = (Path(__file__).parents[1] / 'chain_simple.xml').resolve()
DecaySpectrum._photon_integral.cache_clear()
with openmc.config.patch('chain_file', CHAIN_FILE):
yield
DecaySpectrum._photon_integral.cache_clear()
def test_decay_spectrum_construction():
nuclides = {'I135': 1.5e-3, 'Xe135': 8.2e-4}
d = openmc.stats.DecaySpectrum(nuclides, volume=100.0)
assert d.nuclides == nuclides
assert d.volume == pytest.approx(100.0)
assert len(d) == 2
def test_decay_spectrum_validation():
# nuclides must be a dict
with pytest.raises(TypeError):
openmc.stats.DecaySpectrum(['I135'], volume=1.0)
# densities must be > 0
with pytest.raises(ValueError):
openmc.stats.DecaySpectrum({'I135': -1.0}, volume=1.0)
# volume must be > 0
with pytest.raises(ValueError):
openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=-1.0)
with pytest.raises(ValueError):
openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=0.0)
def test_decay_spectrum_xml_roundtrip():
nuclides = {'I135': 1.5e-3, 'Xe135': 8.2e-4}
d = openmc.stats.DecaySpectrum(nuclides, volume=100.0)
elem = d.to_xml_element('energy')
assert elem.get('type') == 'decay_spectrum'
assert float(elem.get('volume')) == pytest.approx(100.0)
assert elem.findtext('nuclides').split() == list(nuclides)
assert [float(x) for x in elem.findtext('parameters').split()] == pytest.approx(
list(nuclides.values()))
# Round-trip via DecaySpectrum.from_xml_element
d2 = openmc.stats.DecaySpectrum.from_xml_element(elem)
assert d2.nuclides == nuclides
assert d2.volume == pytest.approx(100.0)
# Round-trip via the Univariate dispatcher
d3 = openmc.stats.Univariate.from_xml_element(elem)
assert isinstance(d3, openmc.stats.DecaySpectrum)
assert d3 == d
def test_decay_spectrum_to_distribution(decay_spectrum_chain):
# Single emitting nuclide -> concrete distribution, not None
d = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=10.0)
dist = d.to_distribution()
assert dist is not None
# Result is cached on second call
dist2 = d.to_distribution()
assert dist2 is dist
# Nuclide with no photon source -> None
d_stable = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0)
assert d_stable.to_distribution() is None
# Mixture of emitters -> non-None combined distribution
d_mix = openmc.stats.DecaySpectrum(
{'I135': 1e-3, 'Xe135': 5e-4}, volume=10.0
)
dist_mix = d_mix.to_distribution()
assert dist_mix is not None
def test_decay_spectrum_integral(decay_spectrum_chain):
# For an emitting nuclide, integral should be > 0
d = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=10.0)
assert d.integral() > 0.0
# Proportional to density: doubling density doubles integral
d2 = openmc.stats.DecaySpectrum({'I135': 2e-3}, volume=10.0)
assert d2.integral() == pytest.approx(2.0 * d.integral())
# Proportional to volume
d3 = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=20.0)
assert d3.integral() == pytest.approx(2.0 * d.integral())
# Pure non-emitter -> 0.0
d_stable = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0)
assert d_stable.integral() == pytest.approx(0.0)
def test_decay_spectrum_clip(decay_spectrum_chain):
# Stable / non-emitting nuclides are removed unconditionally
d = openmc.stats.DecaySpectrum(
{'I135': 1e-3, 'Xe135': 5e-4, 'Xe136': 1.0, 'Cs135': 1.0},
volume=10.0,
)
d_clip = d.clip()
assert 'Xe136' not in d_clip.nuclides
assert 'Cs135' not in d_clip.nuclides
assert 'I135' in d_clip.nuclides
assert 'Xe135' in d_clip.nuclides
# Original is unchanged
assert 'Xe136' in d.nuclides
# inplace=True modifies and returns the same object
d_same = d.clip(inplace=True)
assert d_same is d
assert 'Xe136' not in d.nuclides
# A nuclide with negligible emission rate is removed by tolerance clipping.
# U235 has a very small integral (~4e-17 Bq/atom) compared with I135 (~4e-5)
d_tight = openmc.stats.DecaySpectrum(
{'I135': 1e-3, 'U235': 1e-3}, volume=10.0
)
d_tight_clip = d_tight.clip(tolerance=1e-9)
assert 'U235' not in d_tight_clip.nuclides
assert 'I135' in d_tight_clip.nuclides
# All non-emitters -> empty nuclides dict
d_empty = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0)
d_empty.clip(inplace=True)
assert d_empty.nuclides == {}

16
tests/unit_tests/weightwindows/test.py
View file

@ -122,7 +122,8 @@ def model():
return model
def test_weightwindows(model, wws):
@pytest.mark.parametrize("shared_secondary", [False, True])
def test_weightwindows(model, wws, shared_secondary):
ww_files = ('ww_n.txt', 'ww_p.txt')
cwd = Path(__file__).parent.absolute()
@ -131,6 +132,7 @@ def test_weightwindows(model, wws):
with cdtemp(filepaths):
# run once with variance reduction off
model.settings.weight_windows_on = False
model.settings.shared_secondary_bank = shared_secondary
analog_sp = model.run()
os.rename(analog_sp, 'statepoint.analog.h5')
@ -223,7 +225,8 @@ def test_lower_ww_bounds_shape():
assert ww.lower_ww_bounds.shape == (2, 3, 4, 1)
def test_photon_heating(run_in_tmpdir):
@pytest.mark.parametrize("shared_secondary", [False, True])
def test_photon_heating(run_in_tmpdir, shared_secondary):
water = openmc.Material()
water.add_nuclide('H1', 1.0)
water.add_nuclide('O16', 2.0)
@ -246,7 +249,8 @@ def test_photon_heating(run_in_tmpdir):
model.settings.run_mode = 'fixed source'
model.settings.batches = 5
model.settings.particles = 100
model.settings.particles = 101
model.settings.shared_secondary_bank = shared_secondary
tally = openmc.Tally()
tally.scores = ['heating']
@ -260,7 +264,11 @@ def test_photon_heating(run_in_tmpdir):
with openmc.StatePoint(sp_file) as sp:
tally_mean = sp.tallies[tally.id].mean
# these values should be nearly identical
# Note: Our current physics model actually does allow this tally to
# occasionally go slightly negative. However, larger bugs can
# make this more common. We have selected a particle count for
# this test that happens to produce no negative tallies for both
# the shared and non-shared secondary PRNG streams.
assert np.all(tally_mean >= 0)