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

Merge branch 'openmc-dev:develop' into mesh_tally_amalgamation

Browse source
This commit is contained in:
Ebny Walid Ahammed 2026-07-07 14:07:21 -05:00 • committed by GitHub
commit 1f75b45f46
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
203 changed files with 6652 additions and 1540 deletions
Download patch file Download diff file Expand all files Collapse all files

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

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

24
CMakeLists.txt
View file

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

34
cmake/OpenMCConfig.cmake.in
View file

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

11
cmake/OpenMCConfigVersion.cmake.in
View file

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

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

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

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

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

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

@ -98,6 +98,11 @@ sub-elements:
A list of strings representing the nuclide, to define specific
define specific target nuclide collisions to be banked.
.. note::
Electron and positron collision-track events are not associated with
a specific nuclide. If a ``nuclides`` entry is specified, these events
are omitted.
*Default*: None
:reactions:
@ -606,30 +611,30 @@ found in the :ref:`random ray user guide <random_ray>`.
*Default*: None
:adjoint_source:
Specifies an adjoint fixed source for adjoint transport simulations, and
follows the format for :ref:`source_element`. The distributions which make
up the adjoint source are subject to the same restrictions as forward
Specifies an adjoint fixed source for adjoint transport simulations, and
follows the format for :ref:`source_element`. The distributions which make
up the adjoint source are subject to the same restrictions as forward
fixed sources in Random Ray mode.
*Default*: None
:adjoint:
Specifies whether to perform adjoint transport. The default is 'False',
Specifies whether to perform adjoint transport. The default is 'False',
corresponding to forward transport.
*Default*: None
:volume_estimator:
Specifies choice of volume estimator for the random ray solver. Options
Specifies choice of volume estimator for the random ray solver. Options
are 'naive', 'simulation_averaged', or 'hybrid'. The default is 'hybrid'.
*Default*: None
:volume_normalized_flux_tallies:
Specifies whether to normalize flux tallies by volume (bool). The
default is 'False'. When enabled, flux tallies will be reported in units
of cm/cm^3. When disabled, flux tallies will be reported in units of cm
(i.e., total distance traveled by neutrons in the spatial tally
Specifies whether to normalize flux tallies by volume (bool). The
default is 'False'. When enabled, flux tallies will be reported in units
of cm/cm^3. When disabled, flux tallies will be reported in units of cm
(i.e., total distance traveled by neutrons in the spatial tally
region).
*Default*: None
@ -741,14 +746,16 @@ pseudo-random number generator.
*Default*: 1
--------------------
``<stride>`` Element
--------------------
-----------------------------------
``<shared_secondary_bank>`` Element
-----------------------------------
The ``stride`` element is used to specify how many random numbers are allocated
for each source particle history.
*Default*: 152,917
The ``shared_secondary_bank`` element indicates whether to use a shared
secondary particle bank. When enabled, secondary particles are collected into
a global bank, sorted for reproducibility, and load-balanced across MPI ranks
between generations. If not specified, the shared secondary bank is enabled
automatically for fixed-source simulations with weight windows active, and
disabled otherwise.
.. _source_element:
@ -1156,23 +1163,6 @@ based on constraints.
*Default*: 0.05
-------------------------
``<state_point>`` Element
-------------------------
The ``<state_point>`` element indicates at what batches a state point file
should be written. A state point file can be used to restart a run or to get
tally results at any batch. The default behavior when using this tag is to
write out the source bank in the state_point file. This behavior can be
customized by using the ``<source_point>`` element. This element has the
following attributes/sub-elements:
:batches:
A list of integers separated by spaces indicating at what batches a state
point file should be written.
*Default*: Last batch only
--------------------------
``<source_point>`` Element
--------------------------
@ -1222,6 +1212,32 @@ attributes/sub-elements:
*Default*: false
-------------------------
``<state_point>`` Element
-------------------------
The ``<state_point>`` element indicates at what batches a state point file
should be written. A state point file can be used to restart a run or to get
tally results at any batch. The default behavior when using this tag is to
write out the source bank in the state_point file. This behavior can be
customized by using the ``<source_point>`` element. This element has the
following attributes/sub-elements:
:batches:
A list of integers separated by spaces indicating at what batches a state
point file should be written.
*Default*: Last batch only
--------------------
``<stride>`` Element
--------------------
The ``stride`` element is used to specify how many random numbers are allocated
for each source particle history.
*Default*: 152,917
------------------------------
``<surf_source_read>`` Element
------------------------------
@ -1746,11 +1762,11 @@ mesh-based weight windows.
The ratio of the lower to upper weight window bounds.
*Default*: 5.0
For FW-CADIS:
:targets:
A sequence of IDs corresponding to the tallies which cover phase
A sequence of IDs corresponding to the tallies which cover phase
space regions of interest for local variance reduction.
*Default*: None

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

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

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

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

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

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

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

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

8
include/openmc/bank.h
View file

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

5
include/openmc/capi.h
View file

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

24
include/openmc/cell.h
View file

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

3
include/openmc/constants.h
View file

@ -35,7 +35,7 @@ constexpr array<int, 2> VERSION_VOXEL {2, 0};
constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
constexpr array<int, 2> VERSION_PROPERTIES {1, 1};
constexpr array<int, 2> VERSION_WEIGHT_WINDOWS {1, 0};
constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 1};
constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 2};
// ============================================================================
// ADJUSTABLE PARAMETERS
@ -368,6 +368,7 @@ enum class SolverType { MONTE_CARLO, RANDOM_RAY };
enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID };
enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY };
enum class RandomRaySampleMethod { PRNG, HALTON, S2 };
enum class RandomRaySolve { FORWARD, FORWARD_FOR_ADJOINT, ADJOINT };
//==============================================================================
// Geometry Constants

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

4
include/openmc/error.h
View file

@ -64,14 +64,14 @@ void write_message(
int level, const std::string& message, const Params&... fmt_args)
{
if (settings::verbosity >= level) {
write_message(fmt::format(message, fmt_args...));
write_message(fmt::format(fmt::runtime(message), fmt_args...));
}
}
template<typename... Params>
void write_message(const std::string& message, const Params&... fmt_args)
{
write_message(fmt::format(message, fmt_args...));
write_message(fmt::format(fmt::runtime(message), fmt_args...));
}
#ifdef OPENMC_MPI

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

38
include/openmc/geometry.h
View file

@ -3,9 +3,12 @@
#include <cmath>
#include <cstdint>
#include <unordered_map>
#include <vector>
#include "openmc/array.h"
#include "openmc/constants.h"
#include "openmc/random_ray/source_region.h" // For hash_combine
#include "openmc/vector.h"
namespace openmc {
@ -13,6 +16,34 @@ namespace openmc {
class BoundaryInfo;
class GeometryState;
//==============================================================================
//! OverlapKey to store cell and universe data of a single overlap, along with
//! a functor for hashing an OverlapKey into an unordered_map.
//==============================================================================
struct OverlapKey {
int universe_id;
int cell1_id;
int cell2_id;
bool operator==(const OverlapKey& other) const
{
return universe_id == other.universe_id && cell1_id == other.cell1_id &&
cell2_id == other.cell2_id;
}
};
struct OverlapKeyHash {
std::size_t operator()(const OverlapKey& k) const
{
size_t seed = 0;
hash_combine(seed, k.universe_id);
hash_combine(seed, k.cell1_id);
hash_combine(seed, k.cell2_id);
return seed;
}
};
//==============================================================================
// Global variables
//==============================================================================
@ -24,6 +55,10 @@ extern "C" int n_coord_levels; //!< Number of CSG coordinate levels
extern vector<int64_t> overlap_check_count;
// Overlap data structures get cleared every slice_data run
extern vector<OverlapKey> overlap_keys;
extern std::unordered_map<OverlapKey, int, OverlapKeyHash> overlap_key_index;
} // namespace model
//==============================================================================
@ -38,8 +73,7 @@ inline bool coincident(double d1, double d2)
//==============================================================================
//! Check for overlapping cells at a particle's position.
//==============================================================================
bool check_cell_overlap(GeometryState& p, bool error = true);
int check_cell_overlap(GeometryState& p, bool error = true);
//==============================================================================
//! Get the cell instance for a particle at the specified universe level

12
include/openmc/math_functions.h
View file

@ -233,5 +233,17 @@ void get_energy_index(
double standard_normal_cdf(double z);
//==============================================================================
//! Return true if two floating-point values are approximately equal within a
//! combined relative and absolute tolerance.
//!
//! \param a first floating point value
//! \param b second floating point value
//! \param rel_tol relative tolerance
//! \param abs_tol absolute tolerance
//! \return true if a and b are approximately equal, false otherwise
//==============================================================================
bool isclose(double a, double b, double rel_tol, double abs_tol);
} // namespace openmc
#endif // OPENMC_MATH_FUNCTIONS_H

3
include/openmc/nuclide.h
View file

@ -84,6 +84,9 @@ public:
double collapse_rate(int MT, double temperature, span<const double> energy,
span<const double> flux) const;
//! Return a ParticleType object representing this nuclide
ParticleType particle_type() const { return {Z_, A_, metastable_}; }
//============================================================================
// Data members
std::string name_; //!< Name of nuclide, e.g. "U235"

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

124
include/openmc/plot.h
View file

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

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

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

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

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

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

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

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

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

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

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

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

8
include/openmc/weight_windows.h
View file

@ -52,8 +52,12 @@ struct WeightWindow {
double weight_cutoff {DEFAULT_WEIGHT_CUTOFF};
int max_split {10};
//! Whether the weight window is in a valid state
bool is_valid() const { return lower_weight >= 0.0; }
//! Whether the weight window is in a valid state. A non-positive lower
//! bound indicates that no weight window information exists at this
//! location (generators mark such cells with -1, and a lower bound of zero
//! conventionally turns the weight window game off in a cell, as in MCNP
//! wwinp files), in which case no weight window game is played.
bool is_valid() const { return lower_weight > 0.0; }
//! Adjust the weight window by a constant factor
void scale(double factor)

230
openmc/dagmc.py
View file

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

7
openmc/data/angle_distribution.py
View file

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

25
openmc/data/decay.py
View file

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

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

@ -4,50 +4,76 @@ import numpy as np
import openmc.checkvalue as cv
_FILES = {
('icrp74', 'neutron'): Path('icrp74') / 'neutrons.txt',
('icrp74', 'photon'): Path('icrp74') / 'photons.txt',
('icrp116', 'electron'): Path('icrp116') / 'electrons.txt',
('icrp116', 'helium'): Path('icrp116') / 'helium_ions.txt',
('icrp116', 'mu-'): Path('icrp116') / 'negative_muons.txt',
('icrp116', 'pi-'): Path('icrp116') / 'negative_pions.txt',
('icrp116', 'neutron'): Path('icrp116') / 'neutrons.txt',
('icrp116', 'photon'): Path('icrp116') / 'photons.txt',
('icrp116', 'photon kerma'): Path('icrp116') / 'photons_kerma.txt',
('icrp116', 'mu+'): Path('icrp116') / 'positive_muons.txt',
('icrp116', 'pi+'): Path('icrp116') / 'positive_pions.txt',
('icrp116', 'positron'): Path('icrp116') / 'positrons.txt',
('icrp116', 'proton'): Path('icrp116') / 'protons.txt',
_FULL_GEOMETRIES = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO')
_LIMITED_GEOMETRIES = ('AP', 'PA', 'ISO')
_TABLES = {
('icrp74', 'effective', 'neutron'): (
Path('icrp74') / 'neutrons.txt', _FULL_GEOMETRIES),
('icrp74', 'effective', 'photon'): (
Path('icrp74') / 'photons.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'electron'): (
Path('icrp116') / 'electrons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'helium'): (
Path('icrp116') / 'helium_ions.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'mu-'): (
Path('icrp116') / 'negative_muons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'pi-'): (
Path('icrp116') / 'negative_pions.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'neutron'): (
Path('icrp116') / 'neutrons.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'photon'): (
Path('icrp116') / 'photons.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'photon kerma'): (
Path('icrp116') / 'photons_kerma.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'mu+'): (
Path('icrp116') / 'positive_muons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'pi+'): (
Path('icrp116') / 'positive_pions.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'positron'): (
Path('icrp116') / 'positrons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'proton'): (
Path('icrp116') / 'protons.txt', _FULL_GEOMETRIES),
('icrp74', 'ambient', 'neutron'): (
Path('icrp74') / 'neutrons_H10.txt', None),
('icrp74', 'ambient', 'photon'): (
Path('icrp74') / 'photons_H10.txt', None),
}
_DOSE_TABLES = {}
def _load_dose_icrp(data_source: str, particle: str):
"""Load effective dose tables from text files.
def _load_dose_table(data_source: str, dose_quantity: str, particle: str):
"""Load dose tables from text files.
Parameters
----------
data_source : {'icrp74', 'icrp116'}
The dose conversion data source to use
dose_quantity : {'effective', 'ambient'}
Dose quantity to load. 'ambient' corresponds to ambient dose
equivalent H*(10).
particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'}
Incident particle
"""
path = Path(__file__).parent / _FILES[data_source, particle]
key = (data_source, dose_quantity, particle)
path = Path(__file__).parent / _TABLES[key][0]
data = np.loadtxt(path, skiprows=3, encoding='utf-8')
data[:, 0] *= 1e6 # Change energies to eV
_DOSE_TABLES[data_source, particle] = data
_DOSE_TABLES[key] = data
def dose_coefficients(particle, geometry='AP', data_source='icrp116'):
"""Return effective dose conversion coefficients.
def dose_coefficients(
particle, geometry='AP', data_source='icrp116', dose_quantity='effective'
):
"""Return dose conversion coefficients.
This function provides fluence (and air kerma) to effective or ambient dose
(H*(10)) conversion coefficients for various types of external exposures
based on values in ICRP publications. Corrected values found in a
corrigendum are used rather than the values in the original report.
Available libraries include `ICRP Publication 74
This function provides fluence (and air kerma) to effective dose or ambient
dose equivalent (H*(10)) conversion coefficients for various types of
external exposures based on values in ICRP publications. Corrected values
found in a corrigendum are used rather than the values in the original
report. Available libraries include `ICRP Publication 74
<https://doi.org/10.1016/S0146-6453(96)90010-X>` and `ICRP Publication 116
<https://doi.org/10.1016/j.icrp.2011.10.001>`.
@ -63,45 +89,58 @@ def dose_coefficients(particle, geometry='AP', data_source='icrp116'):
particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'}
Incident particle
geometry : {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'}
Irradiation geometry assumed. Refer to ICRP-116 (Section 3.2) for the
meaning of the options here.
Irradiation geometry assumed for effective dose coefficients. Refer to
ICRP-116 (Section 3.2) for the meaning of the options here. This
argument does not apply when ``dose_quantity`` is 'ambient'.
data_source : {'icrp74', 'icrp116'}
The data source for the effective dose conversion coefficients.
The data source for the dose conversion coefficients.
dose_quantity : {'effective', 'ambient'}
Dose quantity to return. 'effective' returns effective dose
coefficients; 'ambient' returns ambient dose equivalent (H*(10))
coefficients.
Returns
-------
energy : numpy.ndarray
Energies at which dose conversion coefficients are given
dose_coeffs : numpy.ndarray
Effective dose coefficients in [pSv cm^2] at provided energies. For
'photon kerma', the coefficients are given in [Sv/Gy].
Dose coefficients in [pSv cm^2] at provided energies. For 'photon
kerma', the coefficients are given in [Sv/Gy].
"""
cv.check_value('geometry', geometry, {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'})
cv.check_value('geometry', geometry, _FULL_GEOMETRIES)
cv.check_value('data_source', data_source, {'icrp74', 'icrp116'})
cv.check_value('dose_quantity', dose_quantity, {'effective', 'ambient'})
if (data_source, particle) not in _FILES:
available_particles = sorted({p for (ds, p) in _FILES if ds == data_source})
key = (data_source, dose_quantity, particle)
if key not in _TABLES:
available_particles = sorted(
p for ds, dq, p in _TABLES
if ds == data_source and dq == dose_quantity
)
msg = (
f"'{particle}' has no dose data in data source {data_source}. "
f"Available particles for {data_source} are: {available_particles}"
f"'{particle}' has no {dose_quantity} dose data in data source "
f"{data_source}. Available particles for {data_source} "
f"with dose quantity {dose_quantity} are: {available_particles}"
)
raise ValueError(msg)
elif (data_source, particle) not in _DOSE_TABLES:
_load_dose_icrp(data_source, particle)
elif key not in _DOSE_TABLES:
_load_dose_table(data_source, dose_quantity, particle)
# Get all data for selected particle
data = _DOSE_TABLES[data_source, particle]
data = _DOSE_TABLES[key]
columns = _TABLES[key][1]
# Determine index for selected geometry
if particle in ('neutron', 'photon', 'proton', 'photon kerma'):
columns = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO')
if columns is None:
if geometry != 'AP':
raise ValueError(
"Irradiation geometry is not defined for ambient dose "
"equivalent coefficients. Use the default geometry='AP'."
)
index = 0
else:
columns = ('AP', 'PA', 'ISO')
index = columns.index(geometry)
index = columns.index(geometry)
# Pull out energy and dose from table
energy = data[:, 0].copy()
dose_coeffs = data[:, index + 1].copy()
return energy, dose_coeffs

50
openmc/data/dose/icrp74/neutrons_H10.txt Normal file
View file

@ -0,0 +1,50 @@
Neutrons: Ambient per fluence, in units of pSv cm², for monoenergetic particles incident.
Energy (MeV) Dose
1.00E-09 6.60
1.00E-08 9.00
2.53E-08 10.6
1.00E-07 12.9
2.00E-07 13.5
5.00E-07 13.6
1.00E-06 13.3
2.00E-06 12.9
5.00E-06 12.0
1.00E-05 11.3
2.00E-05 10.6
5.00E-05 9.90
1.00E-04 9.40
2.00E-04 8.90
5.00E-04 8.30
1.00E-03 7.90
2.00E-03 7.70
5.00E-03 8.00
1.00E-02 10.5
2.00E-02 16.6
3.00E-02 23.7
5.00E-02 41.1
7.00E-02 60.0
1.00E-01 88.0
1.50E-01 132
2.00E-01 170
3.00E-01 233
5.00E-01 322
7.00E-01 375
9.00E-01 400
1 416
1.2 425
2 420
3 412
4 408
5 405
6 400
7 405
8 409
9 420
10 440
12 480
14 520
15 540
16 555
18 570
20 600

28
openmc/data/dose/icrp74/photons_H10.txt Normal file
View file

@ -0,0 +1,28 @@
Photons: Ambient dose (H*10) per fluence, in units of pSv cm²
Energy (MeV) Dose
0.010 0.061
0.015 0.83
0.020 1.05
0.030 0.81
0.040 0.64
0.050 0.55
0.060 0.51
0.080 0.53
0.100 0.61
0.150 0.89
0.200 1.20
0.300 1.80
0.400 2.38
0.500 2.93
0.600 3.44
0.800 4.38
1 5.20
1.5 6.90
2 8.60
3 11.1
4 13.4
5 15.5
6 17.6
8 21.6
10 25.6

34
openmc/data/endf.py
View file

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

7
openmc/data/fission_energy.py
View file

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

10
openmc/data/neutron.py
View file

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

22
openmc/data/photon.py
View file

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

7
openmc/data/reaction.py
View file

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

7
openmc/data/resonance.py
View file

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

3
openmc/data/resonance_covariance.py
View file

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

7
openmc/data/thermal.py
View file

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

12
openmc/deplete/chain.py
View file

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

11
openmc/deplete/d1s.py
View file

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

4
openmc/deplete/independent_operator.py
View file

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

66
openmc/deplete/microxs.py
View file

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

2
openmc/deplete/r2s.py
View file

@ -347,7 +347,7 @@ class R2SManager:
# Run neutron transport and get fluxes and micros. Run via openmc.lib to
# maintain a consistent parallelism strategy with the activation step.
with TemporarySession():
with TemporarySession(output=False):
self.results['fluxes'], self.results['micros'] = get_microxs_and_flux(
self.neutron_model, domains, **micro_kwargs)

4
openmc/examples.py
View file

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

5
openmc/lib/core.py
View file

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

386
openmc/lib/plot.py
View file

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

13
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
@ -2287,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
@ -2322,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
@ -2332,6 +2339,10 @@ class Materials(cv.CheckedList):
for material, flux, energy in zip(
self, multigroup_fluxes, energy_group_structures
):
if material.volume is None:
raise ValueError(
f"Material {material.id} has no volume; cannot deplete"
)
temperature = material.temperature or 293.6
micro_xs = openmc.deplete.MicroXS.from_multigroup_flux(
energies=energy,

225
openmc/mesh.py
View file

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

8
openmc/mgxs/library.py
View file

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

218
openmc/model/model.py
View file

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

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)

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.

12
src/collision_track.cpp
View file

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

227
src/dagmc.cpp
View file

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

8
src/distribution_spatial.cpp
View file

@ -262,9 +262,11 @@ std::pair<Position, double> SphericalIndependent::sample(uint64_t* seed) const
MeshSpatial::MeshSpatial(pugi::xml_node node)
{
if (get_node_value(node, "type", true, true) != "mesh") {
fatal_error(fmt::format(
"Incorrect spatial type '{}' for a MeshSpatial distribution"));
auto spatial_type = get_node_value(node, "type", true, true);
if (spatial_type != "mesh") {
fatal_error(
fmt::format("Incorrect spatial type '{}' for a MeshSpatial distribution",
spatial_type));
}
// No in-tet distributions implemented, could include distributions for the

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

3
src/finalize.cpp
View file

@ -149,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;
@ -219,6 +221,7 @@ int openmc_reset()
settings::cmfd_run = false;
simulation::n_lost_particles = 0;
simulation::simulation_tracks_completed = 0;
return 0;
}

37
src/geometry.cpp
View file

@ -26,16 +26,22 @@ int n_coord_levels;
vector<int64_t> overlap_check_count;
vector<OverlapKey> overlap_keys;
std::unordered_map<OverlapKey, int, OverlapKeyHash> overlap_key_index;
} // namespace model
//==============================================================================
// Non-member functions
//==============================================================================
bool check_cell_overlap(GeometryState& p, bool error)
int check_cell_overlap(GeometryState& p, bool error)
{
int n_coord = p.n_coord();
// If no overlap found, return a nonphysical index
int overlap_index = -1;
// Loop through each coordinate level
for (int j = 0; j < n_coord; j++) {
Universe& univ = *model::universes[p.coord(j).universe()];
@ -44,21 +50,40 @@ bool check_cell_overlap(GeometryState& p, bool error)
for (auto index_cell : univ.cells_) {
Cell& c = *model::cells[index_cell];
if (c.contains(p.coord(j).r(), p.coord(j).u(), p.surface())) {
#pragma omp atomic
++model::overlap_check_count[index_cell];
if (index_cell != p.coord(j).cell()) {
if (error) {
fatal_error(
fmt::format("Overlapping cells detected: {}, {} on universe {}",
c.id_, model::cells[p.coord(j).cell()]->id_, univ.id_));
}
return true;
// With no fatal error (plotter is calling), now adds overlaps and
// ensures order does not matter when making overlap key
int cell_a = model::cells[index_cell]->id_;
int cell_b = model::cells[p.coord(j).cell()]->id_;
int a = std::min(cell_a, cell_b);
int b = std::max(cell_a, cell_b);
OverlapKey key {univ.id_, a, b};
#pragma omp critical(overlap_key_update)
{
auto it = model::overlap_key_index.find(key);
if (it != model::overlap_key_index.end()) {
overlap_index = it->second; // already exists, reuse index
} else {
int idx = int(model::overlap_keys.size());
model::overlap_keys.push_back(key);
model::overlap_key_index[key] = idx;
overlap_index = idx;
}
}
break;
}
#pragma omp atomic
++model::overlap_check_count[index_cell];
}
}
}
return false;
return overlap_index;
}
//==============================================================================

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

59
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;
@ -366,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;
@ -460,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();
@ -505,6 +550,8 @@ void read_separate_xml_files()
read_tallies_xml();
check_pulse_height_compatibility();
// Initialize distribcell_filters
prepare_distribcell();

33
src/lattice.cpp
View file

@ -10,6 +10,7 @@
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include "openmc/string_utils.h"
#include "openmc/vector.h"
#include "openmc/xml_interface.h"
@ -260,25 +261,33 @@ std::pair<double, array<int, 3>> RectLattice::distance(
// Determine the oncoming edge.
double x0 {copysign(0.5 * pitch_[0], u.x)};
double y0 {copysign(0.5 * pitch_[1], u.y)};
double z0;
double d = std::min(
u.x != 0.0 ? (x0 - x) / u.x : INFTY, u.y != 0.0 ? (y0 - y) / u.y : INFTY);
// Evaluate the distance to each oncoming edge independently. Comparing these
// distances directly (rather than reconstructing the crossing position)
// avoids the floating-point cancellation that occurs for large pitches.
double dx = u.x != 0.0 ? (x0 - x) / u.x : INFTY;
double dy = u.y != 0.0 ? (y0 - y) / u.y : INFTY;
double dz = INFTY;
if (is_3d_) {
z0 = copysign(0.5 * pitch_[2], u.z);
d = std::min(d, u.z != 0.0 ? (z0 - z) / u.z : INFTY);
double z0 {copysign(0.5 * pitch_[2], u.z)};
dz = u.z != 0.0 ? (z0 - z) / u.z : INFTY;
}
// Determine which lattice boundaries are being crossed
// The distance to the nearest lattice boundary is the smallest axial
// distance.
double d = std::min({dx, dy, dz});
// Determine which lattice boundaries are being crossed. The axis attaining
// the minimum is exactly equal to d, so at least one translation is always
// set for a finite crossing; a near-equal second axis indicates a corner
// crossing.
array<int, 3> lattice_trans = {0, 0, 0};
if (u.x != 0.0 && std::abs(x + u.x * d - x0) < FP_PRECISION)
if (isclose(d, dx, FP_COINCIDENT, FP_PRECISION))
lattice_trans[0] = copysign(1, u.x);
if (u.y != 0.0 && std::abs(y + u.y * d - y0) < FP_PRECISION)
if (isclose(d, dy, FP_COINCIDENT, FP_PRECISION))
lattice_trans[1] = copysign(1, u.y);
if (is_3d_) {
if (u.z != 0.0 && std::abs(z + u.z * d - z0) < FP_PRECISION)
lattice_trans[2] = copysign(1, u.z);
}
if (is_3d_ && isclose(d, dz, FP_COINCIDENT, FP_PRECISION))
lattice_trans[2] = copysign(1, u.z);
return {d, lattice_trans};
}

8
src/math_functions.cpp
View file

@ -959,4 +959,12 @@ void get_energy_index(
}
}
// Return true if two floating-point values are approximately equal within a
// combined relative and absolute tolerance.
bool isclose(double a, double b, double rel_tol, double abs_tol)
{
return std::abs(a - b) <=
std::max(rel_tol * std::max(std::abs(a), std::abs(b)), abs_tol);
}
} // namespace openmc

26
src/mesh.cpp
View file

@ -6,7 +6,8 @@
#define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers
#include <cmath> // for ceil
#include <cstddef> // for size_t
#include <numeric> // for accumulate
#include <limits>
#include <numeric> // for accumulate
#include <string>
#ifdef _MSC_VER
@ -1080,13 +1081,26 @@ int StructuredMesh::get_bin(Position r) const
int StructuredMesh::n_bins() const
{
return std::accumulate(
shape_.begin(), shape_.begin() + n_dimension_, 1, std::multiplies<>());
// Bin indices are stored as 32-bit ints in the tally system.
int64_t n = 1;
for (int i = 0; i < n_dimension_; ++i)
n *= shape_[i];
if (n > std::numeric_limits<int>::max()) {
fatal_error(fmt::format(
"Mesh {} has too many bins ({}) for 32-bit tally indexing", id_, n));
}
return static_cast<int>(n);
}
int StructuredMesh::n_surface_bins() const
{
return 4 * n_dimension_ * n_bins();
// Surface bin indices are stored as 32-bit ints in the tally system.
int64_t n = static_cast<int64_t>(n_bins()) * 4 * n_dimension_;
if (n > std::numeric_limits<int>::max()) {
fatal_error(fmt::format(
"Mesh {} has too many surface bins ({}) for tally indexing", id_, n));
}
return static_cast<int>(n);
}
tensor::Tensor<double> StructuredMesh::count_sites(
@ -3670,7 +3684,7 @@ Position LibMesh::sample_element(int32_t bin, uint64_t* seed) const
// Get tet vertex coordinates from LibMesh
std::array<Position, 4> tet_verts;
for (int i = 0; i < elem.n_nodes(); i++) {
auto node_ref = elem.node_ref(i);
const auto& node_ref = elem.node_ref(i);
tet_verts[i] = {node_ref(0), node_ref(1), node_ref(2)};
}
// Samples position within tet using Barycentric coordinates
@ -3700,7 +3714,7 @@ int LibMesh::n_vertices() const
Position LibMesh::vertex(int vertex_id) const
{
const auto node_ref = m_->node_ref(vertex_id);
const auto& node_ref = m_->node_ref(vertex_id);
if (length_multiplier_ > 0.0) {
return length_multiplier_ * Position(node_ref(0), node_ref(1), node_ref(2));
} else {

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

181
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()
@ -315,13 +338,14 @@ void Particle::event_cross_surface()
boundary().lattice_translation()[2] != 0) {
// Particle crosses lattice boundary
int i_lattice = coord(boundary().coord_level() - 1).lattice();
bool verbose = settings::verbosity >= 10 || trace();
cross_lattice(*this, boundary(), verbose);
event() = TallyEvent::LATTICE;
// Score cell to cell partial currents
if (!model::active_surface_tallies.empty()) {
auto& lat {*model::lattices[lowest_coord().lattice()]};
auto& lat {*model::lattices[i_lattice]};
bool is_valid;
Direction normal =
lat.get_normal(boundary().lattice_translation(), is_valid);
@ -449,61 +473,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 +550,7 @@ void Particle::event_death()
// Finish particle track output.
if (write_track()) {
write_particle_track(*this);
finalize_particle_track(*this);
}
@ -538,11 +574,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 +905,27 @@ void Particle::write_restart() const
write_dataset(file_id, "id", id());
write_dataset(file_id, "type", type().pdg_number());
// Get source site data for the particle that got lost
int64_t i = current_work();
SourceSite site;
if (settings::run_mode == RunMode::EIGENVALUE) {
// take source data from primary bank for eigenvalue simulation
write_dataset(file_id, "weight", simulation::source_bank[i - 1].wgt);
write_dataset(file_id, "energy", simulation::source_bank[i - 1].E);
write_dataset(file_id, "xyz", simulation::source_bank[i - 1].r);
write_dataset(file_id, "uvw", simulation::source_bank[i - 1].u);
write_dataset(file_id, "time", simulation::source_bank[i - 1].time);
site = simulation::source_bank[i];
} else if (settings::run_mode == RunMode::FIXED_SOURCE &&
settings::use_shared_secondary_bank &&
i < simulation::shared_secondary_bank_read.size()) {
site = simulation::shared_secondary_bank_read[i];
} else if (settings::run_mode == RunMode::FIXED_SOURCE) {
// re-sample using rng random number seed used to generate source particle
int64_t id = (simulation::total_gen + overall_generation() - 1) *
settings::n_particles +
simulation::work_index[mpi::rank] + i;
// Re-sample using the same seed used to generate the source particle.
// current_work() is 0-indexed, compute_particle_id expects 1-indexed.
int64_t id = compute_transport_seed(compute_particle_id(i + 1));
uint64_t seed = init_seed(id, STREAM_SOURCE);
// re-sample source site
auto site = sample_external_source(&seed);
write_dataset(file_id, "weight", site.wgt);
write_dataset(file_id, "energy", site.E);
write_dataset(file_id, "xyz", site.r);
write_dataset(file_id, "uvw", site.u);
write_dataset(file_id, "time", site.time);
site = sample_external_source(&seed);
}
write_dataset(file_id, "weight", site.wgt);
write_dataset(file_id, "energy", site.E);
write_dataset(file_id, "xyz", site.r);
write_dataset(file_id, "uvw", site.u);
write_dataset(file_id, "time", site.time);
// Close file
file_close(file_id);

23
src/particle_restart.cpp
View file

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

25
src/physics.cpp
View file

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

9
src/physics_mg.cpp
View file

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

234
src/plot.cpp
View file

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

20
src/random_ray/flat_source_domain.cpp
View file

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

87
src/random_ray/random_ray_simulation.cpp
View file

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

24
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};
@ -321,7 +322,7 @@ void get_run_parameters(pugi::xml_node node_base)
get_node_value_bool(random_ray_node, "volume_normalized_flux_tallies");
}
if (check_for_node(random_ray_node, "adjoint")) {
FlatSourceDomain::adjoint_ =
FlatSourceDomain::adjoint_requested_ =
get_node_value_bool(random_ray_node, "adjoint");
}
if (check_for_node(random_ray_node, "sample_method")) {
@ -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()

389
src/simulation.cpp
View file

@ -16,6 +16,7 @@
#include "openmc/particle.h"
#include "openmc/photon.h"
#include "openmc/random_lcg.h"
#include "openmc/random_ray/flat_source_domain.h"
#include "openmc/settings.h"
#include "openmc/source.h"
#include "openmc/state_point.h"
@ -86,7 +87,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();
@ -200,9 +201,12 @@ int openmc_simulation_finalize()
if (settings::output_tallies && mpi::master)
write_tallies();
// If weight window generators are present in this simulation,
// write a weight windows file
if (variance_reduction::weight_windows_generators.size() > 0) {
// If weight window generators are present in this simulation, write a
// weight windows file. This is skipped during the forward solve of an
// adjoint (FW-CADIS) run, where only the adjoint-derived weight windows
// are meaningful.
if (variance_reduction::weight_windows_generators.size() > 0 &&
FlatSourceDomain::solve_ != RandomRaySolve::FORWARD_FOR_ADJOINT) {
openmc_weight_windows_export();
}
@ -214,6 +218,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 +272,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 +350,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 +575,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 +599,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 +635,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 +651,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 +665,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 +697,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 +879,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 +889,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 +1037,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 +1046,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

9
src/source.cpp
View file

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

11
src/state_point.cpp
View file

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

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

26
src/tallies/tally.cpp
View file

@ -512,7 +512,7 @@ void Tally::set_strides()
// longest stride.
auto n = filters_.size();
strides_.resize(n, 0);
int stride = 1;
int64_t stride = 1;
for (int i = n - 1; i >= 0; --i) {
strides_[i] = stride;
stride *= model::tally_filters[filters_[i]]->n_bins();
@ -644,8 +644,24 @@ void Tally::set_scores(const vector<std::string>& scores)
break;
case HEATING:
if (settings::photon_transport)
estimator_ = TallyEstimator::COLLISION;
if (settings::photon_transport) {
// Photon heating requires a collision estimator (analog energy
// balance). However, if the tally only scores neutrons, we can keep the
// tracklength estimator since neutron heating uses kerma coefficients
// that support tracklength scoring.
bool neutron_only = false;
for (auto i_filt : filters_) {
auto pf =
dynamic_cast<ParticleFilter*>(model::tally_filters[i_filt].get());
if (pf && pf->particles().size() == 1 &&
pf->particles()[0].is_neutron()) {
neutron_only = true;
break;
}
}
if (!neutron_only)
estimator_ = TallyEstimator::COLLISION;
}
break;
case SCORE_PULSE_HEIGHT: {
@ -871,7 +887,7 @@ void Tally::accumulate()
if (higher_moments_) {
#pragma omp parallel for
// filter bins (specific cell, energy bins)
for (int i = 0; i < results_.shape(0); ++i) {
for (int64_t i = 0; i < results_.shape(0); ++i) {
// score bins (flux, total reaction rate, fission reaction rate, etc.)
for (int j = 0; j < results_.shape(1); ++j) {
double val = results_(i, j, TallyResult::VALUE) * norm;
@ -886,7 +902,7 @@ void Tally::accumulate()
} else {
#pragma omp parallel for
// filter bins (specific cell, energy bins)
for (int i = 0; i < results_.shape(0); ++i) {
for (int64_t i = 0; i < results_.shape(0); ++i) {
// score bins (flux, total reaction rate, fission reaction rate, etc.)
for (int j = 0; j < results_.shape(1); ++j) {
double val = results_(i, j, TallyResult::VALUE) * norm;

27
src/tallies/tally_scoring.cpp
View file

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

4
src/tallies/trigger.cpp
View file

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

7
src/weight_windows.cpp
View file

@ -791,7 +791,7 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
if (method_string == "magic") {
method_ = WeightWindowUpdateMethod::MAGIC;
if (settings::solver_type == SolverType::RANDOM_RAY &&
FlatSourceDomain::adjoint_) {
FlatSourceDomain::adjoint_requested_) {
fatal_error("Random ray weight window generation with MAGIC cannot be "
"done in adjoint mode.");
}
@ -800,7 +800,7 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
if (settings::solver_type != SolverType::RANDOM_RAY) {
fatal_error("FW-CADIS can only be run in random ray solver mode.");
}
FlatSourceDomain::adjoint_ = true;
FlatSourceDomain::adjoint_requested_ = true;
if (check_for_node(node, "targets")) {
FlatSourceDomain::fw_cadis_local_ = true;
targets_ = get_node_array<size_t>(node, "targets");
@ -837,7 +837,8 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
ratio_));
if (ratio_ <= 1.0)
fatal_error(fmt::format("Invalid weight window ratio '{}' (<= 1.0) "
"specified for weight window generation"));
"specified for weight window generation",
ratio_));
// create a matching weight windows object
auto wws = WeightWindows::create();

21
tests/cpp_unit_tests/test_math.cpp
View file

@ -355,3 +355,24 @@ TEST_CASE("Test broaden_wmp_polynomials")
REQUIRE_THAT(ref_val, Catch::Matchers::Approx(test_val));
}
}
TEST_CASE("Test isclose")
{
using openmc::isclose;
// Identical values are always close, regardless of tolerances.
REQUIRE(isclose(1.0, 1.0, 0.0, 0.0));
REQUIRE(isclose(0.0, 0.0, 0.0, 0.0));
// Absolute tolerance governs comparisons near zero.
REQUIRE(isclose(0.0, 1e-15, 0.0, 1e-14));
REQUIRE_FALSE(isclose(0.0, 1e-13, 0.0, 1e-14));
// Relative tolerance scales with the magnitude of the operands.
REQUIRE(isclose(1.0e12, 1.0e12 + 1.0, 1e-12, 0.0));
REQUIRE_FALSE(isclose(1.0e12, 1.0e12 + 10.0, 1e-12, 0.0));
// The looser of the two tolerances wins.
REQUIRE(isclose(1.0, 1.0 + 1e-13, 0.0, 1e-12));
REQUIRE(isclose(1.0e6, 1.0e6 + 1e-4, 1e-9, 0.0));
}

27
tests/cpp_unit_tests/test_tally.cpp
View file

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

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

Binary file not shown.

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

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

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

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

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

Binary file not shown.

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

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

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

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

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

Binary file not shown.

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

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

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