Merge branch 'develop' into virtual_lattice_0.15.2

This commit is contained in:
skywalker_cn 2025-12-18 13:04:39 +00:00
commit 5d3bc29074
527 changed files with 24388 additions and 14341 deletions

View file

@ -75,6 +75,7 @@ jobs:
LIBMESH: ${{ matrix.libmesh }}
NPY_DISABLE_CPU_FEATURES: "AVX512F AVX512_SKX"
OPENBLAS_NUM_THREADS: 1
PYTEST_ADDOPTS: --cov=openmc --cov-report=lcov:coverage-python.lcov
# libfabric complains about fork() as a result of using Python multiprocessing.
# We can work around it with RDMAV_FORK_SAFE=1 in libfabric < 1.13 and with
# FI_EFA_FORK_SAFE=1 in more recent versions.
@ -171,11 +172,37 @@ jobs:
uses: mxschmitt/action-tmate@v3
timeout-minutes: 10
- name: after_success
- name: Generate C++ coverage (gcovr)
shell: bash
run: |
cpp-coveralls -i src -i include -e src/external --exclude-pattern "/usr/*" --dump cpp_cov.json
coveralls --merge=cpp_cov.json --service=github
# Produce LCOV directly from gcov data in the build tree
gcovr \
--root "$GITHUB_WORKSPACE" \
--object-directory "$GITHUB_WORKSPACE/build" \
--filter "$GITHUB_WORKSPACE/src" \
--filter "$GITHUB_WORKSPACE/include" \
--exclude "$GITHUB_WORKSPACE/src/external/.*" \
--exclude "$GITHUB_WORKSPACE/src/include/openmc/external/.*" \
--gcov-ignore-errors source_not_found \
--gcov-ignore-errors output_error \
--gcov-ignore-parse-errors suspicious_hits.warn \
--print-summary \
--lcov -o coverage-cpp.lcov || true
- name: Merge C++ and Python coverage
shell: bash
run: |
# Merge C++ and Python LCOV into a single file for upload
cat coverage-cpp.lcov coverage-python.lcov > coverage.lcov
- name: Upload coverage to Coveralls
if: ${{ hashFiles('coverage.lcov') != '' }}
uses: coverallsapp/github-action@v2
with:
github-token: ${{ secrets.GITHUB_TOKEN }}
parallel: true
flag-name: C++ and Python
path-to-lcov: coverage.lcov
finish:
needs: main
@ -184,5 +211,5 @@ jobs:
- name: Coveralls Finished
uses: coverallsapp/github-action@v2
with:
github-token: ${{ secrets.github_token }}
github-token: ${{ secrets.GITHUB_TOKEN }}
parallel-finished: true

280
AGENTS.md Normal file
View file

@ -0,0 +1,280 @@
# OpenMC AI Coding Agent Instructions
## Project Overview
OpenMC is a Monte Carlo particle transport code for simulating nuclear reactors,
fusion devices, or other systems with neutron/photon radiation. It's a hybrid
C++17/Python codebase where:
- **C++ core** (`src/`, `include/openmc/`) handles the computationally intensive transport simulation
- **Python API** (`openmc/`) provides user-facing model building, post-processing, and depletion capabilities
- **C API bindings** (`openmc/lib/`) wrap the C++ library via ctypes for runtime control
## Architecture & Key Components
### C++ Component Structure
- **Global vectors of unique_ptrs**: Core objects like `model::cells`, `model::universes`, `nuclides` are stored as `vector<unique_ptr<T>>` in nested namespaces (`openmc::model`, `openmc::simulation`, `openmc::settings`, `openmc::data`)
- **Custom container types**: OpenMC provides its own `vector`, `array`, `unique_ptr`, and `make_unique` in the `openmc::` namespace (defined in `vector.h`, `array.h`, `memory.h`). These are currently typedefs to `std::` equivalents but may become custom implementations for accelerator support. Always use `openmc::vector`, not `std::vector`.
- **Geometry systems**:
- **CSG (default)**: Arbitrarily complex Constructive Solid Geometry using `Surface`, `Region`, `Cell`, `Universe`, `Lattice`
- **DAGMC**: CAD-based geometry via Direct Accelerated Geometry Monte Carlo (optional, requires `OPENMC_USE_DAGMC`)
- **Unstructured mesh**: libMesh-based geometry (optional, requires `OPENMC_USE_LIBMESH`)
- **Particle tracking**: `Particle` class with `GeometryState` manages particle transport through geometry
- **Tallies**: Score quantities during simulation via `Filter` and `Tally` objects
- **Random ray solver**: Alternative deterministic method in `src/random_ray/`
- **Optional features**: DAGMC (CAD geometry), libMesh (unstructured mesh), MPI, all controlled by `#ifdef OPENMC_MPI`, etc.
### Python Component Structure
- **ID management**: All geometry objects (Cell, Surface, Material, etc.) inherit from `IDManagerMixin` which auto-assigns unique integer IDs and tracks them via class-level `used_ids` and `next_id`
- **Input validation**: Extensive use of `openmc.checkvalue` module functions (`check_type`, `check_value`, `check_length`) for all setters
- **XML I/O**: Most classes implement `to_xml_element()` and `from_xml_element()` for serialization to OpenMC's XML input format
- **HDF5 output**: Post-simulation data in statepoint files read via `openmc.StatePoint`
- **Depletion**: `openmc.deplete` implements burnup via operator-splitting with various integrators (Predictor, CECM, etc.)
- **Nuclear Data**: `openmc.data` provides programmatic access to nuclear data files (ENDF, ACE, HDF5)
## Critical Build & Test Workflows
### Build Dependencies
- **C++17 compiler**: GCC, Clang, or Intel
- **CMake** (3.16+): Required for configuring and building the C++ library
- **HDF5**: Required for cross section data and output file formats
- **libpng**: Used for generating visualization when OpenMC is run in plotting mode
Without CMake and HDF5, OpenMC cannot be compiled.
### Building the C++ Library
```bash
# Configure with CMake (from build/ directory)
cmake .. -DOPENMC_USE_MPI=ON -DOPENMC_USE_OPENMP=ON -DCMAKE_BUILD_TYPE=RelWithDebInfo
# Available CMake options (all default OFF except OPENMC_USE_OPENMP and OPENMC_BUILD_TESTS):
# -DOPENMC_USE_OPENMP=ON/OFF # OpenMP parallelism
# -DOPENMC_USE_MPI=ON/OFF # MPI support
# -DOPENMC_USE_DAGMC=ON/OFF # CAD geometry support
# -DOPENMC_USE_LIBMESH=ON/OFF # Unstructured mesh
# -DOPENMC_ENABLE_PROFILE=ON/OFF # Profiling flags
# -DOPENMC_ENABLE_COVERAGE=ON/OFF # Coverage analysis
# Build
make -j
# C++ unit tests (uses Catch2)
ctest
```
### Python Development
```bash
# Install in development mode (requires building C++ library first)
pip install -e .
# Python tests (uses pytest)
pytest tests/unit_tests/ # Fast unit tests
pytest tests/regression_tests/ # Full regression suite (requires nuclear data)
```
### Nuclear Data Setup (CRITICAL for Running OpenMC)
Most tests require the NNDC HDF5 nuclear cross-section library.
**Important**: Check if `OPENMC_CROSS_SECTIONS` is already set in the user's
environment before downloading, as many users already have nuclear data
installed. Though do note that if this variable is present that it may point to
different cross section data and that the NNDC data is required for tests to
pass.
**If not already configured, download and setup:**
```bash
# Download NNDC HDF5 cross section library (~800 MB compressed)
wget -q -O - https://anl.box.com/shared/static/teaup95cqv8s9nn56hfn7ku8mmelr95p.xz | tar -C $HOME -xJ
# Set environment variable (add to ~/.bashrc or ~/.zshrc for persistence)
export OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml
```
**Alternative**: Use the provided download script (checks if data exists before downloading):
```bash
bash tools/ci/download-xs.sh # Downloads both NNDC HDF5 and ENDF/B-VII.1 data
```
Without this data, regression tests will fail with "No cross_sections.xml file
found" errors, or, in the case that alternative cross section data is configured
the tests will execute but will not pass. The `cross_sections.xml` file is an
index listing paths to individual HDF5 nuclear data files for each nuclide.
## Testing Expectations
### Environment Requirements
- **Data**: As described above, OpenMC's test suite requires OpenMC to be configured with NNDC data.
- **OpenMP Settings**: OpenMC's tests may fail is more than two OpenMP threads are used. The environment variable `OMP_NUM_THREADS=2` should be set to avoid sporadic test failures.
- **Executable configuration**: The OpenMC executable should compiled with debug symbols enabled.
### C++ Tests
Located in `tests/cpp_unit_tests/`, use Catch2 framework. Run via `ctest` after building with `-DOPENMC_BUILD_TESTS=ON`.
### Python Unit Tests
Located in `tests/unit_tests/`, these are fast, standalone tests that verify Python API functionality without running full simulations. Use standard pytest patterns:
**Categories**:
- **API validation**: Test object creation, property setters/getters, XML serialization (e.g., `test_material.py`, `test_cell.py`, `test_source.py`)
- **Data processing**: Test nuclear data handling, cross sections, depletion chains (e.g., `test_data_neutron.py`, `test_deplete_chain.py`)
- **Library bindings**: Test `openmc.lib` ctypes interface with `model.init_lib()`/`model.finalize_lib()` (e.g., `test_lib.py`)
- **Geometry operations**: Test bounding boxes, containment, lattice generation (e.g., `test_bounding_box.py`, `test_lattice.py`)
**Common patterns**:
- Use fixtures from `tests/unit_tests/conftest.py` (e.g., `uo2`, `water`, `sphere_model`)
- Test invalid inputs with `pytest.raises(ValueError)` or `pytest.raises(TypeError)`
- Use `run_in_tmpdir` fixture for tests that create files
- Tests with `openmc.lib` require calling `model.init_lib()` in try/finally with `model.finalize_lib()`
**Example**:
```python
def test_material_properties():
m = openmc.Material()
m.add_nuclide('U235', 1.0)
assert 'U235' in m.nuclides
with pytest.raises(TypeError):
m.add_nuclide('H1', '1.0') # Invalid type
```
Unit tests should be fast. For tests requiring simulation output, use regression tests instead.
### Python Regression Tests
Regression tests compare OpenMC output against reference data. **Prefer using existing models from `openmc.examples` or those found in tests/unit_tests/conftest.py** (like `pwr_pin_cell()`, `pwr_assembly()`, `slab_mg()`) rather than building from scratch.
**Test Harness Types** (in `tests/testing_harness.py`):
- **PyAPITestHarness**: Standard harness for Python API tests. Compares `inputs_true.dat` (XML hash) and `results_true.dat` (statepoint k-eff and tally values). Requires `model.xml` generation.
- **HashedPyAPITestHarness**: Like PyAPITestHarness but hashes the results for compact comparison
- **TolerantPyAPITestHarness**: For tests with floating-point non-associativity (e.g., random ray solver with single precision). Uses relative tolerance comparisons.
- **WeightWindowPyAPITestHarness**: Compares weight window bounds from `weight_windows.h5`
- **CollisionTrackTestHarness**: Compares collision track data from `collision_track.h5` against `collision_track_true.h5`
- **TestHarness**: Base harness for XML-based tests (no Python model building)
- **PlotTestHarness**: Compares plot output files (PNG or voxel HDF5)
- **CMFDTestHarness**: Specialized for CMFD acceleration tests
- **ParticleRestartTestHarness**: Tests particle restart functionality
Almost all cases use either `PyAPITestHarness` or `HashedPyAPITestHarness`
**Example Test**:
```python
from openmc.examples import pwr_pin_cell
from tests.testing_harness import PyAPITestHarness
def test_my_feature():
model = pwr_pin_cell()
model.settings.particles = 1000 # Modify to exercise feature
harness = PyAPITestHarness('statepoint.10.h5', model)
harness.main()
```
**Workflow**: Create `test.py` and `__init__.py` in `tests/regression_tests/my_test/`, run `pytest --update` to generate reference files (`inputs_true.dat`, `results_true.dat`, etc.), then verify with `pytest` without `--update`. Test results should be generated with a debug build (`-DCMAKE_BUILD_TYPE=Debug`)
**Critical**: When modifying OpenMC code, regenerate affected test references with `pytest --update` and commit updated reference files.
### Test Configuration
`pytest.ini` sets: `python_files = test*.py`, `python_classes = NoThanks` (disables class-based test collection).
### Testing Options
For builds of OpenMC with MPI enabled, the `--mpi` flag should be passed to the test suite to ensure that appropriate tests are executed using two MPI processes.
The entire test suite can be executed with OpenMC running in event-based mode (instead of the default history-based mode) by providing the `--event` flag to the `pytest` command.
## Cross-Language Boundaries
The C API (defined in `include/openmc/capi.h`) exposes C++ functionality to Python via ctypes bindings in `openmc/lib/`. Example:
```cpp
// C++ API in capi.h
extern "C" int openmc_run();
// Python binding in openmc/lib/core.py
_dll.openmc_run.restype = c_int
def run():
_dll.openmc_run()
```
When modifying C++ public APIs, update corresponding ctypes signatures in `openmc/lib/*.py`.
## Code Style & Conventions
### C++ Style (enforced by .clang-format)
OpenMC generally tries to follow C++ core guidelines where possible
(https://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines) and follow
modern C++ practices (e.g. RAII) whenever possible.
- **Naming**:
- Classes: `CamelCase` (e.g., `HexLattice`)
- Functions/methods: `snake_case` (e.g., `get_indices`)
- Variables: `snake_case` with trailing underscore for class members (e.g., `n_particles_`, `energy_`)
- Constants: `UPPER_SNAKE_CASE` (e.g., `SQRT_PI`)
- **Namespaces**: All code in `openmc::` namespace, global state in sub-namespaces
- **Include order**: Related header first, then C/C++ stdlib, third-party libs, local headers
- **Comments**: C++-style (`//`) only, never C-style (`/* */`)
- **Standard**: C++17 features allowed
- **Formatting**: Run `clang-format` (version 15) before committing; install via `tools/dev/install-commit-hooks.sh`
### Python Style
- **PEP8** compliant
- **Docstrings**: numpydoc format for all public functions/methods
- **Type hints**: Use sparingly, primarily for complex signatures
- **Path handling**: Use `pathlib.Path` for filesystem operations, accept `str | os.PathLike` in function arguments
- **Dependencies**: Core dependencies only (numpy, scipy, h5py, pandas, matplotlib, lxml, ipython, uncertainties, setuptools, endf). Other packages must be optional
- **Python version**: Minimum 3.11 (as of Nov 2025)
### ID Management Pattern (Python)
When creating geometry objects, IDs can be auto-assigned or explicit:
```python
# Auto-assigned ID
cell = openmc.Cell() # Gets next available ID
# Explicit ID
cell = openmc.Cell(id=10) # Warning if ID already used
# Reset all IDs (useful in test fixtures)
openmc.reset_auto_ids()
```
### Input Validation Pattern (Python)
All setters use checkvalue functions:
```python
import openmc.checkvalue as cv
@property
def temperature(self):
return self._temperature
@temperature.setter
def temperature(self, temp):
cv.check_type('temperature', temp, Real)
cv.check_greater_than('temperature', temp, 0.0)
self._temperature = temp
```
### Working with HDF5 Files
C++ uses custom HDF5 wrappers in `src/hdf5_interface.cpp`. Python uses h5py directly. Statepoint format version is `VERSION_STATEPOINT` in `include/openmc/constants.h`.
### Conditional Compilation
Check for optional features:
```cpp
#ifdef OPENMC_MPI
// MPI-specific code
#endif
#ifdef OPENMC_DAGMC
// DAGMC-specific code
#endif
```
## Documentation
- **User docs**: Sphinx documentation in `docs/source/` hosted at https://docs.openmc.org
- **C++ docs**: Doxygen-style comments with `\brief`, `\param` tags
- **Python docs**: numpydoc format docstrings
## Common Pitfalls
1. **Forgetting nuclear data**: Tests fail without `OPENMC_CROSS_SECTIONS` environment variable
2. **ID conflicts**: Python objects with duplicate IDs trigger `IDWarning`, use `reset_auto_ids()` between tests
3. **MPI builds**: Code must work with and without MPI; use `#ifdef OPENMC_MPI` guards
4. **Path handling**: Use `pathlib.Path` in new Python code, not `os.path`
5. **Clang-format version**: CI uses version 15; other versions may produce different formatting

View file

@ -1,9 +1,43 @@
cff-version: 1.2.0
message: "If you use this software, please cite it as below."
title: OpenMC
authors:
- family-names: Romano
given-names: Paul K.
orcid: "https://orcid.org/0000-0002-1147-045X"
- family-names: Shriwise
given-names: Patrick C.
orcid: "https://orcid.org/0000-0002-3979-7665"
- family-names: Shimwell
given-names: Jonathan
orcid: "https://orcid.org/0000-0001-6909-0946"
- family-names: Harper
given-names: Sterling
- family-names: Boyd
given-names: Will
- family-names: Nelson
given-names: Adam G.
orcid: "https://orcid.org/0000-0002-3614-0676"
- family-names: Tramm
given-names: John R.
orcid: "https://orcid.org/0000-0002-5397-4402"
- family-names: Ridley
given-names: Gavin
orcid: "https://orcid.org/0000-0003-1635-8042"
- family-names: Johnson
given-names: Andrew
orcid: "https://orcid.org/0000-0003-2125-8775"
- family-names: Peterson
given-names: Ethan E.
orcid: "https://orcid.org/0000-0002-5694-7194"
- family-names: Herman
given-names: Bryan R.
preferred-citation:
authors:
- family-names: Romano
given-names: Paul K.
orcid: "https://orcid.org/0000-0002-1147-045X"
- final-names: Horelik
- family-names: Horelik
given-names: Nicholas E.
- family-names: Herman
given-names: Bryan R.

View file

@ -338,6 +338,7 @@ list(APPEND libopenmc_SOURCES
src/cell.cpp
src/chain.cpp
src/cmfd_solver.cpp
src/collision_track.cpp
src/cross_sections.cpp
src/dagmc.cpp
src/distribution.cpp

View file

@ -1,9 +1,12 @@
get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY)
find_package(fmt REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../fmt)
find_package(pugixml REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../pugixml)
find_package(xtl REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtl)
find_package(xtensor REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtensor)
# Compute the install prefix from this file's location
get_filename_component(_OPENMC_PREFIX "${OpenMC_CMAKE_DIR}/../../.." ABSOLUTE)
find_package(fmt CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
find_package(pugixml CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
find_package(xtl CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
find_package(xtensor CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
if(@OPENMC_USE_DAGMC@)
find_package(DAGMC REQUIRED HINTS @DAGMC_DIR@)
endif()

View file

@ -84,6 +84,17 @@ Functions
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_cell_get_density(int32_t index, const int32_t* instance, double* density)
Get the density of a cell
:param int32_t index: Index in the cells array
:param int32_t* instance: Which instance of the cell. If a null pointer is passed, the density
multiplier of the first instance is returned.
:param double* density: Density of the cell in [g/cm3]
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices)
Set the fill for a cell
@ -113,8 +124,22 @@ Functions
:param double T: Temperature in Kelvin
:param instance: Which instance of the cell. To set the temperature for all
instances, pass a null pointer.
:param set_contained: If the cell is not filled by a material, whether to set the temperatures
of all filled cells
:param bool set_contained: If the cell is not filled by a material, whether
to set the temperatures of all filled cells
:type instance: const int32_t*
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_cell_set_density(index index, double density, const int32_t* instance, bool set_contained)
Set the density of a cell.
:param int32_t index: Index in the cells array
:param double density: Density of the cell in [g/cm3]
:param instance: Which instance of the cell. To set the density multiplier for all
instances, pass a null pointer.
:param bool set_contained: If the cell is not filled by a material, whether
to set the density multiplier of all filled cells
:type instance: const int32_t*
:return: Return status (negative if an error occurred)
:rtype: int

View file

@ -121,9 +121,7 @@ pygments_style = 'tango'
# -- Options for HTML output ---------------------------------------------------
# The theme to use for HTML and HTML Help pages
import sphinx_rtd_theme
html_theme = 'sphinx_rtd_theme'
html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
html_baseurl = "https://docs.openmc.org/en/stable/"
html_logo = '_images/openmc_logo.png'

View file

@ -21,8 +21,8 @@ C++ code in OpenMC must conform to the most recent C++ standard that is fully
supported in the `version of the gcc compiler
<https://gcc.gnu.org/projects/cxx-status.html>`_ that is distributed with the
oldest version of Ubuntu that is still within its `standard support period
<https://ubuntu.com/about/release-cycle>`_. Ubuntu 20.04 LTS will be supported
through April 2025 and is distributed with gcc 9.3.0, which fully supports the
<https://ubuntu.com/about/release-cycle>`_. Ubuntu 22.04 LTS will be supported
through April 2027 and is distributed with gcc 11.4.0, which fully supports the
C++17 standard.
--------------------
@ -31,5 +31,5 @@ CMake Version Policy
Similar to the C++ standard policy, the minimum supported version of CMake
corresponds to whatever version is distributed with the oldest version of Ubuntu
still within its standard support period. Ubuntu 20.04 LTS is distributed with
CMake 3.16.
still within its standard support period. Ubuntu 22.04 LTS is distributed with
CMake 3.22.

View file

@ -0,0 +1,46 @@
.. _io_collision_track:
===========================
Collision Track File Format
===========================
When collision tracking is enabled with ``mcpl=false`` (the default), OpenMC
writes binary data to an HDF5 file named ``collision_track.h5``. The same data
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.0.
**/**
:Attributes:
- **filetype** (*char[]*) -- String indicating the type of file.
For collision-track files the value is ``"collision_track"``.
:Datasets:
- **collision_track_bank** (Compound type) -- Collision information
for each stored event. Each entry in the dataset corresponds to one
collision and contains the following fields:
- ``r`` (*double[3]*) -- Position of the collision in [cm].
- ``u`` (*double[3]*) -- Direction unit vector immediately after the collision.
- ``E`` (*double*) -- Incident particle energy before the collision in [eV].
- ``dE`` (*double*) -- Energy loss over the collision (:math:`E_\text{before} - E_\text{after}`) in [eV].
- ``time`` (*double*) -- Time of the collision in [s].
- ``wgt`` (*double*) -- Particle weight at the collision.
- ``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).
- ``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.
- ``particle`` (*int*) -- Particle type (0=neutron, 1=photon, 2=electron, 3=positron).
- ``parent_id`` (*int64*) -- Unique ID of the parent particle.
- ``progeny_id`` (*int64*) -- Progeny ID of the particle.
In an MPI run, OpenMC writes the combined dataset by gathering collision-track
entries from all ranks before flushing them to disk, so the final file appears
as though it were produced serially.

View file

@ -56,6 +56,27 @@ attributes:
.. _io_chain_reaction:
--------------------
``<source>`` Element
--------------------
The ``<source>`` element represents photon and electron sources associated with
the decay of a nuclide and contains information to construct an
:class:`openmc.stats.Univariate` object that represents this emission as an
energy distribution. This element has the following attributes:
:type:
The type of :class:`openmc.stats.Univariate` source term.
:particle:
The type of particle emitted, e.g., 'photon' or 'electron'
:parameters:
The parameters of the source term, e.g., for a
:class:`openmc.stats.Discrete` source, the energies (in [eV]) at which the
particles are emitted and their relative intensities in [Bq/atom] (in other
words, decay constants).
----------------------
``<reaction>`` Element
----------------------

View file

@ -4,7 +4,7 @@
Depletion Results File Format
=============================
The current version of the depletion results file format is 1.1.
The current version of the depletion results file format is 1.2.
**/**
@ -12,22 +12,20 @@ The current version of the depletion results file format is 1.1.
- **version** (*int[2]*) -- Major and minor version of the
statepoint file format.
:Datasets: - **eigenvalues** (*double[][][2]*) -- k-eigenvalues at each
time/stage. This array has shape (number of timesteps, number of
stages, value). The last axis contains the eigenvalue and the
associated uncertainty
- **number** (*double[][][][]*) -- Total number of atoms. This array
has shape (number of timesteps, number of stages, number of
:Datasets: - **eigenvalues** (*double[][2]*) -- k-eigenvalues at each timestep.
This array has shape (number of timesteps, 2). The second axis
contains the eigenvalue and its associated uncertainty.
- **number** (*double[][][]*) -- Total number of atoms at each
timestep. This array has shape (number of timesteps, number of
materials, number of nuclides).
- **reaction rates** (*double[][][][][]*) -- Reaction rates used to
build depletion matrices. This array has shape (number of
timesteps, number of stages, number of materials, number of
nuclides, number of reactions).
- **reaction rates** (*double[][][][]*) -- Reaction rates at each
timestep. This array has shape (number of timesteps, number of
materials, number of nuclides, number of reactions). Only stored if
write_rates=True.
- **time** (*double[][2]*) -- Time in [s] at beginning/end of each
step.
- **source_rate** (*double[][]*) -- Power in [W] or source rate in
[neutron/sec]. This array has shape (number of timesteps, number
of stages).
- **source_rate** (*double[]*) -- Power in [W] or source rate in
[neutron/sec] for each timestep.
- **depletion time** (*double[]*) -- Average process time in [s]
spent depleting a material across all burnable materials and,
if applicable, MPI processes.

View file

@ -38,11 +38,9 @@ Each ``<surface>`` element can have the following attributes or sub-elements:
:boundary:
The boundary condition for the surface. This can be "transmission",
"vacuum", "reflective", or "periodic". Periodic boundary conditions can
only be applied to x-, y-, and z-planes. Only axis-aligned periodicity is
supported, i.e., x-planes can only be paired with x-planes. Specify which
planes are periodic and the code will automatically identify which planes
are paired together.
"vacuum", "reflective", or "periodic". Specify which planes are
periodic and the code will automatically identify which planes are
paired together.
*Default*: "transmission"

View file

@ -44,6 +44,7 @@ Output Files
statepoint
source
collision_track
summary
properties
depletion_results

View file

@ -4,7 +4,7 @@
Properties File Format
======================
The current version of the properties file format is 1.0.
The current version of the properties file format is 1.1.
**/**
@ -25,6 +25,7 @@ The current version of the properties file format is 1.0.
**/geometry/cells/cell <uid>/**
:Datasets: - **temperature** (*double[]*) -- Temperature of the cell in [K].
- **density** (*double[]*) -- Density of the cell in [g/cm3].
**/materials/**

View file

@ -20,6 +20,85 @@ source neutrons.
*Default*: None
-----------------------------
``<collision_track>`` Element
-----------------------------
The ``<collision_track>`` element indicates to track information about particle
collisions based on a set of criteria and store these events in a file named
``collision_track.h5``. This file records details such as the position of the
interaction, direction of the incoming particle, incident energy and deposited
energy, weight, time of the interaction, and the delayed neutron group (0 for
prompt neutrons). Additional information such as the cell ID, material ID,
universe ID, nuclide ZAID, particle type, and event MT number are also stored.
Users can specify one or more criterion to filter collisions. If no criteria are
specified, it defaults to tracking all collisions across the model.
.. warning::
Storing all collisions can be very memory intensive. For more targeted
tracking, users can employ a variety of parameters such as ``cell_ids``,
``reactions``, ``universe_ids``, ``material_ids``, ``nuclides``, and
``deposited_E_threshold`` to refine the selection of particle interactions
to be banked.
This element can contain one or more of the following attributes or
sub-elements:
:max_collisions:
An integer indicating the maximum number of collisions to be banked per file.
*Default*: 1000
:max_collision_track_files:
An integer indicating the number of collision_track files to be used.
*Default*: 1
:mcpl:
An optional boolean to enable MCPL_-format instead of the native HDF5-based
format. If activated, the output file name and type is changed to
``collision_track.mcpl``.
*Default*: false
.. _MCPL: https://mctools.github.io/mcpl/mcpl.pdf
:cell_ids:
A list of integers representing cell IDs to define specific cells in which
collisions are to be banked.
*Default*: None
:universe_ids:
A list of integers representing the universe IDs to define specific
universes in which collisions are to be banked.
*Default*: None
:material_ids:
A list of integers representing the material IDs to define specific
materials in which collisions are to be banked.
*Default*: None
:nuclides:
A list of strings representing the nuclide, to define specific
define specific target nuclide collisions to be banked.
*Default*: None
:reactions:
A list of integers representing the ENDF-6 format MT numbers or strings
(e.g. (n,fission)) to define specific reaction types to be banked.
*Default*: None
:deposited_E_threshold:
A float defining the minimum deposited energy per collision (in eV) to
trigger banking.
*Default*: 0.0
----------------------------------
``<confidence_intervals>`` Element
----------------------------------
@ -178,6 +257,16 @@ history-based parallelism.
*Default*: false
--------------------------------
``<free_gas_threshold>`` Element
--------------------------------
The ``<free_gas_threshold>`` element specifies the energy multiplier, expressed
in units of :math:`kT`, that determines when the free gas scattering approach is
used for elastic scattering. Values must be positive.
*Default*: 400.0
-----------------------------------
``<generations_per_batch>`` Element
-----------------------------------

View file

@ -149,6 +149,8 @@ The current version of the statepoint file format is 18.1.
tallies will have a value of 0 unless otherwise instructed.
- **multiply_density** (*int*) -- Flag indicating whether reaction
rates should be multiplied by atom density (1) or not (0).
- **higher_moments** (*int*) -- Flag indicating whether
higher-order tally moments are enabled (1) or not (0).
:Datasets: - **n_realizations** (*int*) -- Number of realizations.
- **n_filters** (*int*) -- Number of filters used.

View file

@ -4,7 +4,7 @@
Summary File Format
===================
The current version of the summary file format is 6.0.
The current version of the summary file format is 6.1.
**/**
@ -38,6 +38,7 @@ The current version of the summary file format is 6.0.
is an array if the cell uses distributed materials, otherwise it is
a scalar.
- **temperature** (*double[]*) -- Temperature of the cell in Kelvin.
- **density** (*double[]*) -- Density of the cell in [g/cm3].
- **translation** (*double[3]*) -- Translation applied to the fill
universe. This dataset is present only if fill_type is set to
'universe'.

View file

@ -0,0 +1,362 @@
.. _methods_charged_particle_physics:
========================
Charged Particle Physics
========================
OpenMC neglects the spatial transport of charged particles (electrons and
positrons), assuming they deposit all their energy locally and produce
bremsstrahlung photons at their birth location. This approximation, called
thick-target bremsstrahlung (TTB) approximation is justified by the fact that
charged particles have much shorter stopping ranges compared to neutrons and
photons, especially in high-density materials.
-----------------------------
Charged Particle Interactions
-----------------------------
Bremsstrahlung
--------------
When a charged particle is decelerated in the field of an atom, some of its
kinetic energy is converted into electromagnetic radiation known as
bremsstrahlung, or 'braking radiation'. In each event, an electron or positron
with kinetic energy :math:`T` generates a photon with an energy :math:`E`
between :math:`0` and :math:`T`. Bremsstrahlung is described by a cross section
that is differential in photon energy, in the direction of the emitted photon,
and in the final direction of the charged particle. However, in Monte Carlo
simulations it is typical to integrate over the angular variables to obtain a
single differential cross section with respect to photon energy, which is often
expressed in the form
.. math::
:label: bremsstrahlung-dcs
\frac{d\sigma_{\text{br}}}{dE} = \frac{Z^2}{\beta^2} \frac{1}{E}
\chi(Z, T, \kappa),
where :math:`\kappa = E/T` is the reduced photon energy and :math:`\chi(Z, T,
\kappa)` is the scaled bremsstrahlung cross section, which is experimentally
measured.
Because electrons are attracted to atomic nuclei whereas positrons are
repulsed, the cross section for positrons is smaller, though it approaches that
of electrons in the high energy limit. To obtain the positron cross section, we
multiply :eq:`bremsstrahlung-dcs` by the :math:`\kappa`-independent factor used
in Salvat_,
.. math::
:label: positron-factor
\begin{aligned}
F_{\text{p}}(Z,T) =
& 1 - \text{exp}(-1.2359\times 10^{-1}t + 6.1274\times 10^{-2}t^2 - 3.1516\times 10^{-2}t^3 \\
& + 7.7446\times 10^{-3}t^4 - 1.0595\times 10^{-3}t^5 + 7.0568\times 10^{-5}t^6 \\
& - 1.8080\times 10^{-6}t^7),
\end{aligned}
where
.. math::
:label: positron-factor-t
t = \ln\left(1 + \frac{10^6}{Z^2}\frac{T}{\text{m}_\text{e}c^2} \right).
:math:`F_{\text{p}}(Z,T)` is the ratio of the radiative stopping powers for
positrons and electrons. Stopping power describes the average energy loss per
unit path length of a charged particle as it passes through matter:
.. math::
:label: stopping-power
-\frac{dT}{ds} = n \int E \frac{d\sigma}{dE} dE \equiv S(T),
where :math:`n` is the number density of the material and :math:`d\sigma/dE` is
the cross section differential in energy loss. The total stopping power
:math:`S(T)` can be separated into two components: the radiative stopping
power :math:`S_{\text{rad}}(T)`, which refers to energy loss due to
bremsstrahlung, and the collision stopping power :math:`S_{\text{col}}(T)`,
which refers to the energy loss due to inelastic collisions with bound
electrons in the material that result in ionization and excitation. The
radiative stopping power for electrons is given by
.. math::
:label: radiative-stopping-power
S_{\text{rad}}(T) = n \frac{Z^2}{\beta^2} T \int_0^1 \chi(Z,T,\kappa)
d\kappa.
To obtain the radiative stopping power for positrons,
:eq:`radiative-stopping-power` is multiplied by :eq:`positron-factor`.
While the models for photon interactions with matter described above can safely
assume interactions occur with free atoms, sampling the target atom based on
the macroscopic cross sections, molecular effects cannot necessarily be
disregarded for charged particle treatment. For compounds and mixtures, the
bremsstrahlung cross section is calculated using Bragg's additivity rule as
.. math::
:label: material-bremsstrahlung-dcs
\frac{d\sigma_{\text{br}}}{dE} = \frac{1}{\beta^2 E} \sum_i \gamma_i Z^2_i
\chi(Z_i, T, \kappa),
where the sum is over the constituent elements and :math:`\gamma_i` is the
atomic fraction of the :math:`i`-th element. Similarly, the radiative stopping
power is calculated using Bragg's additivity rule as
.. math::
:label: material-radiative-stopping-power
S_{\text{rad}}(T) = \sum_i w_i S_{\text{rad},i}(T),
where :math:`w_i` is the mass fraction of the :math:`i`-th element and
:math:`S_{\text{rad},i}(T)` is found for element :math:`i` using
:eq:`radiative-stopping-power`. The collision stopping power, however, is a
function of certain quantities such as the mean excitation energy :math:`I` and
the density effect correction :math:`\delta_F` that depend on molecular
properties. These quantities cannot simply be summed over constituent elements
in a compound, but should instead be calculated for the material. The Bethe
formula can be used to find the collision stopping power of the material:
.. math::
:label: material-collision-stopping-power
S_{\text{col}}(T) = \frac{2 \pi r_e^2 m_e c^2}{\beta^2} N_A \frac{Z}{A_M}
[\ln(T^2/I^2) + \ln(1 + \tau/2) + F(\tau) - \delta_F(T)],
where :math:`N_A` is Avogadro's number, :math:`A_M` is the molar mass,
:math:`\tau = T/m_e`, and :math:`F(\tau)` depends on the particle type. For
electrons,
.. math::
:label: F-electron
F_{-}(\tau) = (1 - \beta^2)[1 + \tau^2/8 - (2\tau + 1) \ln2],
while for positrons
.. math::
:label: F-positron
F_{+}(\tau) = 2\ln2 - (\beta^2/12)[23 + 14/(\tau + 2) + 10/(\tau + 2)^2 +
4/(\tau + 2)^3].
The density effect correction :math:`\delta_F` takes into account the reduction
of the collision stopping power due to the polarization of the material the
charged particle is passing through by the electric field of the particle.
It can be evaluated using the method described by Sternheimer_, where the
equation for :math:`\delta_F` is
.. math::
:label: density-effect-correction
\delta_F(\beta) = \sum_{i=1}^n f_i \ln[(l_i^2 + l^2)/l_i^2] -
l^2(1-\beta^2).
Here, :math:`f_i` is the oscillator strength of the :math:`i`-th transition,
given by :math:`f_i = n_i/Z`, where :math:`n_i` is the number of electrons in
the :math:`i`-th subshell. The frequency :math:`l` is the solution of the
equation
.. math::
:label: density-effect-l
\frac{1}{\beta^2} - 1 = \sum_{i=1}^{n} \frac{f_i}{\bar{\nu}_i^2 + l^2},
where :math:`\bar{v}_i` is defined as
.. math::
:label: density-effect-nubar
\bar{\nu}_i = h\nu_i \rho / h\nu_p.
The plasma energy :math:`h\nu_p` of the medium is given by
.. math::
:label: plasma-frequency
h\nu_p = \sqrt{\frac{(hc)^2 r_e \rho_m N_A Z}{\pi A}},
where :math:`A` is the atomic weight and :math:`\rho_m` is the density of the
material. In :eq:`density-effect-nubar`, :math:`h\nu_i` is the oscillator
energy, and :math:`\rho` is an adjustment factor introduced to give agreement
between the experimental values of the oscillator energies and the mean
excitation energy. The :math:`l_i` in :eq:`density-effect-correction` are
defined as
.. math::
:label: density-effect-li
\begin{aligned}
l_i &= (\bar{\nu}_i^2 + 2/3f_i)^{1/2} ~~~~&\text{for}~~ \bar{\nu}_i > 0 \\
l_n &= f_n^{1/2} ~~~~&\text{for}~~ \bar{\nu}_n = 0,
\end{aligned}
where the second case applies to conduction electrons. For a conductor,
:math:`f_n` is given by :math:`n_c/Z`, where :math:`n_c` is the effective
number of conduction electrons, and :math:`v_n = 0`. The adjustment factor
:math:`\rho` is determined using the equation for the mean excitation energy:
.. math::
:label: mean-excitation-energy
\ln I = \sum_{i=1}^{n-1} f_i \ln[(h\nu_i\rho)^2 + 2/3f_i(h\nu_p)^2]^{1/2} +
f_n \ln (h\nu_pf_n^{1/2}).
.. _ttb:
Thick-Target Bremsstrahlung Approximation
+++++++++++++++++++++++++++++++++++++++++
Since charged particles lose their energy on a much shorter distance scale than
neutral particles, not much error should be introduced by neglecting to
transport electrons. However, the bremsstrahlung emitted from high energy
electrons and positrons can travel far from the interaction site. Thus, even
without a full electron transport mode it is necessary to model bremsstrahlung.
We use a thick-target bremsstrahlung (TTB) approximation based on the models in
Salvat_ and Kaltiaisenaho_ for generating bremsstrahlung photons, which assumes
the charged particle loses all its energy in a single homogeneous material
region.
To model bremsstrahlung using the TTB approximation, we need to know the number
of photons emitted by the charged particle and the energy distribution of the
photons. These quantities can be calculated using the continuous slowing down
approximation (CSDA). The CSDA assumes charged particles lose energy
continuously along their trajectory with a rate of energy loss equal to the
total stopping power, ignoring fluctuations in the energy loss. The
approximation is useful for expressing average quantities that describe how
charged particles slow down in matter. For example, the CSDA range approximates
the average path length a charged particle travels as it slows to rest:
.. math::
:label: csda-range
R(T) = \int^T_0 \frac{dT'}{S(T')}.
Actual path lengths will fluctuate around :math:`R(T)`. The average number of
photons emitted per unit path length is given by the inverse bremsstrahlung
mean free path:
.. math::
:label: inverse-bremsstrahlung-mfp
\lambda_{\text{br}}^{-1}(T,E_{\text{cut}})
= n\int_{E_{\text{cut}}}^T\frac{d\sigma_{\text{br}}}{dE}dE
= n\frac{Z^2}{\beta^2}\int_{\kappa_{\text{cut}}}^1\frac{1}{\kappa}
\chi(Z,T,\kappa)d\kappa.
The lower limit of the integral in :eq:`inverse-bremsstrahlung-mfp` is non-zero
because the bremsstrahlung differential cross section diverges for small photon
energies but is finite for photon energies above some cutoff energy
:math:`E_{\text{cut}}`. The mean free path
:math:`\lambda_{\text{br}}^{-1}(T,E_{\text{cut}})` is used to calculate the
photon number yield, defined as the average number of photons emitted with
energy greater than :math:`E_{\text{cut}}` as the charged particle slows down
from energy :math:`T` to :math:`E_{\text{cut}}`. The photon number yield is
given by
.. math::
:label: photon-number-yield
Y(T,E_{\text{cut}}) = \int^{R(T)}_{R(E_{\text{cut}})}
\lambda_{\text{br}}^{-1}(T',E_{\text{cut}})ds = \int_{E_{\text{cut}}}^T
\frac{\lambda_{\text{br}}^{-1}(T',E_{\text{cut}})}{S(T')}dT'.
:math:`Y(T,E_{\text{cut}})` can be used to construct the energy spectrum of
bremsstrahlung photons: the number of photons created with energy between
:math:`E_1` and :math:`E_2` by a charged particle with initial kinetic energy
:math:`T` as it comes to rest is given by :math:`Y(T,E_1) - Y(T,E_2)`.
To simulate the emission of bremsstrahlung photons, the total stopping power
and bremsstrahlung differential cross section for positrons and electrons must
be calculated for a given material using :eq:`material-bremsstrahlung-dcs` and
:eq:`material-radiative-stopping-power`. These quantities are used to build the
tabulated bremsstrahlung energy PDF and CDF for that material for each incident
energy :math:`T_k` on the energy grid. The following algorithm is then applied
to sample the photon energies:
1. For an incident charged particle with energy :math:`T`, sample the number of
emitted photons as
.. math::
N = \lfloor Y(T,E_{\text{cut}}) + \xi_1 \rfloor.
2. Rather than interpolate the PDF between indices :math:`k` and :math:`k+1`
for which :math:`T_k < T < T_{k+1}`, which is computationally expensive, use
the composition method and sample from the PDF at either :math:`k` or
:math:`k+1`. Using linear interpolation on a logarithmic scale, the PDF can
be expressed as
.. math::
p_{\text{br}}(T,E) = \pi_k p_{\text{br}}(T_k,E) + \pi_{k+1}
p_{\text{br}}(T_{k+1},E),
where the interpolation weights are
.. math::
\pi_k = \frac{\ln T_{k+1} - \ln T}{\ln T_{k+1} - \ln T_k},~~~
\pi_{k+1} = \frac{\ln T - \ln T_k}{\ln T_{k+1} - \ln T_k}.
Sample either the index :math:`i = k` or :math:`i = k+1` according to the
point probabilities :math:`\pi_{k}` and :math:`\pi_{k+1}`.
3. Determine the maximum value of the CDF :math:`P_{\text{br,max}}`.
3. Sample the photon energies using the inverse transform method with the
tabulated CDF :math:`P_{\text{br}}(T_i, E)` i.e.,
.. math::
E = E_j \left[ (1 + a_j) \frac{\xi_2 P_{\text{br,max}} -
P_{\text{br}}(T_i, E_j)} {E_j p_{\text{br}}(T_i, E_j)} + 1
\right]^{\frac{1}{1 + a_j}}
where the interpolation factor :math:`a_j` is given by
.. math::
a_j = \frac{\ln p_{\text{br}}(T_i,E_{j+1}) - \ln p_{\text{br}}(T_i,E_j)}
{\ln E_{j+1} - \ln E_j}
and :math:`P_{\text{br}}(T_i, E_j) \le \xi_2 P_{\text{br,max}} \le
P_{\text{br}}(T_i, E_{j+1})`.
We ignore the range of the electron or positron, i.e., the bremsstrahlung
photons are produced in the same location that the charged particle was
created. The direction of the photons is assumed to be the same as the
direction of the incident charged particle, which is a reasonable approximation
at higher energies when the bremsstrahlung radiation is emitted at small
angles.
Electron-Positron Annihilation
------------------------------
When a positron collides with an electron, both particles are annihilated and
generally two photons with equal energy are created. If the kinetic energy of
the positron is high enough, the two photons can have different energies, and
the higher-energy photon is emitted preferentially in the direction of flight
of the positron. It is also possible to produce a single photon if the
interaction occurs with a bound electron, and in some cases three (or, rarely,
even more) photons can be emitted. However, the annihilation cross section is
largest for low-energy positrons, and as the positron energy decreases, the
angular distribution of the emitted photons becomes isotropic.
In OpenMC, we assume the most likely case in which a low-energy positron (which
has already lost most of its energy to bremsstrahlung radiation) interacts with
an electron which is free and at rest. Two photons with energy equal to the
electron rest mass energy :math:`m_e c^2 = 0.511` MeV are emitted isotropically
in opposite directions.
.. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf
.. _Salvat: https://doi.org/10.1787/32da5043-en
.. _Sternheimer: https://doi.org/10.1103/PhysRevB.26.6067

View file

@ -14,6 +14,7 @@ Theory and Methodology
random_numbers
neutron_physics
photon_physics
charged_particles_physics
tallies
eigenvalue
depletion
@ -21,4 +22,4 @@ Theory and Methodology
parallelization
cmfd
variance_reduction
random_ray
random_ray

View file

@ -290,7 +290,10 @@ create and store fission sites for the following generation. First, the average
number of prompt and delayed neutrons must be determined to decide whether the
secondary neutrons will be prompt or delayed. This is important because delayed
neutrons have a markedly different spectrum from prompt neutrons, one that has a
lower average energy of emission. The total number of neutrons emitted
lower average energy of emission. Furthermore, in simulations where tracking
time of neutrons is important, we need to consider the emission time delay of
the secondary neutrons, which is dependent on the decay constant of the
delayed neutron precursor. The total number of neutrons emitted
:math:`\nu_t` is given as a function of incident energy in the ENDF format. Two
representations exist for :math:`\nu_t`. The first is a polynomial of order
:math:`N` with coefficients :math:`c_0,c_1,\dots,c_N`. If :math:`\nu_t` has this
@ -306,8 +309,8 @@ interpolation law. The number of prompt neutrons released per fission event
:math:`\nu_p` is also given as a function of incident energy and can be
specified in a polynomial or tabular format. The number of delayed neutrons
released per fission event :math:`\nu_d` can only be specified in a tabular
format. In practice, we only need to determine :math:`nu_t` and
:math:`nu_d`. Once these have been determined, we can calculated the delayed
format. In practice, we only need to determine :math:`\nu_t` and
:math:`\nu_d`. Once these have been determined, we can calculate the delayed
neutron fraction
.. math::
@ -335,8 +338,14 @@ neutrons. Otherwise, we produce :math:`\lfloor \nu \rfloor + 1` neutrons. Then,
for each fission site produced, we sample the outgoing angle and energy
according to the algorithms given in :ref:`sample-angle` and
:ref:`sample-energy` respectively. If the neutron is to be born delayed, then
there is an extra step of sampling a delayed neutron precursor group since they
each have an associated secondary energy distribution.
there is an extra step of sampling a delayed neutron precursor group to get the
associated secondary energy distribution and the decay constant
:math:`\lambda`, which is needed to sample the emission delay time :math:`t_d`:
.. math::
:label: sample-delay-time
t_d = -\frac{\ln \xi}{\lambda}.
The sampled outgoing angle and energy of fission neutrons along with the
position of the collision site are stored in an array called the fission

View file

@ -667,342 +667,6 @@ and Auger electrons:
5. Repeat from step 1 for vacancy left by the transition electron.
Electron-Positron Annihilation
------------------------------
When a positron collides with an electron, both particles are annihilated and
generally two photons with equal energy are created. If the kinetic energy of
the positron is high enough, the two photons can have different energies, and
the higher-energy photon is emitted preferentially in the direction of flight
of the positron. It is also possible to produce a single photon if the
interaction occurs with a bound electron, and in some cases three (or, rarely,
even more) photons can be emitted. However, the annihilation cross section is
largest for low-energy positrons, and as the positron energy decreases, the
angular distribution of the emitted photons becomes isotropic.
In OpenMC, we assume the most likely case in which a low-energy positron (which
has already lost most of its energy to bremsstrahlung radiation) interacts with
an electron which is free and at rest. Two photons with energy equal to the
electron rest mass energy :math:`m_e c^2 = 0.511` MeV are emitted isotropically
in opposite directions.
Bremsstrahlung
--------------
When a charged particle is decelerated in the field of an atom, some of its
kinetic energy is converted into electromagnetic radiation known as
bremsstrahlung, or 'braking radiation'. In each event, an electron or positron
with kinetic energy :math:`T` generates a photon with an energy :math:`E`
between :math:`0` and :math:`T`. Bremsstrahlung is described by a cross section
that is differential in photon energy, in the direction of the emitted photon,
and in the final direction of the charged particle. However, in Monte Carlo
simulations it is typical to integrate over the angular variables to obtain a
single differential cross section with respect to photon energy, which is often
expressed in the form
.. math::
:label: bremsstrahlung-dcs
\frac{d\sigma_{\text{br}}}{dE} = \frac{Z^2}{\beta^2} \frac{1}{E}
\chi(Z, T, \kappa),
where :math:`\kappa = E/T` is the reduced photon energy and :math:`\chi(Z, T,
\kappa)` is the scaled bremsstrahlung cross section, which is experimentally
measured.
Because electrons are attracted to atomic nuclei whereas positrons are
repulsed, the cross section for positrons is smaller, though it approaches that
of electrons in the high energy limit. To obtain the positron cross section, we
multiply :eq:`bremsstrahlung-dcs` by the :math:`\kappa`-independent factor used
in Salvat_,
.. math::
:label: positron-factor
\begin{aligned}
F_{\text{p}}(Z,T) =
& 1 - \text{exp}(-1.2359\times 10^{-1}t + 6.1274\times 10^{-2}t^2 - 3.1516\times 10^{-2}t^3 \\
& + 7.7446\times 10^{-3}t^4 - 1.0595\times 10^{-3}t^5 + 7.0568\times 10^{-5}t^6 \\
& - 1.8080\times 10^{-6}t^7),
\end{aligned}
where
.. math::
:label: positron-factor-t
t = \ln\left(1 + \frac{10^6}{Z^2}\frac{T}{\text{m}_\text{e}c^2} \right).
:math:`F_{\text{p}}(Z,T)` is the ratio of the radiative stopping powers for
positrons and electrons. Stopping power describes the average energy loss per
unit path length of a charged particle as it passes through matter:
.. math::
:label: stopping-power
-\frac{dT}{ds} = n \int E \frac{d\sigma}{dE} dE \equiv S(T),
where :math:`n` is the number density of the material and :math:`d\sigma/dE` is
the cross section differential in energy loss. The total stopping power
:math:`S(T)` can be separated into two components: the radiative stopping
power :math:`S_{\text{rad}}(T)`, which refers to energy loss due to
bremsstrahlung, and the collision stopping power :math:`S_{\text{col}}(T)`,
which refers to the energy loss due to inelastic collisions with bound
electrons in the material that result in ionization and excitation. The
radiative stopping power for electrons is given by
.. math::
:label: radiative-stopping-power
S_{\text{rad}}(T) = n \frac{Z^2}{\beta^2} T \int_0^1 \chi(Z,T,\kappa)
d\kappa.
To obtain the radiative stopping power for positrons,
:eq:`radiative-stopping-power` is multiplied by :eq:`positron-factor`.
While the models for photon interactions with matter described above can safely
assume interactions occur with free atoms, sampling the target atom based on
the macroscopic cross sections, molecular effects cannot necessarily be
disregarded for charged particle treatment. For compounds and mixtures, the
bremsstrahlung cross section is calculated using Bragg's additivity rule as
.. math::
:label: material-bremsstrahlung-dcs
\frac{d\sigma_{\text{br}}}{dE} = \frac{1}{\beta^2 E} \sum_i \gamma_i Z^2_i
\chi(Z_i, T, \kappa),
where the sum is over the constituent elements and :math:`\gamma_i` is the
atomic fraction of the :math:`i`-th element. Similarly, the radiative stopping
power is calculated using Bragg's additivity rule as
.. math::
:label: material-radiative-stopping-power
S_{\text{rad}}(T) = \sum_i w_i S_{\text{rad},i}(T),
where :math:`w_i` is the mass fraction of the :math:`i`-th element and
:math:`S_{\text{rad},i}(T)` is found for element :math:`i` using
:eq:`radiative-stopping-power`. The collision stopping power, however, is a
function of certain quantities such as the mean excitation energy :math:`I` and
the density effect correction :math:`\delta_F` that depend on molecular
properties. These quantities cannot simply be summed over constituent elements
in a compound, but should instead be calculated for the material. The Bethe
formula can be used to find the collision stopping power of the material:
.. math::
:label: material-collision-stopping-power
S_{\text{col}}(T) = \frac{2 \pi r_e^2 m_e c^2}{\beta^2} N_A \frac{Z}{A_M}
[\ln(T^2/I^2) + \ln(1 + \tau/2) + F(\tau) - \delta_F(T)],
where :math:`N_A` is Avogadro's number, :math:`A_M` is the molar mass,
:math:`\tau = T/m_e`, and :math:`F(\tau)` depends on the particle type. For
electrons,
.. math::
:label: F-electron
F_{-}(\tau) = (1 - \beta^2)[1 + \tau^2/8 - (2\tau + 1) \ln2],
while for positrons
.. math::
:label: F-positron
F_{+}(\tau) = 2\ln2 - (\beta^2/12)[23 + 14/(\tau + 2) + 10/(\tau + 2)^2 +
4/(\tau + 2)^3].
The density effect correction :math:`\delta_F` takes into account the reduction
of the collision stopping power due to the polarization of the material the
charged particle is passing through by the electric field of the particle.
It can be evaluated using the method described by Sternheimer_, where the
equation for :math:`\delta_F` is
.. math::
:label: density-effect-correction
\delta_F(\beta) = \sum_{i=1}^n f_i \ln[(l_i^2 + l^2)/l_i^2] -
l^2(1-\beta^2).
Here, :math:`f_i` is the oscillator strength of the :math:`i`-th transition,
given by :math:`f_i = n_i/Z`, where :math:`n_i` is the number of electrons in
the :math:`i`-th subshell. The frequency :math:`l` is the solution of the
equation
.. math::
:label: density-effect-l
\frac{1}{\beta^2} - 1 = \sum_{i=1}^{n} \frac{f_i}{\bar{\nu}_i^2 + l^2},
where :math:`\bar{v}_i` is defined as
.. math::
:label: density-effect-nubar
\bar{\nu}_i = h\nu_i \rho / h\nu_p.
The plasma energy :math:`h\nu_p` of the medium is given by
.. math::
:label: plasma-frequency
h\nu_p = \sqrt{\frac{(hc)^2 r_e \rho_m N_A Z}{\pi A}},
where :math:`A` is the atomic weight and :math:`\rho_m` is the density of the
material. In :eq:`density-effect-nubar`, :math:`h\nu_i` is the oscillator
energy, and :math:`\rho` is an adjustment factor introduced to give agreement
between the experimental values of the oscillator energies and the mean
excitation energy. The :math:`l_i` in :eq:`density-effect-correction` are
defined as
.. math::
:label: density-effect-li
\begin{aligned}
l_i &= (\bar{\nu}_i^2 + 2/3f_i)^{1/2} ~~~~&\text{for}~~ \bar{\nu}_i > 0 \\
l_n &= f_n^{1/2} ~~~~&\text{for}~~ \bar{\nu}_n = 0,
\end{aligned}
where the second case applies to conduction electrons. For a conductor,
:math:`f_n` is given by :math:`n_c/Z`, where :math:`n_c` is the effective
number of conduction electrons, and :math:`v_n = 0`. The adjustment factor
:math:`\rho` is determined using the equation for the mean excitation energy:
.. math::
:label: mean-excitation-energy
\ln I = \sum_{i=1}^{n-1} f_i \ln[(h\nu_i\rho)^2 + 2/3f_i(h\nu_p)^2]^{1/2} +
f_n \ln (h\nu_pf_n^{1/2}).
.. _ttb:
Thick-Target Bremsstrahlung Approximation
+++++++++++++++++++++++++++++++++++++++++
Since charged particles lose their energy on a much shorter distance scale than
neutral particles, not much error should be introduced by neglecting to
transport electrons. However, the bremsstrahlung emitted from high energy
electrons and positrons can travel far from the interaction site. Thus, even
without a full electron transport mode it is necessary to model bremsstrahlung.
We use a thick-target bremsstrahlung (TTB) approximation based on the models in
Salvat_ and Kaltiaisenaho_ for generating bremsstrahlung photons, which assumes
the charged particle loses all its energy in a single homogeneous material
region.
To model bremsstrahlung using the TTB approximation, we need to know the number
of photons emitted by the charged particle and the energy distribution of the
photons. These quantities can be calculated using the continuous slowing down
approximation (CSDA). The CSDA assumes charged particles lose energy
continuously along their trajectory with a rate of energy loss equal to the
total stopping power, ignoring fluctuations in the energy loss. The
approximation is useful for expressing average quantities that describe how
charged particles slow down in matter. For example, the CSDA range approximates
the average path length a charged particle travels as it slows to rest:
.. math::
:label: csda-range
R(T) = \int^T_0 \frac{dT'}{S(T')}.
Actual path lengths will fluctuate around :math:`R(T)`. The average number of
photons emitted per unit path length is given by the inverse bremsstrahlung
mean free path:
.. math::
:label: inverse-bremsstrahlung-mfp
\lambda_{\text{br}}^{-1}(T,E_{\text{cut}})
= n\int_{E_{\text{cut}}}^T\frac{d\sigma_{\text{br}}}{dE}dE
= n\frac{Z^2}{\beta^2}\int_{\kappa_{\text{cut}}}^1\frac{1}{\kappa}
\chi(Z,T,\kappa)d\kappa.
The lower limit of the integral in :eq:`inverse-bremsstrahlung-mfp` is non-zero
because the bremsstrahlung differential cross section diverges for small photon
energies but is finite for photon energies above some cutoff energy
:math:`E_{\text{cut}}`. The mean free path
:math:`\lambda_{\text{br}}^{-1}(T,E_{\text{cut}})` is used to calculate the
photon number yield, defined as the average number of photons emitted with
energy greater than :math:`E_{\text{cut}}` as the charged particle slows down
from energy :math:`T` to :math:`E_{\text{cut}}`. The photon number yield is
given by
.. math::
:label: photon-number-yield
Y(T,E_{\text{cut}}) = \int^{R(T)}_{R(E_{\text{cut}})}
\lambda_{\text{br}}^{-1}(T',E_{\text{cut}})ds = \int_{E_{\text{cut}}}^T
\frac{\lambda_{\text{br}}^{-1}(T',E_{\text{cut}})}{S(T')}dT'.
:math:`Y(T,E_{\text{cut}})` can be used to construct the energy spectrum of
bremsstrahlung photons: the number of photons created with energy between
:math:`E_1` and :math:`E_2` by a charged particle with initial kinetic energy
:math:`T` as it comes to rest is given by :math:`Y(T,E_1) - Y(T,E_2)`.
To simulate the emission of bremsstrahlung photons, the total stopping power
and bremsstrahlung differential cross section for positrons and electrons must
be calculated for a given material using :eq:`material-bremsstrahlung-dcs` and
:eq:`material-radiative-stopping-power`. These quantities are used to build the
tabulated bremsstrahlung energy PDF and CDF for that material for each incident
energy :math:`T_k` on the energy grid. The following algorithm is then applied
to sample the photon energies:
1. For an incident charged particle with energy :math:`T`, sample the number of
emitted photons as
.. math::
N = \lfloor Y(T,E_{\text{cut}}) + \xi_1 \rfloor.
2. Rather than interpolate the PDF between indices :math:`k` and :math:`k+1`
for which :math:`T_k < T < T_{k+1}`, which is computationally expensive, use
the composition method and sample from the PDF at either :math:`k` or
:math:`k+1`. Using linear interpolation on a logarithmic scale, the PDF can
be expressed as
.. math::
p_{\text{br}}(T,E) = \pi_k p_{\text{br}}(T_k,E) + \pi_{k+1}
p_{\text{br}}(T_{k+1},E),
where the interpolation weights are
.. math::
\pi_k = \frac{\ln T_{k+1} - \ln T}{\ln T_{k+1} - \ln T_k},~~~
\pi_{k+1} = \frac{\ln T - \ln T_k}{\ln T_{k+1} - \ln T_k}.
Sample either the index :math:`i = k` or :math:`i = k+1` according to the
point probabilities :math:`\pi_{k}` and :math:`\pi_{k+1}`.
3. Determine the maximum value of the CDF :math:`P_{\text{br,max}}`.
3. Sample the photon energies using the inverse transform method with the
tabulated CDF :math:`P_{\text{br}}(T_i, E)` i.e.,
.. math::
E = E_j \left[ (1 + a_j) \frac{\xi_2 P_{\text{br,max}} -
P_{\text{br}}(T_i, E_j)} {E_j p_{\text{br}}(T_i, E_j)} + 1
\right]^{\frac{1}{1 + a_j}}
where the interpolation factor :math:`a_j` is given by
.. math::
a_j = \frac{\ln p_{\text{br}}(T_i,E_{j+1}) - \ln p_{\text{br}}(T_i,E_j)}
{\ln E_{j+1} - \ln E_j}
and :math:`P_{\text{br}}(T_i, E_j) \le \xi_2 P_{\text{br,max}} \le
P_{\text{br}}(T_i, E_{j+1})`.
We ignore the range of the electron or positron, i.e., the bremsstrahlung
photons are produced in the same location that the charged particle was
created. The direction of the photons is assumed to be the same as the
direction of the incident charged particle, which is a reasonable approximation
at higher energies when the bremsstrahlung radiation is emitted at small
angles.
.. _photon_production:
@ -1070,5 +734,3 @@ emitted photon.
.. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf
.. _Salvat: https://doi.org/10.1787/32da5043-en
.. _Sternheimer: https://doi.org/10.1103/PhysRevB.26.6067

View file

@ -387,6 +387,101 @@ of this is that the longer you run a simulation, the better you know your
results. Therefore, by running a simulation long enough, it is possible to
reduce the stochastic uncertainty to arbitrarily low levels.
Skewness
++++++++
The `skewness`_ of a population quantifies the asymmetry of the probability
distribution around its mean. Positive and negative skewness indicate a
longer/heavier right and left tail respectively. Let :math:`x_1,\ldots,x_n` be
the per-realization values for a bin, with sample mean :math:`\bar{x}` and
sample central moments:
.. math::
m_k \;=\; \frac{1}{n}\sum_{i=1}^{n}\bigl(x_i-\bar{x}\bigr)^k.
OpenMC reports the *adjusted Fisher-Pearson skewness* (defined for :math:`n \ge
3`), which is commonly used in many statistical packages:
.. math::
G_1 \;=\; \frac{\sqrt{n \cdot (n-1)}}{\,n-2\,}\cdot\frac{m_3}{m_2^{3/2}}.
where :math:`m_2` and :math:`m_3` correspond to the biased sample second and
third central moment respectively.
Kurtosis
++++++++
The `kurtosis`_ of a population quantifies tail weight (also called tailedness)
of the probability distribution relative to a normal distribution. Positive
excess kurtosis indicates *heavier tails* whereas negative excess kurtosis
indicates *lighter tails*. Kurtosis is especially useful for identifying bins
where occasional extreme scores dominate uncertainty. OpenMC reports the
*adjusted excess kurtosis* (defined for :math:`n \ge 4`):
.. math::
G_2 \;=\; \frac{(n-1)}{(n-2)(n-3)}
\left[(n+1)\,\frac{m_4}{m_2^{2}} \;-\; 3(n-1)\right].
where :math:`m_2` and :math:`m_4` correspond to the biased sample second and
fourth central moment respectively. For a perfectly normal distribution, the
excess kurtosis is :math:`0`.
Variance of Variance
++++++++++++++++++++
The variance of the variance (also known as the coefficient of variation
squared) measures *stability of the sample variance* :math:`s^2` and, by
extension, the reliability of reported relative errors. High VOV means that
error bars themselves are noisy—often due to heavy tails, skewness, or too few
realizations.
.. math::
VOV = \frac{s^2(s_{\bar{X}}^2)}{s_{\bar{X}}^4 } = \frac{m_4}{m_2^2} - \frac{1}{n}
where :math:`s_{\bar{X}}^2` is the estimated variance of the mean and
:math:`s^2(s_{\bar{X}}^2)` is the estimated variance in :math:`s_{\bar{X}}^2`.
The MCNP manual suggests a hard threshold such that :math:`VOV < 0.1` to improve
the probability of forming a reliable confidence interval. However, OpenMC does
not enforce an universal cut-off because the suitability of any single threshold
depends strongly on problem specifics (estimator choice, variance-reduction
settings, tally binning, or even effective sample size).
Normality Tests (D'Agostino-Pearson)
++++++++++++++++++++++++++++++++++++
These normality test verify the hypothesis that fluctuations are *approximately
normal*, a working assumption behind many Monte Carlo diagnostics and
`confidence-interval heuristics`_. Tests are provided for: (i) skewness-only,
(ii) kurtosis-only, and (iii) the *omnibus* combination. OpenMC uses the
finite-sample-adjusted skewness :math:`G_1` and excess kurtosis :math:`G_2`
above to construct standardized normal scores :math:`Z_1` (from :math:`G_1`) and
:math:`Z_2` (from :math:`G_2`) via the D'Agostino-Pearson transformations. The
omnibus statistic is
.. math::
K^2 \;=\; Z_1^{\,2} \;+\; Z_2^{\,2}
\;\sim\; \chi^2_{(2)} \quad \text{under } H_0:\ \text{normality}.
OpenMC reports :math:`Z_1`, :math:`Z_2`, :math:`K^2`, and their p-values when
prerequisites are met (skewness for :math:`n\ge 3`, kurtosis and omnibus for
:math:`n\ge 4`). Given a user-chosen significance level :math:`\alpha` (default
is :math:`0.05`), reject :math:`H_0` if :math:`\text{p-value}<\alpha`; otherwise
fail to reject. OpenMC leaves the interpretation to the user, who should
consider VOV together with skewness, kurtosis, and normality tests results when
judging whether reported confidence intervals are credible for their application
[#norm-tests]_.
.. [#norm-tests]
Higher-moments accumulation must be enabled with ``higher_moments = True``
for running these diagnostics including the skewness, kurtosis, and normality
tests.
Figure of Merit
+++++++++++++++
@ -405,14 +500,16 @@ defined as
.. math::
:label: relative_error
r = \frac{s_\bar{X}}{\bar{x}}.
r = \frac{s_{\bar{X}}}{\bar{x}}.
Based on this definition, one can see that a higher FOM is desirable. The FOM is
useful as a comparative tool. For example, if a variance reduction technique is
being applied to a simulation, the FOM with variance reduction can be compared
to the FOM without variance reduction to ascertain whether the reduction in
variance outweighs the potential increase in execution time (e.g., due to
particle splitting).
particle splitting). It is important to note that MCNP reports the FOM using CPU
time (wall-clock time multiplied by the number of threads/cores), whereas OpenMC
reports the FOM using only the wall-clock time :math:`t`.
Confidence Intervals
++++++++++++++++++++
@ -521,6 +618,8 @@ improve the estimate of the percentile.
.. rubric:: References
.. _confidence-interval heuristics: https://doi.org/10.1080/00031305.1990.10475751
.. _following approximation: https://doi.org/10.1080/03610918708812641
.. _Bessel's correction: https://en.wikipedia.org/wiki/Bessel's_correction
@ -541,6 +640,10 @@ improve the estimate of the percentile.
.. _converges in distribution: https://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_distribution
.. _skewness: https://en.wikipedia.org/wiki/Skewness
.. _kurtosis: https://en.wikipedia.org/wiki/Kurtosis
.. _confidence intervals: https://en.wikipedia.org/wiki/Confidence_interval
.. _Student's t-distribution: https://en.wikipedia.org/wiki/Student%27s_t-distribution

View file

@ -176,7 +176,8 @@ Geometry Plotting
:nosignatures:
:template: myclass.rst
openmc.Plot
openmc.SlicePlot
openmc.VoxelPlot
openmc.WireframeRayTracePlot
openmc.SolidRayTracePlot
openmc.Plots
@ -216,6 +217,9 @@ Post-processing
:nosignatures:
:template: myfunction.rst
openmc.read_collision_track_file
openmc.read_collision_track_hdf5
openmc.read_collision_track_mcpl
openmc.voxel_to_vtk
The following classes and functions are used for functional expansion reconstruction.

View file

@ -287,6 +287,16 @@ the following abstract base classes:
abc.SIIntegrator
abc.DepSystemSolver
R2S Automation
--------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
R2SManager
D1S Functions
-------------

View file

@ -0,0 +1,226 @@
====================
What's New in 0.15.3
====================
.. currentmodule:: openmc
-------
Summary
-------
This release of OpenMC includes many bug fixes, performance improvements, and
several notable new features. The major highlights of this release include a new
:class:`~openmc.deplete.R2SManager` class that automates the workflow for
rigorous 2-step (R2S) shutdown dose rate calculations, the ability to collect
higher moments for tally results that can be used to test normality, a new
uncertainty-aware criticality search method, a new collision tracking feature
that enables detailed tracking of particle interactions, support for distributed
cell densities, and several new tally filters. The random ray solver also
continues to receive significant updates, including automatic setup
capabilities, improved geometry handling, and better weight window support.
Depletion capabilities have been expanded with thermochemical redox control,
external transfer rates, and improved performance.
------------------------------------
Compatibility Notes and Deprecations
------------------------------------
MCPL has been changed from a build-time dependency to a runtime optional
dependency, which means OpenMC will attempt to load the MCPL library at
runtime when needed rather than requiring it at build time.
The ``openmc.mgxs.Library.add_to_tallies_file`` method has been renamed to
:meth:`openmc.mgxs.Library.add_to_tallies`.
------------
New Features
------------
- A new collision tracking feature enables detailed tracking of particle
interactions (`#3417 <https://github.com/openmc-dev/openmc/pull/3417>`_)
- Added :meth:`~openmc.model.Model.keff_search` method for automated criticality
searches (`#3569 <https://github.com/openmc-dev/openmc/pull/3569>`_)
- Introduced automated workflow for mesh- or cell-based R2S calculations
(`#3508 <https://github.com/openmc-dev/openmc/pull/3508>`_)
- Ability to source electron/positrons directly for charged particle
simulations (`#3404 <https://github.com/openmc-dev/openmc/pull/3404>`_)
- Multi-group capability for kinetics parameter calculations with Iterated
Fission Probability (`#3425
<https://github.com/openmc-dev/openmc/pull/3425>`_)
- Introduced a new :class:`openmc.MeshMaterialFilter` class (`#3406
<https://github.com/openmc-dev/openmc/pull/3406>`_)
- Added support for distributed cell densities (`#3546
<https://github.com/openmc-dev/openmc/pull/3546>`_)
- Implemented a :class:`openmc.WeightWindowsList` class that enables export to
HDF5 (`#3456 <https://github.com/openmc-dev/openmc/pull/3456>`_)
- Added :meth:`openmc.Material.mean_free_path` method (`#3469
<https://github.com/openmc-dev/openmc/pull/3469>`_)
- Introduced :func:`openmc.lib.TemporarySession` context manager (`#3475
<https://github.com/openmc-dev/openmc/pull/3475>`_)
- Added material depletion function for tracking individual material depletion
(`#3420 <https://github.com/openmc-dev/openmc/pull/3420>`_)
- Added methods on :class:`~openmc.Material` class for waste disposal rating /
classification (`#3366 <https://github.com/openmc-dev/openmc/pull/3366>`_,
`#3376 <https://github.com/openmc-dev/openmc/pull/3376>`_)
- Support for thermochemical redox control transfer rates in depletion
(`#2783 <https://github.com/openmc-dev/openmc/pull/2783>`_)
- Support for external transfer rates source term in depletion (`#3088
<https://github.com/openmc-dev/openmc/pull/3088>`_)
- Added combing capability for fission site sampling and delayed neutron
emission time (`#2992 <https://github.com/openmc-dev/openmc/pull/2992>`_)
- Ability to specify reference direction for azimuthal angle in
:class:`~openmc.stats.PolarAzimuthal` distribution (`#3582
<https://github.com/openmc-dev/openmc/pull/3582>`_)
- Allow spatial constraints on element sources within
:class:`~openmc.MeshSource` (`#3431
<https://github.com/openmc-dev/openmc/pull/3431>`_)
- Added VTK HDF (.vtkhdf) format support for writing VTK data (`#3252
<https://github.com/openmc-dev/openmc/pull/3252>`_)
- Implemented filter weight capability (`#3345
<https://github.com/openmc-dev/openmc/pull/3345>`_)
- Optionally collect higher moments for tallies (`#3363
<https://github.com/openmc-dev/openmc/pull/3363>`_)
- Several random ray solver enhancements:
- Random Ray AutoMagic Setup for automatic configuration (`#3351 <https://github.com/openmc-dev/openmc/pull/3351>`_)
- Point source locator for random ray mode (`#3360 <https://github.com/openmc-dev/openmc/pull/3360>`_)
- Support for DAGMC geometries (`#3374 <https://github.com/openmc-dev/openmc/pull/3374>`_)
- Optimized mapping of source regions to tallies (`#3465 <https://github.com/openmc-dev/openmc/pull/3465>`_)
- Base source region refactor (`#3576 <https://github.com/openmc-dev/openmc/pull/3576>`_)
---------------------------
Bug Fixes and Small Changes
---------------------------
- Add two MPI barriers in R2S workflow (`#3646 <https://github.com/openmc-dev/openmc/pull/3646>`_)
- Fix a few warnings, rename add_to_tallies_file (`#3639 <https://github.com/openmc-dev/openmc/pull/3639>`_)
- Fix typo in DAGMC lost particle test (`#3634 <https://github.com/openmc-dev/openmc/pull/3634>`_)
- Avoid multiprocessing Pool when running depletion tests with MPI (`#3633 <https://github.com/openmc-dev/openmc/pull/3633>`_)
- Support MPI parallelism in R2SManager (`#3632 <https://github.com/openmc-dev/openmc/pull/3632>`_)
- Update documentation for particle tracks (`#3627 <https://github.com/openmc-dev/openmc/pull/3627>`_)
- Adding variance of variance and normality tests for tally statistics (`#3454 <https://github.com/openmc-dev/openmc/pull/3454>`_)
- Avoid divide-by-zero in ``from_multigroup_flux`` when flux is zero (`#3624 <https://github.com/openmc-dev/openmc/pull/3624>`_)
- Write particle states as separate lines in track VTK files (`#3628 <https://github.com/openmc-dev/openmc/pull/3628>`_)
- Reset DAGMC history when reviving from source (`#3601 <https://github.com/openmc-dev/openmc/pull/3601>`_)
- Add energy group structure: SCALE-999 (`#3564 <https://github.com/openmc-dev/openmc/pull/3564>`_)
- Fix bug in normalization of tally results with no_reduce (`#3619 <https://github.com/openmc-dev/openmc/pull/3619>`_)
- Enable nuclide filters with get_decay_photon_energy (`#3614 <https://github.com/openmc-dev/openmc/pull/3614>`_)
- Update ``check_type`` calls to accept both ``str`` and ``os.PathLike`` objects (`#3618 <https://github.com/openmc-dev/openmc/pull/3618>`_)
- Speed up ``apply_time_correction`` by reducing file I/O and deepcopies (`#3617 <https://github.com/openmc-dev/openmc/pull/3617>`_)
- FW-CADIS Disregard Max Realizations Setting (`#3616 <https://github.com/openmc-dev/openmc/pull/3616>`_)
- Random Ray Geometry Debug Mode Fix (`#3615 <https://github.com/openmc-dev/openmc/pull/3615>`_)
- Don't write reaction rates in depletion results by default (`#3609 <https://github.com/openmc-dev/openmc/pull/3609>`_)
- Allow Path objects in MGXSLibrary.export_to_hdf5 (`#3608 <https://github.com/openmc-dev/openmc/pull/3608>`_)
- Clip mixture distributions based on mean times integral (`#3603 <https://github.com/openmc-dev/openmc/pull/3603>`_)
- Allow V0 in atomic_mass function (for ENDF/B-VII.0 data) (`#3607 <https://github.com/openmc-dev/openmc/pull/3607>`_)
- Re-run flaky tests when needed (`#3604 <https://github.com/openmc-dev/openmc/pull/3604>`_)
- Ability to load mesh objects from weight_windows.h5 file (`#3598 <https://github.com/openmc-dev/openmc/pull/3598>`_)
- Switch to using coveralls github action for reporting (`#3594 <https://github.com/openmc-dev/openmc/pull/3594>`_)
- Add user setting for free gas threshold (`#3593 <https://github.com/openmc-dev/openmc/pull/3593>`_)
- Speed up time correction factors (`#3592 <https://github.com/openmc-dev/openmc/pull/3592>`_)
- Fix caching issue when using NCrystal materials (`#3538 <https://github.com/openmc-dev/openmc/pull/3538>`_)
- Fix random ray source region mesh export when using model.export_to_xml() (`#3579 <https://github.com/openmc-dev/openmc/pull/3579>`_)
- Ensure weight_windows_file information is read from XML (`#3587 <https://github.com/openmc-dev/openmc/pull/3587>`_)
- Add missing documentation on <source> in depletion chain file format (`#3590 <https://github.com/openmc-dev/openmc/pull/3590>`_)
- Adding tally filter type option to statepoint get_tally (`#3584 <https://github.com/openmc-dev/openmc/pull/3584>`_)
- Optional separation of mesh-material-volume calc from get_homogenized_materials (`#3581 <https://github.com/openmc-dev/openmc/pull/3581>`_)
- Fix IFP implementation (`#3580 <https://github.com/openmc-dev/openmc/pull/3580>`_)
- Remove several TODOs related to C++17 support (`#3574 <https://github.com/openmc-dev/openmc/pull/3574>`_)
- Fix performance regression in libMesh unstructured mesh tallies (`#3577 <https://github.com/openmc-dev/openmc/pull/3577>`_)
- Update find_package calls in OpenMCConfig.cmake (`#3572 <https://github.com/openmc-dev/openmc/pull/3572>`_)
- Ensure ``n_dimension_`` attribute is set for unstructured meshes (`#3575 <https://github.com/openmc-dev/openmc/pull/3575>`_)
- Allow newer Sphinx version and fix docbuild warnings (`#3571 <https://github.com/openmc-dev/openmc/pull/3571>`_)
- Fixed a bug when combining TimeFilter, MeshFilter, and tracklength estimator (`#3525 <https://github.com/openmc-dev/openmc/pull/3525>`_)
- PowerLaw raises an error if sampling interval contains negative values (`#3542 <https://github.com/openmc-dev/openmc/pull/3542>`_)
- depletion: fix performance of chain matrix construction (`#3567 <https://github.com/openmc-dev/openmc/pull/3567>`_)
- Do not apply boundary conditions when initialized in volume calculation mode (`#3562 <https://github.com/openmc-dev/openmc/pull/3562>`_)
- Bump up tolerance for flaky activation test (`#3560 <https://github.com/openmc-dev/openmc/pull/3560>`_)
- Fixed a bug in plotting cross sections with S(a,b) data (`#3558 <https://github.com/openmc-dev/openmc/pull/3558>`_)
- Change test order to run unit tests first (`#3533 <https://github.com/openmc-dev/openmc/pull/3533>`_)
- adding ecco 33 (`#3556 <https://github.com/openmc-dev/openmc/pull/3556>`_)
- Refactor endf_data to be a fixture (`#3539 <https://github.com/openmc-dev/openmc/pull/3539>`_)
- Revert "fix broken CI" (`#3554 <https://github.com/openmc-dev/openmc/pull/3554>`_)
- fix broken CI (`#3551 <https://github.com/openmc-dev/openmc/pull/3551>`_)
- Leverage particle.move_distance in event advance (`#3544 <https://github.com/openmc-dev/openmc/pull/3544>`_)
- fix tests that accidentaly got broken (`#3543 <https://github.com/openmc-dev/openmc/pull/3543>`_)
- not printing nuclides with 0 percent to terminal (option 2 ) (`#3448 <https://github.com/openmc-dev/openmc/pull/3448>`_)
- Fix a bug in time cutoff behavior (`#3526 <https://github.com/openmc-dev/openmc/pull/3526>`_)
- Avoid duplicate materials written to XML (`#3536 <https://github.com/openmc-dev/openmc/pull/3536>`_)
- Use cached property for openmc.data.Decay.sources (`#3535 <https://github.com/openmc-dev/openmc/pull/3535>`_)
- more helpful error message for dose_coefficients (`#3534 <https://github.com/openmc-dev/openmc/pull/3534>`_)
- Adding 616 group structure (`#3531 <https://github.com/openmc-dev/openmc/pull/3531>`_)
- Remove unused special accessors for tallies (`#3527 <https://github.com/openmc-dev/openmc/pull/3527>`_)
- Consistent XML parsing using functions from _xml module (`#3517 <https://github.com/openmc-dev/openmc/pull/3517>`_)
- Add stat:sum field to MCPL files for proper weight normalization (`#3522 <https://github.com/openmc-dev/openmc/pull/3522>`_)
- Remove reorder_attributes from openmc._xml (`#3519 <https://github.com/openmc-dev/openmc/pull/3519>`_)
- fixed a bug in MeshMaterialFilter.from_volumes (`#3520 <https://github.com/openmc-dev/openmc/pull/3520>`_)
- Fixed a bug in distribcell offsets logic (`#3424 <https://github.com/openmc-dev/openmc/pull/3424>`_)
- Add test for FW-CADIS based WW generation on a DAGMC model (`#3504 <https://github.com/openmc-dev/openmc/pull/3504>`_)
- Fix for Weight Window Scaling Bug (`#3511 <https://github.com/openmc-dev/openmc/pull/3511>`_)
- Fix: ``materials``, ``plots``, and ``tallies`` cannot be passed as lists (`#3513 <https://github.com/openmc-dev/openmc/pull/3513>`_)
- Allow already-initialized openmc.lib in TemporarySession (`#3505 <https://github.com/openmc-dev/openmc/pull/3505>`_)
- Update DAGMC and libMesh precompiler definitions (`#3510 <https://github.com/openmc-dev/openmc/pull/3510>`_)
- Avoid adding ParentNuclideFilter twice when calling prepare_tallies (`#3506 <https://github.com/openmc-dev/openmc/pull/3506>`_)
- Enabling MCPL source files to be read when using surf_source_read (`#3472 <https://github.com/openmc-dev/openmc/pull/3472>`_)
- Boundary info accessors (`#3496 <https://github.com/openmc-dev/openmc/pull/3496>`_)
- automatically finding appropriate dimension when making regular mesh from domain (`#3468 <https://github.com/openmc-dev/openmc/pull/3468>`_)
- Add accessor methods for LocalCoord (`#3494 <https://github.com/openmc-dev/openmc/pull/3494>`_)
- Make MCPL a Runtime Optional Dependency (`#3429 <https://github.com/openmc-dev/openmc/pull/3429>`_)
- Use auto-chunking for StepResult HDF5 writing (`#3498 <https://github.com/openmc-dev/openmc/pull/3498>`_)
- Provide a way to get ID maps from plot parameters on the Model class (`#3481 <https://github.com/openmc-dev/openmc/pull/3481>`_)
- Update OSX install instructions to point to x64 platform (`#3501 <https://github.com/openmc-dev/openmc/pull/3501>`_)
- Update conda install instructions for macOS Apple silicon (`#3488 <https://github.com/openmc-dev/openmc/pull/3488>`_)
- Only show warning if in restart mode (`#3478 <https://github.com/openmc-dev/openmc/pull/3478>`_)
- Add flag to CMakeLists to use submodules instead of searching (`#3480 <https://github.com/openmc-dev/openmc/pull/3480>`_)
- Added citation metadata file (`#3409 <https://github.com/openmc-dev/openmc/pull/3409>`_)
- fix zam parsing (`#3484 <https://github.com/openmc-dev/openmc/pull/3484>`_)
- Support flux collapse method in ``get_microxs_and_flux`` (`#3466 <https://github.com/openmc-dev/openmc/pull/3466>`_)
- Stabilize Adjoint Source (`#3476 <https://github.com/openmc-dev/openmc/pull/3476>`_)
- Refactor and Harden Configuration Management (`#3461 <https://github.com/openmc-dev/openmc/pull/3461>`_)
- Updated Docs to Not Give Specific Python Version Requirement (`#3473 <https://github.com/openmc-dev/openmc/pull/3473>`_)
- Parallelization of Weight Window Update (`#3467 <https://github.com/openmc-dev/openmc/pull/3467>`_)
- Limit Random Ray Weight Window Generation to Final Batch (`#3464 <https://github.com/openmc-dev/openmc/pull/3464>`_)
- Fix Dockerfile DAGMC build (`#3463 <https://github.com/openmc-dev/openmc/pull/3463>`_)
- Fix Weight Window Infinite Loop Bug (`#3457 <https://github.com/openmc-dev/openmc/pull/3457>`_)
- Weight Window Birth Scaling (`#3459 <https://github.com/openmc-dev/openmc/pull/3459>`_)
- Adding checks to geometry.plot to avoid material name overlaps (`#3458 <https://github.com/openmc-dev/openmc/pull/3458>`_)
- Fixing crash when calling Geometry.plot when DAGMCUniverse in geometry (`#3455 <https://github.com/openmc-dev/openmc/pull/3455>`_)
- fixing expansion of elemental Ta bug (`#3443 <https://github.com/openmc-dev/openmc/pull/3443>`_)
- Prevent Adjoint Sources from Trending towards Infinity (`#3449 <https://github.com/openmc-dev/openmc/pull/3449>`_)
- adding plot function to DAGMCUnvierse (`#3451 <https://github.com/openmc-dev/openmc/pull/3451>`_)
- Allow specifying number of equiprobable angles for thermal scattering data generation (`#3346 <https://github.com/openmc-dev/openmc/pull/3346>`_)
- Change Dockerfile from debian:bookworm-slim to ubuntu:24.04 (`#3442 <https://github.com/openmc-dev/openmc/pull/3442>`_)
- Fix Resetting of Auto IDs When Generating MGXS (`#3437 <https://github.com/openmc-dev/openmc/pull/3437>`_)
- Allowing chain_file to be chain object to save reloading time (`#3436 <https://github.com/openmc-dev/openmc/pull/3436>`_)
- update units for flux (`#3441 <https://github.com/openmc-dev/openmc/pull/3441>`_)
- Fix raytrace infinite loop (`#3423 <https://github.com/openmc-dev/openmc/pull/3423>`_)
- Apply Max Number of Events Check to Random Rays (`#3438 <https://github.com/openmc-dev/openmc/pull/3438>`_)
- Add user setting for source rejection fraction (`#3433 <https://github.com/openmc-dev/openmc/pull/3433>`_)
- Adding fix and tests for spherical mesh as spatial distribution (`#3428 <https://github.com/openmc-dev/openmc/pull/3428>`_)
- Random Ray Missed Cell Policy Change for Adjoint Mode (`#3434 <https://github.com/openmc-dev/openmc/pull/3434>`_)
- Random Ray External Source Plotting Fix (`#3430 <https://github.com/openmc-dev/openmc/pull/3430>`_)
- Avoid negative heating values during pair production and bremsstrahlung (`#3426 <https://github.com/openmc-dev/openmc/pull/3426>`_)
- Fix no serialization of periodic_surface_id bug (`#3421 <https://github.com/openmc-dev/openmc/pull/3421>`_)
- Update _get_start_data to always grab the beginning of timestep time (`#3414 <https://github.com/openmc-dev/openmc/pull/3414>`_)
- Fixed a bug in charged particle energy deposition (`#3416 <https://github.com/openmc-dev/openmc/pull/3416>`_)
- Fix bug where the same mesh is written multiple times to settings.xml (`#3418 <https://github.com/openmc-dev/openmc/pull/3418>`_)
- small typo - spelling of Debian (`#3411 <https://github.com/openmc-dev/openmc/pull/3411>`_)
- added test for dagmc geometry plot (`#3375 <https://github.com/openmc-dev/openmc/pull/3375>`_)
- Random Ray Misc Memory Error Fixes (`#3405 <https://github.com/openmc-dev/openmc/pull/3405>`_)
- added type hints to model file (`#3399 <https://github.com/openmc-dev/openmc/pull/3399>`_)
- Apply resolve paths to path values in ``config`` (`#3400 <https://github.com/openmc-dev/openmc/pull/3400>`_)
- Fixing an incorrect computation of CDF of bremsstrahlung photons (`#3396 <https://github.com/openmc-dev/openmc/pull/3396>`_)
- Fix weight modification for uniform source sampling (`#3395 <https://github.com/openmc-dev/openmc/pull/3395>`_)
- Updates to VTK data checks (`#3371 <https://github.com/openmc-dev/openmc/pull/3371>`_)
- Map Compton subshell data to atomic relaxation data (`#3392 <https://github.com/openmc-dev/openmc/pull/3392>`_)
- Skip atomic relaxation if binding energy is larger than photon energy (`#3391 <https://github.com/openmc-dev/openmc/pull/3391>`_)
- Fix extremely large yields from Bremsstrahlung (`#3386 <https://github.com/openmc-dev/openmc/pull/3386>`_)
- corrected tally name in D1S example (`#3383 <https://github.com/openmc-dev/openmc/pull/3383>`_)
- Install MCPL using same build type as OpenMC in CI (`#3388 <https://github.com/openmc-dev/openmc/pull/3388>`_)
- using reduce chain level to remove need for reduce chain (`#3377 <https://github.com/openmc-dev/openmc/pull/3377>`_)
- Fix negative distances from bins_crossed for CylindricalMesh (`#3370 <https://github.com/openmc-dev/openmc/pull/3370>`_)
- Add check for equal value bins in an EnergyFilter (`#3372 <https://github.com/openmc-dev/openmc/pull/3372>`_)
- Fix for Issue Loading MGXS Data Files with LLVM 20 or Newer (`#3368 <https://github.com/openmc-dev/openmc/pull/3368>`_)
- Report plot ID instead of index for unsupported plot types in random ray mode (`#3361 <https://github.com/openmc-dev/openmc/pull/3361>`_)
- Handle Missing Tags in Versioning by Setting Default to 0 (`#3359 <https://github.com/openmc-dev/openmc/pull/3359>`_)
- added kg units to doc string in results class (`#3358 <https://github.com/openmc-dev/openmc/pull/3358>`_)

View file

@ -7,6 +7,7 @@ Release Notes
.. toctree::
:maxdepth: 1
0.15.3
0.15.2
0.15.1
0.15.0

View file

@ -6,42 +6,189 @@ Decay Sources
Through the :ref:`depletion <usersguide_depletion>` capabilities in OpenMC, it
is possible to simulate radiation emitted from the decay of activated materials.
For fusion energy systems, this is commonly done using what is known as the
`rigorous 2-step <https://doi.org/10.1016/S0920-3796(02)00144-8>`_ (R2S) method.
In this method, a neutron transport calculation is used to determine the neutron
flux and reaction rates over a cell- or mesh-based spatial discretization of the
model. Then, the neutron flux in each discrete region is used to predict the
activated material composition using a depletion solver. Finally, a photon
transport calculation with a source based on the activity and energy spectrum of
the activated materials is used to determine a desired physical response (e.g.,
a dose rate) at one or more locations of interest.
For fusion energy systems, this is commonly done using either the `rigorous
2-step <https://doi.org/10.1016/S0920-3796(02)00144-8>`_ (R2S) method or the
`direct 1-step <https://doi.org/10.1016/S0920-3796(01)00188-0>`_ (D1S) method.
In the R2S method, a neutron transport calculation is used to determine the
neutron flux and reaction rates over a cell- or mesh-based spatial
discretization of the model. Then, the neutron flux in each discrete region is
used to predict the activated material composition using a depletion solver.
Finally, a photon transport calculation with a source based on the activity and
energy spectrum of the activated materials is used to determine a desired
physical response (e.g., a dose rate) at one or more locations of interest.
OpenMC includes automation for both the R2S and D1S methods as described in the
following sections.
Once a depletion simulation has been completed in OpenMC, the intrinsic decay
source can be determined as follows. First the activated material composition
can be determined using the :class:`openmc.deplete.Results` object. Indexing an
instance of this class with the timestep index returns a
:class:`~openmc.deplete.StepResult` object, which itself has a
:meth:`~openmc.deplete.StepResult.get_material` method. Once the activated
:class:`~openmc.Material` has been obtained, the
:meth:`~openmc.Material.get_decay_photon_energy` method will give the energy
spectrum of the decay photon source. The integral of the spectrum also indicates
the intensity of the source in units of [Bq]. Altogether, the workflow looks as
follows::
Rigorous 2-Step (R2S) Calculations
==================================
OpenMC includes an :class:`openmc.deplete.R2SManager` class that fully automates
cell- and mesh-based R2S calculations. Before we describe this class, it is
useful to understand the basic mechanics of how an R2S calculation works.
Generally, it involves the following steps:
1. The :meth:`openmc.deplete.get_microxs_and_flux` function is called to run a
neutron transport calculation that determines fluxes and microscopic cross
sections in each activation region.
2. The :class:`openmc.deplete.IndependentOperator` and
:class:`openmc.deplete.PredictorIntegrator` classes are used to carry out a
depletion (activation) calculation in order to determine predicted material
compositions based on a set of timesteps and source rates.
3. The activated material composition is determined using the
:class:`openmc.deplete.Results` class. Indexing an instance of this class
with the timestep index returns a :class:`~openmc.deplete.StepResult` object,
which itself has a :meth:`~openmc.deplete.StepResult.get_material` method
returning an activated material.
4. The :meth:`openmc.Material.get_decay_photon_energy` method is used to obtain
the energy spectrum of the decay photon source. The integral of the spectrum
also indicates the intensity of the source in units of [Bq].
5. A new photon source is defined using one of OpenMC's source classes with the
energy distribution set equal to the object returned by the
:meth:`openmc.Material.get_decay_photon_energy` method. The source is then
assigned to a photon :class:`~openmc.Model`.
6. A photon transport calculation is run with ``model.run()``.
Altogether, the workflow looks as follows::
# Run neutron transport calculation
fluxes, micros = openmc.deplete.get_microxs_and_flux(model, domains)
# Run activation calculation
op = openmc.deplete.IndependentOperator(mats, fluxes, micros)
timesteps = ...
source_rates = ...
integrator = openmc.deplete.Integrator(op, timesteps, source_rates)
integrator.integrate()
# Get decay photon source at last timestep
results = openmc.deplete.Results("depletion_results.h5")
# Get results at last timestep
step = results[-1]
# Get activated material composition for ID=1
activated_mat = step.get_material('1')
# Determine photon source
photon_energy = activated_mat.get_decay_photon_energy()
photon_source = openmc.IndependentSource(
space=...,
energy=photon_energy,
particle='photon',
strength=photon_energy.integral()
)
By default, the :meth:`~openmc.Material.get_decay_photon_energy` method will
eliminate spectral lines with very low intensity, but this behavior can be
configured with the ``clip_tolerance`` argument.
# Run photon transport calculation
model.settings.source = photon_source
model.run()
Note that by default, the :meth:`~openmc.Material.get_decay_photon_energy`
method will eliminate spectral lines with very low intensity, but this behavior
can be configured with the ``clip_tolerance`` argument.
Cell-based R2S
--------------
In practice, users do not need to manually go through each of the steps in an R2S
calculation described above. The :class:`~openmc.deplete.R2SManager` fully
automates the execution of neutron transport, depletion, decay source
generation, and photon transport. For a cell-based R2S calculation, once you
have a :class:`~openmc.Model` that has been defined, simply create an instance
of :class:`~openmc.deplete.R2SManager` by passing the model and a list of cells
to activate::
r2s = openmc.deplete.R2SManager(model, [cell1, cell2, cell3])
Note that the ``volume`` attribute must be set for any cell that is to be
activated. The :class:`~openmc.deplete.R2SManager` class allows you to
optionally specify a separate photon model; if not given as an argument, it will
create a shallow copy of the original neutron model (available as the
``neutron_model`` attribute) and store it in the ``photon_model`` attribute. We
can use this to define tallies specific to the photon model::
dose_tally = openmc.Tally()
...
r2s.photon_model.tallies = [dose_tally]
Next, define the timesteps and source rates for the activation calculation::
timesteps = [(3.0, 'd'), (5.0, 'h')]
source_rates = [1e12, 0.0]
In this case, the model is irradiated for 3 days with a source rate of
:math:`10^{12}` neutron/sec and then the source is turned off and the activated
materials are allowed to decay for 5 hours. These parameters should be passed to
the :meth:`~openmc.deplete.R2SManager.run` method to execute the full R2S
calculation. Before we can do that though, for a cell-based calculation, the one
other piece of information that is needed is bounding boxes of the activated
cells::
bounding_boxes = {
cell1.id: cell1.bounding_box,
cell2.id: cell2.bounding_box,
cell3.id: cell3.bounding_box
}
Note that calling the ``bounding_box`` attribute may not work for all
constructive solid geometry regions (for example, a cell that uses a
non-axis-aligned plane). In these cases, the bounding box will need to be
specified manually. Once you have a set of bounding boxes, the R2S calculation
can be run::
r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes)
If not specified otherwise, a photon transport calculation is run at each time
in the depletion schedule. That means in the case above, we would see three
photon transport calculations. To specify specific times at which photon
transport calculations should be run, pass the ``photon_time_indices`` argument.
For example, if we wanted to run a photon transport calculation only on the last
time (after the 5 hour decay), we would run::
r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes,
photon_time_indices=[2])
After an R2S calculation has been run, the :class:`~openmc.deplete.R2SManager`
instance will have a ``results`` dictionary that allows you to directly access
results from each of the steps. It will also write out all the output files into
a directory that is named "r2s_<timestamp>/". The ``output_dir`` argument to the
:meth:`~openmc.deplete.R2SManager.run` method enables you to override the
default output directory name if desired.
The :meth:`~openmc.deplete.R2SManager.run` method actually runs three
lower-level methods under the hood::
r2s.step1_neutron_transport(...)
r2s.step2_activation(...)
r2s.step3_photon_transport(...)
For users looking for more control over the calculation, these lower-level
methods can be used in lieu of the :meth:`openmc.deplete.R2SManager.run` method.
Mesh-based R2S
--------------
Executing a mesh-based R2S calculation looks nearly identical to the cell-based
R2S workflow described above. The only difference is that instead of passing a
list of cells to the ``domains`` argument of
:class:`~openmc.deplete.R2SManager`, you need to define a mesh object and pass
that instead. This might look like the following::
# Define a regular Cartesian mesh
mesh = openmc.RegularMesh()
mesh.lower_left = (-50., -50., 0.)
mesh.upper_right = (50., 50., 75.)
mesh.dimension = (10, 10, 5)
r2s = openmc.deplete.R2SManager(model, mesh)
Executing the R2S calculation is then performed by adding photon tallies and
calling the :meth:`~openmc.deplete.R2SManager.run` method with the appropriate
timesteps and source rates. Note that in this case we do not need to define cell
volumes or bounding boxes as is required for a cell-based R2S calculation.
Instead, during the neutron transport step, OpenMC will run a raytracing
calculation to determine material volume fractions within each mesh element
using the :meth:`openmc.MeshBase.material_volumes` method. Arguments to this
method can be customized via the ``mat_vol_kwargs`` argument to the
:meth:`~openmc.deplete.R2SManager.run` method. Most often, this would involve
customizing the number of rays traced to obtain better estimates of volumes. As
an example, if we wanted to run the raytracing calculation with 10 million rays,
we would run::
r2s.run(timesteps, source_rates, mat_vol_kwargs={'n_samples': 10_000_000})
Direct 1-Step (D1S) Calculations
================================

View file

@ -182,6 +182,8 @@ boundary condition.
Periodic boundary conditions can be applied to pairs of planar surfaces.
If there are only two periodic surfaces they will be matched automatically.
Otherwise it is necessary to specify pairs explicitly using the
:attr:`Surface.periodic_surface` attribute as in the following example::
@ -192,7 +194,7 @@ Otherwise it is necessary to specify pairs explicitly using the
Both rotational and translational periodic boundary conditions are specified in
the same fashion. If both planes have the same normal vector, a translational
periodicity is assumed; rotational periodicity is assumed otherwise. Currently,
only rotations about the :math:`z`-axis are supported.
rotations must be about the :math:`x`-, :math:`y`-, or :math:`z`-axis.
For a rotational periodic BC, the normal vectors of each surface must point
inwards---towards the valid geometry. For example, a :class:`XPlane` and
@ -530,6 +532,89 @@ UWUW and OpenMC material ID space will cause an error. To automatically resolve
these ID overlaps, ``auto_ids`` can be set to ``True`` to append the UWUW
material IDs to the OpenMC material ID space.
Material overrides and differentiation
--------------------------------------
Programmatic access to DAGMC cell information for material overrides
and differentiation requires synchronization of the DAGMC universe
representation across Python and C-API::
model.init_lib()
model.sync_dagmc_universes()
model.finalize_lib()
Upon completion of these steps, the :attr:`DAGMCUniverse.cells` attribute will
be populated with :class:`DAGMCCell` proxy objects that represent the cells
defined in the DAGMC model. The :class:`DAGMCCell` objects will have
:class:`openmc.Material`'s' applied according to the assignments upon
initialization of the model. These materials can be replaced in the same manner
as :class:`openmc.Cell` objects to override material assignments in the DAGMC
model.
Depletion with DAGMC geometry
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The synchronization of :class:`openmc.DAGMCUniverse`'s is important for
depletion calculations using DAGMC geometry when materials need to be
differentiated to perform material burnup independently in each DAGMC cell. See
:meth:`openmc.model.Model.differentiate_mats`.
Material overrides
~~~~~~~~~~~~~~~~~~
OpenMC supports overriding material assignments defined inside a DAGMC HDF5
model so that CAD-assigned materials can be replaced by :class:`openmc.Material`
objects. This is useful when the CAD geometry provides the shape but OpenMC
materials (specific nuclide content, densities, or depletion behavior) are
required.
Replacing materials by name
^^^^^^^^^^^^^^^^^^^^^^^^^^^
If a DAGMC file includes material name tags, you can replace all cells that
reference a particular name with an :class:`openmc.Material` using
:meth:`~openmc.DAGMCUniverse.replace_material_assignment`::
import openmc
dag_univ = openmc.DAGMCUniverse('dagmc.h5m')
fuel = openmc.Material(name='fuel')
fuel.add_nuclide('U235', 0.05)
fuel.add_nuclide('U238', 0.95)
fuel.set_density('g/cm3', 10.5)
dag_univ.replace_material_assignment('Fuel', fuel)
This lets you keep CAD geometry while adopting OpenMC material definitions.
Per-cell material overrides
^^^^^^^^^^^^^^^^^^^^^^^^^^^
To assign overrides without initializing :class:`openmc.Model`, the
:meth:`openmc.DAGMCUniverse.add_material_override` method can be used to assign
materials to particular DAGMC cells. The method accepts either an integer cell
ID::
dag_univ = openmc.DAGMCUniverse('dagmc.h5m')
enriched = openmc.Material(name='fuel_enriched')
enriched.add_nuclide('U235', 0.10)
enriched.add_nuclide('U238', 0.90)
enriched.set_density('g/cm3', 10.5)
dag_univ.add_material_override(1, enriched)
In the case that the :class:`openmc.DAGMCUniverse` has already been synchronized,
a :class:`openmc.DAGMCCell` object can also be provide to assign the material.
Overrides are written to the `<material_overrides>` element of the
:ref:`<dagmc_universe> <dagmc_element>` XML element so the C++ core can apply
them on initialization.
.. _Direct Accelerated Geometry Monte Carlo: https://svalinn.github.io/DAGMC/
.. _University of Wisconsin Unified Workflow: https://svalinn.github.io/DAGMC/usersguide/uw2.html

View file

@ -67,6 +67,23 @@ are needed to compute kinetics parameters in OpenMC:
Obtaining kinetics parameters
-----------------------------
The ``Model`` class can be used to automatically generate all IFP tallies using
the Python API with :attr:`openmc.Settings.ifp_n_generation` greater than 0 and
the :meth:`openmc.Model.add_ifp_kinetics_tallies` method::
model = openmc.Model(geometry, settings=settings)
model.add_kinetics_parameters_tallies(num_groups=6) # Add 6 precursor groups
Alternatively, each of the tallies can be manually defined using group-wise or
total :math:`\beta_{\text{eff}}` specified by providing a 6-group
:class:`openmc.DelayedGroupFilter`::
beta_tally = openmc.Tally(name="group-beta-score")
beta_tally.scores = ["ifp-beta-numerator"]
# Add DelayedGroupFilter to enable group-wise tallies
beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, 7)))]
Here is an example showing how to declare the three available IFP scores in a
single tally::
@ -95,6 +112,12 @@ for ``ifp-denominator``:
\beta_{\text{eff}} = \frac{S_{\text{ifp-beta-numerator}}}{S_{\text{ifp-denominator}}}
The kinetics parameters can be retrieved directly from a statepoint file using
the :meth:`openmc.StatePoint.ifp_results` method::
with openmc.StatePoint(output_path) as sp:
generation_time, beta_eff = sp.get_kinetics_parameters()
.. only:: html
.. rubric:: References
@ -107,4 +130,4 @@ for ``ifp-denominator``:
of the Iterated Fission Probability Method in OpenMC to Compute Adjoint-Weighted
Kinetics Parameters", International Conference on Mathematics and Computational
Methods Applied to Nuclear Science and Engineering (M&C 2025), Denver, April 27-30,
2025 (to be presented).
2025.

View file

@ -6,13 +6,14 @@ Geometry Visualization
.. currentmodule:: openmc
OpenMC is capable of producing two-dimensional slice plots of a geometry as well
as three-dimensional voxel plots using the geometry plotting :ref:`run mode
<usersguide_run_modes>`. The geometry plotting mode relies on the presence of a
:ref:`plots.xml <io_plots>` file that indicates what plots should be created. To
create this file, one needs to create one or more :class:`openmc.Plot`
instances, add them to a :class:`openmc.Plots` collection, and then use the
:class:`Plots.export_to_xml` method to write the ``plots.xml`` file.
OpenMC is capable of producing two-dimensional slice plots of a geometry,
three-dimensional voxel plots, and three-dimensional raytrace plots using the
geometry plotting :ref:`run mode <usersguide_run_modes>`. The geometry plotting
mode relies on the presence of a :ref:`plots.xml <io_plots>` file that indicates
what plots should be created. To create this file, one needs to create one or
more instances of the various plot classes described below, add them to a
:class:`openmc.Plots` collection, and then use the :class:`Plots.export_to_xml`
method to write the ``plots.xml`` file.
-----------
Slice Plots
@ -21,15 +22,14 @@ Slice Plots
.. image:: ../_images/atr.png
:width: 300px
By default, when an instance of :class:`openmc.Plot` is created, it indicates
that a 2D slice plot should be made. You can specify the origin of the plot
(:attr:`Plot.origin`), the width of the plot in each direction
(:attr:`Plot.width`), the number of pixels to use in each direction
(:attr:`Plot.pixels`), and the basis directions for the plot. For example, to
create a :math:`x` - :math:`z` plot centered at (5.0, 2.0, 3.0) with a width of
(50., 50.) and 400x400 pixels::
The :class:`openmc.SlicePlot` class indicates that a 2D slice plot should be
made. You can specify the origin of the plot (:attr:`SlicePlot.origin`), the
width of the plot in each direction (:attr:`SlicePlot.width`), the number of
pixels to use in each direction (:attr:`SlicePlot.pixels`), and the basis
directions for the plot. For example, to create a :math:`x` - :math:`z` plot
centered at (5.0, 2.0, 3.0) with a width of (50., 50.) and 400x400 pixels::
plot = openmc.Plot()
plot = openmc.SlicePlot()
plot.basis = 'xz'
plot.origin = (5.0, 2.0, 3.0)
plot.width = (50., 50.)
@ -47,7 +47,7 @@ that location.
By default, a unique color will be assigned to each cell in the geometry. If you
want your plot to be colored by material instead, change the
:attr:`Plot.color_by` attribute::
:attr:`SlicePlot.color_by` attribute::
plot.color_by = 'material'
@ -68,8 +68,8 @@ particular cells/materials should be given colors of your choosing::
Note that colors can be given as RGB tuples or by a string indicating a valid
`SVG color <https://www.w3.org/TR/SVG11/types.html#ColorKeywords>`_.
When you're done creating your :class:`openmc.Plot` instances, you need to then
assign them to a :class:`openmc.Plots` collection and export it to XML::
When you're done creating your :class:`openmc.SlicePlot` instances, you need to
then assign them to a :class:`openmc.Plots` collection and export it to XML::
plots = openmc.Plots([plot1, plot2, plot3])
plots.export_to_xml()
@ -97,13 +97,11 @@ Voxel Plots
.. image:: ../_images/3dba.png
:width: 200px
The :class:`openmc.Plot` class can also be told to generate a 3D voxel plot
instead of a 2D slice plot. Simply change the :attr:`Plot.type` attribute to
'voxel'. In this case, the :attr:`Plot.width` and :attr:`Plot.pixels` attributes
should be three items long, e.g.::
The :class:`openmc.VoxelPlot` class enables the generation of a 3D voxel plot
instead of a 2D slice plot. In this case, the :attr:`VoxelPlot.width` and
:attr:`VoxelPlot.pixels` attributes should be three items long, e.g.::
vox_plot = openmc.Plot()
vox_plot.type = 'voxel'
vox_plot = openmc.VoxelPlot()
vox_plot.width = (100., 100., 50.)
vox_plot.pixels = (400, 400, 200)

View file

@ -68,17 +68,17 @@ generation, and particle number of the desired particle. For example, to create
a track file for particle 4 of batch 1 and generation 2::
settings = openmc.Settings()
settings.track = (1, 2, 4)
settings.track = [(1, 2, 4)]
To specify multiple particles, the length of the iterable should be a multiple
of three, e.g., if we wanted particles 3 and 4 from batch 1 and generation 2::
To specify multiple particles, specify a list of tuples, e.g., if we wanted
particles 3 and 4 from batch 1 and generation 2::
settings.track = (1, 2, 3, 1, 2, 4)
settings.track = [(1, 2, 3), (1, 2, 4)]
After running OpenMC, the working directory will contain a file of the form
"track_(batch #)_(generation #)_(particle #).h5" for each particle tracked.
These track files can be converted into VTK poly data files with the
:class:`openmc.Tracks` class.
After running OpenMC (now, without the ``-t`` argument), the working directory
will contain a file named `tracks.h5`, which contains a collection of particle
tracks. These track files can be converted into VTK poly data files or
matplotlib plots with the :class:`openmc.Tracks` class.
----------------------
Source Site Processing

View file

@ -644,7 +644,8 @@ model to use these multigroup cross sections. An example is given below::
nparticles=2000,
overwrite_mgxs_library=False,
mgxs_path="mgxs.h5",
correction=None
correction=None,
source_energy=None
)
The most important parameter to set is the ``method`` parameter, which can be
@ -706,6 +707,31 @@ generation and use an existing library file.
with a :math:`\rho` default value of 1.0, which can be adjusted with the
``settings.random_ray['diagonal_stabilization_rho']`` parameter.
When generating MGXS data with either the ``stochastic_slab`` or
``infinite_medium`` methods, by default the simulation will use a uniform source
distribution spread evenly over all energy groups. This ensures that all energy
groups receive tallies and therefore produce non-zero total multigroup cross
sections. Additionally, the function will convert any sources in the model into
simplified spatial sources that retain the original energy distributions. If
sources are present, they will be used 99% of the time to sample source energies
during MGXS generation. The other 1% of the time, energies will be sampled
uniformly over all energy groups to ensure that all groups receive some tallies.
However, the user may wish to specify a different source energy spectrum (for
instance, if they are using a FileSource, such that the energy distribution
cannot be extracted from the python source object). This can be done by
providing a :class:`openmc.stats.Univariate` distribution as the
``source_energy`` parameter of the :meth:`openmc.Model.convert_to_multigroup`
method. If provided, it will override any sources present in the model and will
be used 99% of the time to sample source energies during MGXS generation. The
other 1% of the time, energies will be sampled uniformly over all energy groups
to ensure that all groups receive some tallies.
For instance, a D-D fusion simulation may involve a complex file source. In this
case, the user may wish to provide a discrete 2.45 MeV energy source
distribution for MGXS generation as::
source_energy = openmc.stats.delta_function(2.45e6)
Ultimately, the methods described above are all just approximations.
Approximations in the generated MGXS data will fundamentally limit the potential
accuracy of the random ray solver. However, the methods described above are all
@ -765,7 +791,7 @@ energy decomposition::
# Create a "tallies.xml" file for the MGXS Library
tallies = openmc.Tallies()
mgxs_lib.add_to_tallies_file(tallies, merge=True)
mgxs_lib.add_to_tallies(tallies, merge=True)
# Export
tallies.export_to_xml()
@ -1105,11 +1131,10 @@ given below:
tallies.export_to_xml()
# Create voxel plot
plot = openmc.Plot()
plot = openmc.VoxelPlot()
plot.origin = [0, 0, 0]
plot.width = [2*pitch, 2*pitch, 1]
plot.pixels = [1000, 1000, 1]
plot.type = 'voxel'
# Instantiate a Plots collection and export to XML
plots = openmc.Plots([plot])
@ -1189,11 +1214,10 @@ given below:
tallies.export_to_xml()
# Create voxel plot
plot = openmc.Plot()
plot = openmc.VoxelPlot()
plot.origin = [0, 0, 0]
plot.width = [2*pitch, 2*pitch, 1]
plot.pixels = [1000, 1000, 1]
plot.type = 'voxel'
# Instantiate a Plots collection and export to XML
plots = openmc.Plots([plot])

View file

@ -48,6 +48,7 @@ flags:
restart file
-s, --threads N Run with *N* OpenMP threads
-t, --track Write tracks for all particles (up to max_tracks)
-q, --verbosity V Set the output verbosity to *V*
-v, --version Show version information
-h, --help Show help message

View file

@ -756,6 +756,62 @@ instance, whereas the :meth:`openmc.Track.filter` method returns a new
track_files = [f"tracks_p{rank}.h5" for rank in range(32)]
openmc.Tracks.combine(track_files, "tracks.h5")
Collision Track File
---------------------
OpenMC can generate a collision track file that contains detailed collision
information (position, direction, energy, deposited energy, time, weight, cell
ID, material ID, universe ID, nuclide ZAID, particle type, particle delayed
group and particle ID) for each particle collision depending on user-defined
parameters. To invoke this feature, set the
:attr:`~openmc.Settings.collision_track` attribute as shown in this example::
settings.collision_track = {
"max_collisions": 300,
"reactions": ["(n,fission)", "(n,2n)"],
"material_ids": [1,2],
"nuclides": ["U238", "O16"],
"cell_ids": [5, 12]
}
In this example, collision track information is written to the
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.
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:
>>> data = openmc.read_collision_track_file('collision_track.h5')
>>> data.dtype
dtype([('r', [('x', '<f8'), ('y', '<f8'), ('z', '<f8')]),
('u', [('x', '<f8'), ('y', '<f8'), ('z', '<f8')]), ('E', '<f8'),
('dE', '<f8'), ('time', '<f8'), ('wgt', '<f8'), ('event_mt', '<i4'),
('delayed_group', '<i4'), ('cell_id', '<i4'), ('nuclide_id', '<i4'),
('material_id', '<i4'), ('universe_id', '<i4'), ('n_collision', '<i4'),
('particle', '<i4'), ('parent_id', '<i8'), ('progeny_id', '<i8')])
The full list of fields is as follows:
:r: Position (each direction in [cm])
:u: Direction
:E: Energy in [eV]
:dE: Energy deposited during collision in [eV]
:time: Time in [s]
:wgt: Weight of the particle
:event_mt: Reaction MT number
:delayed_group: Delayed group of the particle
:cell_id: Cell ID
:nuclide_id: Nuclide ID (10000×Z + 10×A + M)
:material_id: Material ID
:universe_id: Universe ID
:n_collision: Number of collision suffered by the particle
:particle: Particle type
:parent_id: Source particle ID
:progeny_id: Progeny ID
-----------------------
Restarting a Simulation
-----------------------

View file

@ -51,7 +51,7 @@ With the :class:`~openmc.WeightWindowGenerator` instance added to the
:attr:`~openmc.Settings`, the rest of the problem can be defined as normal. When
running, note that the second iteration and beyond may be several orders of
magnitude slower than the first. As the weight windows are applied in each
iteration, particles may be agressively split, resulting in a large number of
iteration, particles may be aggressively split, resulting in a large number of
secondary (split) particles being generated per initial source particle. This is
not necessarily a bad thing, as the split particles are much more efficient at
exploring low flux regions of phase space as compared to initial particles.
@ -133,8 +133,7 @@ random ray mode can be found in the :ref:`Random Ray User Guide <random_ray>`.
# we used for source region decomposition
wwg = openmc.WeightWindowGenerator(
method='fw_cadis',
mesh=mesh,
max_realizations=settings.batches
mesh=mesh
)
# Add generator to openmc.settings object
@ -162,7 +161,7 @@ solver, the Python input just needs to load the h5 file::
settings.weight_window_checkpoints = {'collision': True, 'surface': True}
settings.survival_biasing = False
settings.weight_windows = openmc.WeightWindowsList.from_hdf5('weight_windows.h5')
settings.weight_windows_file = "weight_windows.h5"
settings.weight_windows_on = True
The :class:`~openmc.WeightWindowGenerator` instance is not needed to load an

View file

@ -128,14 +128,14 @@ settings_file.export_to_xml()
# Exporting to OpenMC plots.xml file
###############################################################################
plot_xy = openmc.Plot(plot_id=1)
plot_xy = openmc.SlicePlot(plot_id=1)
plot_xy.filename = 'plot_xy'
plot_xy.origin = [0, 0, 0]
plot_xy.width = [6, 6]
plot_xy.pixels = [400, 400]
plot_xy.color_by = 'material'
plot_yz = openmc.Plot(plot_id=2)
plot_yz = openmc.SlicePlot(plot_id=2)
plot_yz.filename = 'plot_yz'
plot_yz.basis = 'yz'
plot_yz.origin = [0, 0, 0]

View file

@ -135,7 +135,7 @@ settings_file.export_to_xml()
# Exporting to OpenMC plots.xml file
###############################################################################
plot = openmc.Plot(plot_id=1)
plot = openmc.SlicePlot(plot_id=1)
plot.origin = [0, 0, 0]
plot.width = [4, 4]
plot.pixels = [400, 400]

View file

@ -128,7 +128,7 @@ settings_file.export_to_xml()
# Exporting to OpenMC plots.xml file
###############################################################################
plot = openmc.Plot(plot_id=1)
plot = openmc.SlicePlot(plot_id=1)
plot.origin = [0, 0, 0]
plot.width = [4, 4]
plot.pixels = [400, 400]

View file

@ -0,0 +1,101 @@
import matplotlib.pyplot as plt
import numpy as np
import openmc
###############################################################################
# Create materials for the problem
uo2 = openmc.Material(name="UO2 fuel at 2.4% wt enrichment")
uo2.set_density("g/cm3", 10.29769)
uo2.add_element("U", 1.0, enrichment=2.4)
uo2.add_element("O", 2.0)
helium = openmc.Material(name="Helium for gap")
helium.set_density("g/cm3", 0.001598)
helium.add_element("He", 2.4044e-4)
zircaloy = openmc.Material(name="Zircaloy 4")
zircaloy.set_density("g/cm3", 6.55)
zircaloy.add_element("Sn", 0.014, "wo")
zircaloy.add_element("Fe", 0.00165, "wo")
zircaloy.add_element("Cr", 0.001, "wo")
zircaloy.add_element("Zr", 0.98335, "wo")
borated_water = openmc.Material(name="Borated water")
borated_water.set_density("g/cm3", 0.740582)
borated_water.add_element("B", 2.0e-4) # 3x the original pincell
borated_water.add_element("H", 5.0e-2)
borated_water.add_element("O", 2.4e-2)
borated_water.add_s_alpha_beta("c_H_in_H2O")
###############################################################################
# Define problem geometry
# Create cylindrical surfaces
fuel_or = openmc.ZCylinder(r=0.39218, name="Fuel OR")
clad_ir = openmc.ZCylinder(r=0.40005, name="Clad IR")
clad_or = openmc.ZCylinder(r=0.45720, name="Clad OR")
# Create a region represented as the inside of a rectangular prism
pitch = 1.25984
box = openmc.model.RectangularPrism(pitch, pitch, boundary_type="reflective")
# Create cells, mapping materials to regions
fuel = openmc.Cell(fill=uo2, region=-fuel_or)
gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
water = openmc.Cell(fill=borated_water, region=+clad_or & -box)
# Create a model and assign geometry
model = openmc.Model()
model.geometry = openmc.Geometry([fuel, gap, clad, water])
###############################################################################
# Define problem settings
# Set the mode
model.settings.run_mode = "fixed source"
# Indicate how many batches and particles to run
model.settings.batches = 10
model.settings.particles = 10000
# Set time cutoff (we only care about t < 100 seconds, see tally below)
model.settings.cutoff = {"time_neutron": 100}
# Create the neutron pulse source (by default, isotropic direction, t=0)
space = openmc.stats.Point() # At the origin (0, 0, 0)
energy = openmc.stats.delta_function(14.1e6) # At 14.1 MeV
model.settings.source = openmc.IndependentSource(space=space, energy=energy)
###############################################################################
# Define tallies
# Create time filter
t_grid = np.insert(np.logspace(-6, 2, 100), 0, 0.0)
time_filter = openmc.TimeFilter(t_grid)
# Tally for total neutron density in time
density_tally = openmc.Tally(name="Density")
density_tally.filters = [time_filter]
density_tally.scores = ["inverse-velocity"]
# Add tallies to model
model.tallies = openmc.Tallies([density_tally])
# Run the model
model.run(apply_tally_results=True)
# Bin-averaged result
density_mean = density_tally.mean.ravel() / np.diff(t_grid)
# Plot particle density versus time
fig, ax = plt.subplots()
ax.stairs(density_mean, t_grid)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("Time [s]")
ax.set_ylabel("Total density")
ax.grid()
plt.show()

View file

@ -192,11 +192,10 @@ tallies.export_to_xml()
# Exporting to OpenMC plots.xml file
###############################################################################
plot = openmc.Plot()
plot = openmc.VoxelPlot()
plot.origin = [0, 0, 0]
plot.width = [pitch, pitch, pitch]
plot.pixels = [1000, 1000, 1]
plot.type = 'voxel'
# Instantiate a Plots collection and export to XML
plots = openmc.Plots([plot])

View file

@ -20,6 +20,8 @@ extern vector<SourceSite> source_bank;
extern SharedArray<SourceSite> surf_source_bank;
extern SharedArray<CollisionTrackSite> collision_track_bank;
extern SharedArray<SourceSite> fission_bank;
extern vector<vector<int>> ifp_source_delayed_group_bank;

105
include/openmc/bank_io.h Normal file
View file

@ -0,0 +1,105 @@
#ifndef OPENMC_BANK_IO_H
#define OPENMC_BANK_IO_H
#include "hdf5.h"
#include "openmc/message_passing.h"
#include "openmc/span.h"
#include "openmc/vector.h"
#include <algorithm>
#ifdef OPENMC_MPI
#include <mpi.h>
#endif
namespace openmc {
template<typename SiteType>
void write_bank_dataset(
const char* dataset_name, hid_t group_id, span<SiteType> bank,
const vector<int64_t>& bank_index, hid_t membanktype, hid_t filebanktype
#ifdef OPENMC_MPI
,
MPI_Datatype mpi_dtype
#endif
)
{
int64_t dims_size = bank_index.back();
int64_t count_size = bank_index[mpi::rank + 1] - bank_index[mpi::rank];
#ifdef PHDF5
hsize_t dims[] {static_cast<hsize_t>(dims_size)};
hid_t dspace = H5Screate_simple(1, dims, nullptr);
hid_t dset = H5Dcreate(group_id, dataset_name, filebanktype, dspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
hsize_t count[] {static_cast<hsize_t>(count_size)};
hid_t memspace = H5Screate_simple(1, count, nullptr);
hsize_t start[] {static_cast<hsize_t>(bank_index[mpi::rank])};
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
H5Dwrite(dset, membanktype, memspace, dspace, plist, bank.data());
H5Sclose(dspace);
H5Sclose(memspace);
H5Dclose(dset);
H5Pclose(plist);
#else
if (mpi::master) {
hsize_t dims[] {static_cast<hsize_t>(dims_size)};
hid_t dspace = H5Screate_simple(1, dims, nullptr);
hid_t dset = H5Dcreate(group_id, dataset_name, filebanktype, dspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
#ifdef OPENMC_MPI
vector<SiteType> temp_bank {bank.begin(), bank.end()};
#endif
for (int i = 0; i < mpi::n_procs; ++i) {
hsize_t count[] {static_cast<hsize_t>(bank_index[i + 1] - bank_index[i])};
hid_t memspace = H5Screate_simple(1, count, nullptr);
#ifdef OPENMC_MPI
if (i > 0) {
MPI_Recv(bank.data(), count[0], mpi_dtype, i, i, mpi::intracomm,
MPI_STATUS_IGNORE);
}
#endif
hid_t dspace_rank = H5Dget_space(dset);
hsize_t start[] {static_cast<hsize_t>(bank_index[i])};
H5Sselect_hyperslab(
dspace_rank, H5S_SELECT_SET, start, nullptr, count, nullptr);
H5Dwrite(
dset, membanktype, memspace, dspace_rank, H5P_DEFAULT, bank.data());
H5Sclose(memspace);
H5Sclose(dspace_rank);
}
H5Dclose(dset);
#ifdef OPENMC_MPI
std::copy(temp_bank.begin(), temp_bank.end(), bank.begin());
#endif
}
#ifdef OPENMC_MPI
else {
if (!bank.empty()) {
MPI_Send(
bank.data(), bank.size(), mpi_dtype, 0, mpi::rank, mpi::intracomm);
}
}
#endif
#endif
}
} // namespace openmc
#endif // OPENMC_BANK_IO_H

View file

@ -138,18 +138,26 @@ protected:
//==============================================================================
//! A BC that rotates particles about a global axis.
//
//! Currently only rotations about the z-axis are supported.
//! Only rotations about the x, y, and z axes are supported.
//==============================================================================
class RotationalPeriodicBC : public PeriodicBC {
public:
RotationalPeriodicBC(int i_surf, int j_surf);
enum PeriodicAxis { x, y, z };
RotationalPeriodicBC(int i_surf, int j_surf, PeriodicAxis axis);
double compute_periodic_rotation(
double rise_1, double run_1, double rise_2, double run_2) const;
void handle_particle(Particle& p, const Surface& surf) const override;
protected:
//! Angle about the axis by which particle coordinates will be rotated
double angle_;
//! Ensure that choice of axes is right handed. axis_1_idx_ corresponds to the
//! independent axis and axis_2_idx_ corresponds to the dependent axis in the
//! 2D plane perpendicular to the planes' axis of rotation
int zero_axis_idx_;
int axis_1_idx_;
int axis_2_idx_;
};
} // namespace openmc

View file

@ -17,6 +17,8 @@ int openmc_cell_get_fill(
int openmc_cell_get_id(int32_t index, int32_t* id);
int openmc_cell_get_temperature(
int32_t index, const int32_t* instance, double* T);
int openmc_cell_get_density(
int32_t index, const int32_t* instance, double* rho);
int openmc_cell_get_translation(int32_t index, double xyz[]);
int openmc_cell_get_rotation(int32_t index, double rot[], size_t* n);
int openmc_cell_get_name(int32_t index, const char** name);
@ -27,6 +29,8 @@ int openmc_cell_set_fill(
int openmc_cell_set_id(int32_t index, int32_t id);
int openmc_cell_set_temperature(
int32_t index, double T, const int32_t* instance, bool set_contained = false);
int openmc_cell_set_density(int32_t index, double rho, const int32_t* instance,
bool set_contained = false);
int openmc_cell_set_translation(int32_t index, const double xyz[]);
int openmc_cell_set_rotation(int32_t index, const double rot[], size_t rot_len);
int openmc_dagmc_universe_get_cell_ids(

View file

@ -216,6 +216,18 @@ public:
//! \return Temperature in [K]
double temperature(int32_t instance = -1) const;
//! Get the density multiplier of a cell instance
//! \param[in] instance Instance index. If -1 is given, the density multiplier
//! for the first instance is returned.
//! \return Density multiplier
double density_mult(int32_t instance = -1) const;
//! Get the density of a cell instance in g/cm3
//! \param[in] instance Instance index. If -1 is given, the density
//! for the first instance is returned.
//! \return Density in [g/cm3]
double density(int32_t instance = -1) const;
//! Set the temperature of a cell instance
//! \param[in] T Temperature in [K]
//! \param[in] instance Instance index. If -1 is given, the temperature for
@ -226,6 +238,16 @@ public:
void set_temperature(
double T, int32_t instance = -1, bool set_contained = false);
//! Set the density of a cell instance
//! \param[in] density Density [g/cm3]
//! \param[in] instance Instance index. If -1 is given, the density
//! for all instances is set.
//! \param[in] set_contained If this cell is not filled with a material,
//! collect all contained cells with material fills and set their
//! densities.
void set_density(
double density, int32_t instance = -1, bool set_contained = false);
int32_t n_instances() const;
//! Set the rotation matrix of a cell instance
@ -341,6 +363,9 @@ public:
//! T. The units are sqrt(eV).
vector<double> sqrtkT_;
//! \brief Unitless density multiplier(s) within this cell.
vector<double> density_mult_;
//! \brief Neighboring cells in the same universe.
NeighborList neighbors_;

View file

@ -0,0 +1,23 @@
#ifndef OPENMC_COLLISION_TRACK_H
#define OPENMC_COLLISION_TRACK_H
#include <string>
namespace openmc {
class Particle;
//! Reserve space in the collision track bank according to user settings.
void collision_track_reserve_bank();
//! Write collision track data to disk when the bank is full or the batch ends.
void collision_track_flush_bank();
//! Record the current particle as a collision-track entry when applicable.
//!
//! \param particle Particle whose collision should be recorded if eligible
void collision_track_record(Particle& particle);
} // namespace openmc
#endif // OPENMC_COLLISION_TRACK_H

View file

@ -28,12 +28,13 @@ constexpr int HDF5_VERSION[] {3, 0};
constexpr array<int, 2> VERSION_STATEPOINT {18, 1};
constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
constexpr array<int, 2> VERSION_TRACK {3, 0};
constexpr array<int, 2> VERSION_SUMMARY {6, 0};
constexpr array<int, 2> VERSION_SUMMARY {6, 1};
constexpr array<int, 2> VERSION_VOLUME {1, 0};
constexpr array<int, 2> VERSION_VOXEL {2, 0};
constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
constexpr array<int, 2> VERSION_PROPERTIES {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, 0};
// ============================================================================
// ADJUSTABLE PARAMETERS
@ -68,6 +69,11 @@ constexpr double MIN_HITS_PER_BATCH {1.5};
// prevent extremely large adjoint source terms from being generated.
constexpr double ZERO_FLUX_CUTOFF {1e-22};
// The minimum macroscopic cross section value considered non-void for the
// random ray solver. Materials with any group with a cross section below this
// value will be converted to pure void.
constexpr double MINIMUM_MACRO_XS {1e-6};
// ============================================================================
// MATH AND PHYSICAL CONSTANTS
@ -286,7 +292,7 @@ enum class MgxsType {
// ============================================================================
// TALLY-RELATED CONSTANTS
enum class TallyResult { VALUE, SUM, SUM_SQ, SIZE };
enum class TallyResult { VALUE, SUM, SUM_SQ, SUM_THIRD, SUM_FOURTH };
enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT };

View file

@ -51,6 +51,8 @@ public:
Distribution* phi() const { return phi_.get(); }
private:
Direction v_ref_ {1.0, 0.0, 0.0}; //!< reference direction
Direction w_ref_;
UPtrDist mu_; //!< Distribution of polar angle
UPtrDist phi_; //!< Distribution of azimuthal angle
};

View file

@ -37,6 +37,12 @@ void adjust_indices();
void assign_temperatures();
//==============================================================================
//! Finalize densities (compute density multipliers).
//==============================================================================
void finalize_cell_densities();
//==============================================================================
//! \brief Obtain a list of temperatures that each nuclide/thermal scattering
//! table appears at in the model. Later, this list is used to determine the

View file

@ -100,8 +100,8 @@ void read_llong(hid_t obj_id, const char* name, long long* buffer, bool indep);
void read_string(
hid_t obj_id, const char* name, size_t slen, char* buffer, bool indep);
void read_tally_results(
hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results);
void read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
hsize_t n_results, double* results);
void write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims,
const char* name, const double* buffer);
void write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims,
@ -114,9 +114,9 @@ void write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
void write_llong(hid_t group_id, int ndim, const hsize_t* dims,
const char* name, const long long* buffer, bool indep);
void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
const char* name, char const* buffer, bool indep);
void write_tally_results(
hid_t group_id, hsize_t n_filter, hsize_t n_score, const double* results);
const char* name, const char* buffer, bool indep);
void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
hsize_t n_results, const double* results);
} // extern "C"
//==============================================================================

View file

@ -68,15 +68,14 @@ vector<T> _ifp(const T& value, const vector<T>& data)
//!
//! Add the IFP information in the IFP banks using the same index
//! as the one used to append the fission site to the fission bank.
//! The information stored are the delayed group number and lifetime
//! of the neutron that created the fission event.
//! Multithreading protection is guaranteed by the index returned by the
//! thread_safe_append call in physics.cpp.
//!
//! Needs to be done after the delayed group is found.
//!
//! \param[in] p Particle
//! \param[in] site Fission site
//! \param[in] idx Bank index from the thread_safe_append call in physics.cpp
void ifp(const Particle& p, const SourceSite& site, int64_t idx);
void ifp(const Particle& p, int64_t idx);
//! Resize the IFP banks used in the simulation
void resize_simulation_ifp_banks();

View file

@ -99,6 +99,13 @@ public:
//----------------------------------------------------------------------------
// Accessors
//! Get the atom density in [atom/b-cm]
//! \return Density in [atom/b-cm]
double atom_density(int32_t i, double rho_multiplier = 1.0) const
{
return atom_density_(i) * rho_multiplier;
}
//! Get density in [atom/b-cm]
//! \return Density in [atom/b-cm]
double density() const { return density_; }

View file

@ -38,6 +38,21 @@ vector<SourceSite> mcpl_source_sites(std::string path);
void write_mcpl_source_point(const char* filename, span<SourceSite> source_bank,
const vector<int64_t>& bank_index);
//! Write an MCPL collision track file
//!
//! This function writes collision track data to an MCPL file. Additional
//! collision-specific metadata (such as energy deposition, material info, etc.)
//! is stored in the file header as blob data.
//!
//! \param[in] filename Path to MCPL file
//! \param[in] collision_track_bank Vector of CollisionTrackSites to write to
//! file for this MPI rank.
//! \param[in] bank_index Pointer to vector of site index ranges over all
//! MPI ranks.
void write_mcpl_collision_track(const char* filename,
span<CollisionTrackSite> collision_track_bank,
const vector<int64_t>& bank_index);
//! Check if MCPL functionality is available
bool is_mcpl_interface_available();

View file

@ -132,8 +132,14 @@ public:
// Constructors and destructor
Mesh() = default;
Mesh(pugi::xml_node node);
Mesh(hid_t group);
virtual ~Mesh() = default;
// Factory method for creating meshes from either an XML node or HDF5 group
template<typename T>
static const std::unique_ptr<Mesh>& create(
T dataset, const std::string& mesh_type, const std::string& mesh_library);
// Methods
//! Perform any preparation needed to support point location within the mesh
virtual void prepare_for_point_location() {};
@ -258,6 +264,7 @@ class StructuredMesh : public Mesh {
public:
StructuredMesh() = default;
StructuredMesh(pugi::xml_node node) : Mesh {node} {};
StructuredMesh(hid_t group) : Mesh {group} {};
virtual ~StructuredMesh() = default;
using MeshIndex = std::array<int, 3>;
@ -423,6 +430,7 @@ class PeriodicStructuredMesh : public StructuredMesh {
public:
PeriodicStructuredMesh() = default;
PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {};
PeriodicStructuredMesh(hid_t group) : StructuredMesh {group} {};
Position local_coords(const Position& r) const override
{
@ -442,6 +450,7 @@ public:
// Constructors
RegularMesh() = default;
RegularMesh(pugi::xml_node node);
RegularMesh(hid_t group);
// Overridden methods
int get_index_in_direction(double r, int i) const override;
@ -481,6 +490,8 @@ public:
//! Return the volume for a given mesh index
double volume(const MeshIndex& ijk) const override;
int set_grid();
// Data members
double volume_frac_; //!< Volume fraction of each mesh element
double element_volume_; //!< Volume of each mesh element
@ -492,6 +503,7 @@ public:
// Constructors
RectilinearMesh() = default;
RectilinearMesh(pugi::xml_node node);
RectilinearMesh(hid_t group);
// Overridden methods
int get_index_in_direction(double r, int i) const override;
@ -534,6 +546,7 @@ public:
// Constructors
CylindricalMesh() = default;
CylindricalMesh(pugi::xml_node node);
CylindricalMesh(hid_t group);
// Overridden methods
virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
@ -598,6 +611,7 @@ public:
// Constructors
SphericalMesh() = default;
SphericalMesh(pugi::xml_node node);
SphericalMesh(hid_t group);
// Overridden methods
virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
@ -666,9 +680,9 @@ class UnstructuredMesh : public Mesh {
public:
// Constructors
UnstructuredMesh() {};
UnstructuredMesh() { n_dimension_ = 3; };
UnstructuredMesh(pugi::xml_node node);
UnstructuredMesh(const std::string& filename);
UnstructuredMesh(hid_t group);
static const std::string mesh_type;
virtual std::string get_mesh_type() const override;
@ -775,6 +789,7 @@ public:
// Constructors
MOABMesh() = default;
MOABMesh(pugi::xml_node);
MOABMesh(hid_t group);
MOABMesh(const std::string& filename, double length_multiplier = 1.0);
MOABMesh(std::shared_ptr<moab::Interface> external_mbi);
@ -944,6 +959,7 @@ class LibMesh : public UnstructuredMesh {
public:
// Constructors
LibMesh(pugi::xml_node node);
LibMesh(hid_t group);
LibMesh(const std::string& filename, double length_multiplier = 1.0);
LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0);
@ -991,25 +1007,26 @@ public:
libMesh::MeshBase* mesh_ptr() const { return m_; };
protected:
// Methods
//! Translate a bin value to an element reference
virtual const libMesh::Elem& get_element_from_bin(int bin) const;
//! Translate an element pointer to a bin index
virtual int get_bin_from_element(const libMesh::Elem* elem) const;
libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set
//!< during intialization
private:
void initialize() override;
void set_mesh_pointer_from_filename(const std::string& filename);
void build_eqn_sys();
// Methods
//! Translate a bin value to an element reference
const libMesh::Elem& get_element_from_bin(int bin) const;
//! Translate an element pointer to a bin index
int get_bin_from_element(const libMesh::Elem* elem) const;
// Data members
unique_ptr<libMesh::MeshBase> unique_m_ =
nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is
//!< created inside OpenMC
libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set
//!< during intialization
vector<unique_ptr<libMesh::PointLocatorBase>>
pl_; //!< per-thread point locators
unique_ptr<libMesh::EquationSystems>
@ -1023,8 +1040,34 @@ private:
libMesh::BoundingBox bbox_; //!< bounding box of the mesh
libMesh::dof_id_type
first_element_id_; //!< id of the first element in the mesh
};
class AdaptiveLibMesh : public LibMesh {
public:
// Constructor
AdaptiveLibMesh(
libMesh::MeshBase& input_mesh, double length_multiplier = 1.0);
// Overridden methods
int n_bins() const override;
void add_score(const std::string& var_name) override;
void set_score_data(const std::string& var_name, const vector<double>& values,
const vector<double>& std_dev) override;
void write(const std::string& filename) const override;
protected:
// Overridden methods
int get_bin_from_element(const libMesh::Elem* elem) const override;
const libMesh::Elem& get_element_from_bin(int bin) const override;
private:
// Data members
const libMesh::dof_id_type num_active_; //!< cached number of active elements
const bool adaptive_; //!< whether this mesh has adaptivity enabled or not
std::vector<libMesh::dof_id_type>
bin_to_elem_map_; //!< mapping bin indices to dof indices for active
//!< elements
@ -1043,6 +1086,11 @@ private:
//! \param[in] root XML node
void read_meshes(pugi::xml_node root);
//! Read meshes from an HDF5 file
//
//! \param[in] group HDF5 group ("meshes" group)
void read_meshes(hid_t group);
//! Write mesh data to an HDF5 group
//
//! \param[in] group HDF5 group

View file

@ -18,6 +18,7 @@ extern bool master;
#ifdef OPENMC_MPI
extern MPI_Datatype source_site;
extern MPI_Datatype collision_track_site;
extern MPI_Comm intracomm;
#endif

View file

@ -164,8 +164,8 @@ namespace data {
// Minimum/maximum transport energy for each particle type. Order corresponds to
// that of the ParticleType enum
extern array<double, 2> energy_min;
extern array<double, 2> energy_max;
extern array<double, 4> energy_min;
extern array<double, 4> energy_max;
//! Minimum temperature in [K] that nuclide data is available at
extern double temperature_min;

View file

@ -56,6 +56,25 @@ struct SourceSite {
int64_t progeny_id;
};
struct CollisionTrackSite {
Position r;
Direction u;
double E;
double dE;
double time {0.0};
double wgt {1.0};
int event_mt {0};
int delayed_group {0};
int cell_id {0};
int nuclide_id;
int material_id {0};
int universe_id {0};
int n_collision {0};
ParticleType particle;
int64_t parent_id;
int64_t progeny_id;
};
//! State of a particle used for particle track files
struct TrackState {
Position r; //!< Position in [cm]
@ -154,9 +173,10 @@ struct NuclideMicroXS {
// Energy and temperature last used to evaluate these cross sections. If
// these values have changed, then the cross sections must be re-evaluated.
double last_E {0.0}; //!< Last evaluated energy
double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant
//!< * temperature (eV))
double last_E {0.0}; //!< Last evaluated energy
double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant
//!< * temperature (eV))
double ncrystal_xs {-1.0}; //!< NCrystal cross section
};
//==============================================================================
@ -377,47 +397,25 @@ public:
#ifdef OPENMC_DAGMC_ENABLED
// DagMC state variables
moab::DagMC::RayHistory& history()
{
return history_;
}
Direction& last_dir()
{
return last_dir_;
}
moab::DagMC::RayHistory& history() { return history_; }
Direction& last_dir() { return last_dir_; }
#endif
// material of current and last cell
int& material()
{
return material_;
}
const int& material() const
{
return material_;
}
int& material_last()
{
return material_last_;
}
const int& material_last() const
{
return material_last_;
}
int& material() { return material_; }
const int& material() const { return material_; }
int& material_last() { return material_last_; }
const int& material_last() const { return material_last_; }
// temperature of current and last cell
double& sqrtkT()
{
return sqrtkT_;
}
const double& sqrtkT() const
{
return sqrtkT_;
}
double& sqrtkT_last()
{
return sqrtkT_last_;
}
double& sqrtkT() { return sqrtkT_; }
const double& sqrtkT() const { return sqrtkT_; }
double& sqrtkT_last() { return sqrtkT_last_; }
// density multiplier of the current and last cell
double& density_mult() { return density_mult_; }
const double& density_mult() const { return density_mult_; }
double& density_mult_last() { return density_mult_last_; }
private:
int64_t id_ {-1}; //!< Unique ID
@ -447,6 +445,9 @@ private:
double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV
double sqrtkT_last_ {0.0}; //!< last temperature
double density_mult_ {1.0}; //!< density multiplier
double density_mult_last_ {1.0}; //!< last density multiplier
double collision_distance_ {INFTY};
#ifdef OPENMC_DAGMC_ENABLED
@ -653,6 +654,7 @@ public:
int& event_mt() { return event_mt_; } // MT number of collision
const int& event_mt() const { return event_mt_; }
int& delayed_group() { return delayed_group_; } // delayed group
const int& delayed_group() const { return delayed_group_; }
const int& parent_nuclide() const { return parent_nuclide_; }
int& parent_nuclide() { return parent_nuclide_; } // Parent nuclide

View file

@ -10,13 +10,6 @@
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
// Monoatomic ideal-gas scattering treatment threshold
constexpr double FREE_GAS_THRESHOLD {400.0};
//==============================================================================
// Non-member functions
//==============================================================================

View file

@ -27,8 +27,9 @@ public:
//----------------------------------------------------------------------------
// Methods
virtual void update_neutron_source(double k_eff);
double compute_k_eff(double k_eff_old) const;
virtual void update_single_neutron_source(SourceRegionHandle& srh);
virtual void update_all_neutron_sources();
void compute_k_eff();
virtual void normalize_scalar_flux_and_volumes(
double total_active_distance_per_iteration);
@ -41,7 +42,7 @@ public:
void output_to_vtk() const;
void convert_external_sources();
void count_external_source_regions();
void set_adjoint_sources(const vector<double>& forward_flux);
void set_adjoint_sources();
void flux_swap();
virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const;
double compute_fixed_source_normalization_factor() const;
@ -54,9 +55,8 @@ public:
bool is_target_void);
void apply_mesh_to_cell_and_children(int32_t i_cell, int32_t mesh_idx,
int32_t target_material_id, bool is_target_void);
void prepare_base_source_regions();
SourceRegionHandle get_subdivided_source_region_handle(
int64_t sr, int mesh_bin, Position r, double dist, Direction u);
SourceRegionKey sr_key, Position r, Direction u);
void finalize_discovered_source_regions();
void apply_transport_stabilization();
int64_t n_source_regions() const
@ -67,6 +67,10 @@ public:
{
return source_regions_.n_source_regions() * negroups_;
}
int64_t lookup_base_source_region_idx(const GeometryState& p) const;
SourceRegionKey lookup_source_region_key(const GeometryState& p) const;
int64_t lookup_mesh_bin(int64_t sr, Position r) const;
int lookup_mesh_idx(int64_t sr) const;
//----------------------------------------------------------------------------
// Static Data members
@ -86,6 +90,7 @@ public:
//----------------------------------------------------------------------------
// Public Data members
double k_eff_ {1.0}; // Eigenvalue
bool mapped_all_tallies_ {false}; // If all source regions have been visited
int64_t n_external_source_regions_ {0}; // Total number of source regions with
@ -110,14 +115,6 @@ public:
// The abstract container holding all source region-specific data
SourceRegionContainer source_regions_;
// Base source region container. When source region subdivision via mesh
// is in use, this container holds the original (non-subdivided) material
// filled cell instance source regions. These are useful as they can be
// initialized with external source and mesh domain information ahead of time.
// Then, dynamically discovered source regions can be initialized by cloning
// their base region.
SourceRegionContainer base_source_regions_;
// Parallel hash map holding all source regions discovered during
// a single iteration. This is a threadsafe data structure that is cleaned
// out after each iteration and stored in the "source_regions_" container.
@ -134,8 +131,17 @@ public:
// Map that relates a SourceRegionKey to the external source index. This map
// is used to check if there are any point sources within a subdivided source
// region at the time it is discovered.
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
point_source_map_;
std::unordered_map<SourceRegionKey, vector<int>, SourceRegionKey::HashFunctor>
external_point_source_map_;
// Map that relates a base source region index to the external source index.
// This map is used to check if there are any volumetric sources within a
// subdivided source region at the time it is discovered.
std::unordered_map<int64_t, vector<int>> external_volumetric_source_map_;
// Map that relates a base source region index to a mesh index. This map
// is used to check which subdivision mesh is present in a source region.
std::unordered_map<int64_t, int> mesh_map_;
// If transport corrected MGXS data is being used, there may be negative
// in-group scattering cross sections that can result in instability in MOC
@ -147,12 +153,11 @@ protected:
//----------------------------------------------------------------------------
// Methods
void apply_external_source_to_source_region(
Discrete* discrete, double strength_factor, SourceRegionHandle& srh);
void apply_external_source_to_cell_instances(int32_t i_cell,
Discrete* discrete, double strength_factor, int target_material_id,
const vector<int32_t>& instances);
void apply_external_source_to_cell_and_children(int32_t i_cell,
Discrete* discrete, double strength_factor, int32_t target_material_id);
int src_idx, SourceRegionHandle& srh);
void apply_external_source_to_cell_instances(int32_t i_cell, int src_idx,
int target_material_id, const vector<int32_t>& instances);
void apply_external_source_to_cell_and_children(
int32_t i_cell, int src_idx, int32_t target_material_id);
virtual void set_flux_to_flux_plus_source(int64_t sr, double volume, int g);
void set_flux_to_source(int64_t sr, int g);
virtual void set_flux_to_old_flux(int64_t sr, int g);
@ -165,6 +170,9 @@ protected:
simulation_volume_; // Total physical volume of the simulation domain, as
// defined by the 3D box of the random ray source
double
fission_rate_; // The system's fission rate (per cm^3), in eigenvalue mode
// Volumes for each tally and bin/score combination. This intermediate data
// structure is used when tallying quantities that must be normalized by
// volume (i.e., flux). The vector is index by tally index, while the inner 2D

View file

@ -20,7 +20,7 @@ class LinearSourceDomain : public FlatSourceDomain {
public:
//----------------------------------------------------------------------------
// Methods
void update_neutron_source(double k_eff) override;
void update_single_neutron_source(SourceRegionHandle& srh) override;
void normalize_scalar_flux_and_volumes(
double total_active_distance_per_iteration) override;

View file

@ -48,7 +48,6 @@ public:
static double distance_active_; // Active ray length
static unique_ptr<Source> ray_source_; // Starting source for ray sampling
static RandomRaySourceShape source_shape_; // Flag for linear source
static bool mesh_subdivision_enabled_; // Flag for mesh subdivision
static RandomRaySampleMethod sample_method_; // Flag for sampling method
//----------------------------------------------------------------------------

View file

@ -21,11 +21,7 @@ public:
// Methods
void compute_segment_correction_factors();
void apply_fixed_sources_and_mesh_domains();
void prepare_fixed_sources_adjoint(vector<double>& forward_flux,
SourceRegionContainer& forward_source_regions,
SourceRegionContainer& forward_base_source_regions,
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>&
forward_source_region_map);
void prepare_fixed_sources_adjoint();
void simulate();
void output_simulation_results() const;
void instability_check(
@ -45,9 +41,6 @@ private:
// Contains all flat source region data
unique_ptr<FlatSourceDomain> domain_;
// Random ray eigenvalue
double k_eff_ {1.0};
// Tracks the average FSR miss rate for analysis and reporting
double avg_miss_rate_ {0.0};

View file

@ -308,7 +308,6 @@ public:
//----------------------------------------------------------------------------
// Constructors
SourceRegion(int negroups, bool is_linear);
SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr);
SourceRegion() = default;
//----------------------------------------------------------------------------

View file

@ -32,6 +32,24 @@ enum class IFPParameter {
GenerationTime,
};
struct CollisionTrackConfig {
bool mcpl_write {false}; //!< Write collision tracks using MCPL?
std::unordered_set<int>
cell_ids; //!< Cell ids where collisions will be written
std::unordered_set<int>
mt_numbers; //!< MT Numbers where collisions will be written
std::unordered_set<int>
universe_ids; //!< Universe IDs where collisions will be written
std::unordered_set<int>
material_ids; //!< Material IDs where collisions will be written
std::unordered_set<std::string>
nuclides; //!< Nuclides where collisions will be written
double deposited_energy_threshold {0.0}; //!< Minimum deposited energy [eV]
int64_t max_collisions {
1000}; //!< Maximum events recorded per collision track file
int64_t max_files {1}; //!< Maximum number of collision track files
};
//==============================================================================
// Global variable declarations
//==============================================================================
@ -41,6 +59,7 @@ namespace settings {
// Boolean flags
extern bool assume_separate; //!< assume tallies are spatially separate?
extern bool check_overlaps; //!< check overlaps in geometry?
extern bool collision_track; //!< flag to use collision track feature?
extern bool confidence_intervals; //!< use confidence intervals for results?
extern bool
create_fission_neutrons; //!< create fission neutrons (fixed source)?
@ -145,11 +164,15 @@ extern std::unordered_set<int>
statepoint_batch; //!< Batches when state should be written
extern std::unordered_set<int>
source_write_surf_id; //!< Surface ids where sources will be written
extern CollisionTrackConfig collision_track_config;
extern double source_rejection_fraction; //!< Minimum fraction of source sites
//!< that must be accepted
extern double free_gas_threshold; //!< Threshold multiplier for free gas
//!< scattering treatment
extern int
max_history_splits; //!< maximum number of particle splits for weight windows
extern int max_secondaries; //!< maximum number of secondaries in the bank
extern int64_t ssw_max_particles; //!< maximum number of particles to be
//!< banked on surfaces per process
extern int64_t ssw_max_files; //!< maximum number of surface source files

View file

@ -22,6 +22,7 @@ constexpr int STATUS_EXIT_ON_TRIGGER {2};
namespace simulation {
extern int ct_current_file; //!< current collision track file index
extern "C" int current_batch; //!< current batch
extern "C" int current_gen; //!< current fission generation
extern "C" bool initialized; //!< has simulation been initialized?

View file

@ -106,6 +106,8 @@ public:
bool writable() const { return writable_; }
bool higher_moments() const { return higher_moments_; }
//----------------------------------------------------------------------------
// Other methods.
@ -190,6 +192,9 @@ private:
//! Whether to multiply by atom density for reaction rates
bool multiply_density_ {true};
//! Whether to accumulate higher moments (third and fourth)
bool higher_moments_ {false};
int64_t index_;
};
@ -203,11 +208,14 @@ extern vector<unique_ptr<Tally>> tallies;
extern vector<int> active_tallies;
extern vector<int> active_analog_tallies;
extern vector<int> active_tracklength_tallies;
extern vector<int> active_timed_tracklength_tallies;
extern vector<int> active_collision_tallies;
extern vector<int> active_meshsurf_tallies;
extern vector<int> active_surface_tallies;
extern vector<int> active_pulse_height_tallies;
extern vector<int> pulse_height_cells;
extern vector<double> time_grid;
} // namespace model
namespace simulation {
@ -239,6 +247,13 @@ void read_tallies_xml(pugi::xml_node root);
//! batch to a new random variable
void accumulate_tallies();
//! Determine distance to next time boundary
//
//! \param time Current time of particle
//! \param speed Speed of particle
//! \return Distance to next time boundary (or INFTY if none)
double distance_to_time_boundary(double time, double speed);
//! Determine which tallies should be active
void setup_active_tallies();

View file

@ -91,6 +91,16 @@ void score_analog_tally_mg(Particle& p);
//! \param distance The distance in [cm] traveled by the particle
void score_tracklength_tally(Particle& p, double distance);
//! Score time filtered tallies using a tracklength estimate of the flux.
//
//! This is triggered at every event (surface crossing, lattice crossing, or
//! collision) and thus cannot be done for tallies that require post-collision
//! information.
//
//! \param p The particle being tracked
//! \param total_distance The distance in [cm] traveled by the particle
void score_timed_tracklength_tally(Particle& p, double total_distance);
//! Score surface or mesh-surface tallies for particle currents.
//
//! \param p The particle being tracked

View file

@ -27,10 +27,10 @@ public:
double heating;
};
Interpolation interp_; //!< interpolation type
int inelastic_flag_; //!< inelastic competition flag
int absorption_flag_; //!< other absorption flag
bool multiply_smooth_; //!< multiply by smooth cross section?
Interpolation interp_; //!< interpolation type
int inelastic_flag_; //!< inelastic competition flag
int absorption_flag_; //!< other absorption flag
bool multiply_smooth_; //!< multiply by smooth cross section?
vector<double> energy_; //!< incident energies
auto n_energy() const { return energy_.size(); }

View file

@ -39,6 +39,9 @@ Use \fIN\fP OpenMP threads.
.B "\-t\fR, \fP\-\-track"
Write tracks for all particles (up to max_tracks).
.TP
.BI \-q " V" "\fR,\fP \-\-verbosity" " V"
Set the output verbosity to \fIV\fP.
.TP
.B "\-v\fR, \fP\-\-version"
Show version information.
.TP

View file

@ -37,7 +37,7 @@ from openmc.tracks import *
from .config import *
# Import a few names from the model module
from openmc.model import Model
from openmc.model import Model, SearchResult
from . import examples

View file

@ -72,6 +72,10 @@ class Cell(IDManagerMixin):
temperature : float or iterable of float
Temperature of the cell in Kelvin. Multiple temperatures can be given
to give each distributed cell instance a unique temperature.
density : float or iterable of float
Density of the cell in [g/cm3]. Multiple densities can be given to give
each distributed cell instance a unique density. Densities set here will
override the density set on materials used to fill the cell.
translation : Iterable of float
If the cell is filled with a universe, this array specifies a vector
that is used to translate (shift) the universe.
@ -109,6 +113,7 @@ class Cell(IDManagerMixin):
self._rotation = None
self._rotation_matrix = None
self._temperature = None
self._density = None
self._translation = None
self._paths = None
self._num_instances = None
@ -146,6 +151,7 @@ class Cell(IDManagerMixin):
if self.fill_type == 'material':
string += '\t{0: <15}=\t{1}\n'.format('Temperature',
self.temperature)
string += '\t{0: <15}=\t{1}\n'.format('Density', self.density)
string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation)
string += '{: <16}=\t{}\n'.format('\tVolume', self.volume)
@ -262,6 +268,30 @@ class Cell(IDManagerMixin):
else:
self._temperature = temperature
@property
def density(self):
return self._density
@density.setter
def density(self, density):
# Make sure densities are greater than zero
cv.check_type('cell density', density, (Iterable, Real), none_ok=True)
if isinstance(density, Iterable):
cv.check_type('cell density', density, Iterable, Real)
for rho in density:
cv.check_greater_than('cell density', rho, 0.0, True)
elif isinstance(density, Real):
cv.check_greater_than('cell density', density, 0.0, True)
# If this cell is filled with a universe or lattice, propagate
# densities to all cells contained. Otherwise, simply assign it.
if self.fill_type in ('universe', 'lattice'):
for c in self.get_all_cells().values():
if c.fill_type == 'material':
c._density = density
else:
self._density = density
@property
def translation(self):
return self._translation
@ -530,6 +560,8 @@ class Cell(IDManagerMixin):
clone.volume = self.volume
if self.temperature is not None:
clone.temperature = self.temperature
if self.density is not None:
clone.density = self.density
if self.translation is not None:
clone.translation = self.translation
if self.rotation is not None:
@ -662,6 +694,12 @@ class Cell(IDManagerMixin):
else:
element.set("temperature", str(self.temperature))
if self.density is not None:
if isinstance(self.density, Iterable):
element.set("density", ' '.join(str(t) for t in self.density))
else:
element.set("density", str(self.density))
if self.translation is not None:
element.set("translation", ' '.join(map(str, self.translation)))
@ -723,10 +761,13 @@ class Cell(IDManagerMixin):
c.temperature = temperature
else:
c.temperature = temperature[0]
density = get_elem_list(elem, 'density', float)
if density is not None:
c.density = density if len(density) > 1 else density[0]
v = get_text(elem, 'volume')
if v is not None:
c.volume = float(v)
for key in ('temperature', 'rotation', 'translation'):
for key in ('temperature', 'density', 'rotation', 'translation'):
values = get_elem_list(elem, key, float)
if values is not None:
if key == 'rotation' and len(values) == 9:

View file

@ -37,7 +37,7 @@ def check_type(name, value, expected_type, expected_iter_type=None, *, none_ok=F
[t.__name__ for t in expected_type]))
else:
msg = (f'Unable to set "{name}" to "{value}" which is not of type "'
f'{expected_type.__name__}"')
f'{expected_type}"')
raise TypeError(msg)
if expected_iter_type:

View file

@ -302,6 +302,8 @@ class DAGMCUniverse(openmc.UniverseBase):
dagmc_element = ET.Element('dagmc_universe')
dagmc_element.set('id', str(self.id))
if self.name:
dagmc_element.set('name', self.name)
if self.auto_geom_ids:
dagmc_element.set('auto_geom_ids', 'true')
if self.auto_mat_ids:

View file

@ -324,7 +324,7 @@ def atomic_mass(isotope):
# isotopes of their element (e.g. C0), calculate the atomic mass as
# the sum of the atomic mass times the natural abundance of the isotopes
# that make up the element.
for element in ['C', 'Zn', 'Pt', 'Os', 'Tl']:
for element in ['C', 'Zn', 'Pt', 'Os', 'Tl', 'V']:
isotope_zero = element.lower() + '0'
_ATOMIC_MASS[isotope_zero] = 0.
for iso, abundance in isotopes(element):

View file

@ -13,7 +13,7 @@ import openmc.checkvalue as cv
from openmc.exceptions import DataError
from openmc.mixin import EqualityMixin
from openmc.stats import Discrete, Tabular, Univariate, combine_distributions
from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
from .data import ATOMIC_NUMBER, gnds_name
from .function import INTERPOLATION_SCHEME
from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
@ -126,9 +126,7 @@ class FissionProductYields(EqualityMixin):
for j in range(n_products):
Z, A = divmod(int(values[4*j]), 1000)
isomeric_state = int(values[4*j + 1])
name = ATOMIC_SYMBOL[Z] + str(A)
if isomeric_state > 0:
name += f'_m{isomeric_state}'
name = gnds_name(Z, A, isomeric_state)
yield_j = ufloat(values[4*j + 2], values[4*j + 3])
yields[name] = yield_j
@ -256,10 +254,7 @@ class DecayMode(EqualityMixin):
A += delta_A
Z += delta_Z
if self._daughter_state > 0:
return f'{ATOMIC_SYMBOL[Z]}{A}_m{self._daughter_state}'
else:
return f'{ATOMIC_SYMBOL[Z]}{A}'
return gnds_name(Z, A, self._daughter_state)
@property
def parent(self):
@ -348,10 +343,7 @@ class Decay(EqualityMixin):
self.nuclide['atomic_number'] = Z
self.nuclide['mass_number'] = A
self.nuclide['isomeric_state'] = metastable
if metastable > 0:
self.nuclide['name'] = f'{ATOMIC_SYMBOL[Z]}{A}_m{metastable}'
else:
self.nuclide['name'] = f'{ATOMIC_SYMBOL[Z]}{A}'
self.nuclide['name'] = gnds_name(Z, A, metastable)
self.nuclide['mass'] = items[1] # AWR
self.nuclide['excited_state'] = items[2] # State of the original nuclide
self.nuclide['stable'] = (items[4] == 1) # Nucleus stability flag
@ -591,7 +583,7 @@ def decay_photon_energy(nuclide: str) -> Univariate | None:
openmc.stats.Univariate or None
Distribution of energies in [eV] of photons emitted from decay, or None
if no photon source exists. Note that the probabilities represent
intensities, given as [Bq].
intensities, given as [Bq/atom] (in other words, decay constants).
"""
if not _DECAY_PHOTON_ENERGY:
chain_file = openmc.config.get('chain_file')

View file

@ -11,7 +11,7 @@ import h5py
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
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)
from .fission_energy import FissionEnergyRelease
@ -678,11 +678,7 @@ class IncidentNeutron(EqualityMixin):
temperature = ev.target['temperature']
# Determine name
element = ATOMIC_SYMBOL[atomic_number]
if metastable > 0:
name = f'{element}{mass_number}_m{metastable}'
else:
name = f'{element}{mass_number}'
name = gnds_name(atomic_number, mass_number, metastable)
# Instantiate incident neutron data
data = cls(name, atomic_number, mass_number, metastable,

View file

@ -17,6 +17,7 @@ from .stepresult import *
from .results import *
from .integrators import *
from .transfer_rates import *
from .r2s import *
from . import abc
from . import cram
from . import helpers

View file

@ -12,6 +12,7 @@ from copy import deepcopy
from inspect import signature
from numbers import Real, Integral
from pathlib import Path
from textwrap import dedent
import time
from typing import Optional, Union, Sequence
from warnings import warn
@ -526,7 +527,7 @@ class Integrator(ABC):
r"""Abstract class for solving the time-integration for depletion
"""
_params = r"""
_params = dedent(r"""
Parameters
----------
operator : openmc.deplete.abc.TransportOperator
@ -617,7 +618,7 @@ class Integrator(ABC):
.. versionadded:: 0.15.3
"""
""")
def __init__(
self,
@ -630,17 +631,7 @@ class Integrator(ABC):
solver: str = "cram48",
continue_timesteps: bool = False,
):
# Check number of stages previously used
if operator.prev_res is not None:
res = operator.prev_res[-1]
if res.data.shape[0] != self._num_stages:
raise ValueError(
"{} incompatible with previous restart calculation. "
"Previous scheme used {} intermediate solutions, while "
"this uses {}".format(
self.__class__.__name__, res.data.shape[0],
self._num_stages))
elif continue_timesteps:
if continue_timesteps and operator.prev_res is None:
raise ValueError("Continuation run requires passing prev_results.")
self.operator = operator
self.chain = operator.chain
@ -774,12 +765,8 @@ class Integrator(ABC):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
n_end : list of numpy.ndarray
Concentrations at end of timestep
"""
@property
@ -810,9 +797,9 @@ class Integrator(ABC):
"""Get beginning of step concentrations, reaction rates from restart"""
res = self.operator.prev_res[-1]
# Depletion methods expect list of arrays
bos_conc = list(res.data[0])
rates = res.rates[0]
k = ufloat(res.k[0, 0], res.k[0, 1])
bos_conc = list(res.data)
rates = res.rates
k = ufloat(res.k[0], res.k[1])
if res.source_rate != 0.0:
# Scale reaction rates by ratio of source rates
@ -854,7 +841,8 @@ class Integrator(ABC):
self,
final_step: bool = True,
output: bool = True,
path: PathLike = 'depletion_results.h5'
path: PathLike = 'depletion_results.h5',
write_rates: bool = False
):
"""Perform the entire depletion process across all steps
@ -873,6 +861,11 @@ class Integrator(ABC):
Path to file to write. Defaults to 'depletion_results.h5'.
.. versionadded:: 0.15.0
write_rates : bool, optional
Whether reaction rates should be written to the results file for
each step. Defaults to ``False`` to reduce file size.
.. versionadded:: 0.15.3
"""
with change_directory(self.operator.output_dir):
n = self.operator.initial_condition()
@ -889,18 +882,22 @@ class Integrator(ABC):
n, res = self._get_bos_data_from_restart(source_rate, n)
# Solve Bateman equations over time interval
proc_time, n_list, res_list = self(n, res.rates, dt, source_rate, i)
proc_time, n_end = self(n, res.rates, dt, source_rate, i)
# Insert BOS concentration, transport results
n_list.insert(0, n)
res_list.insert(0, res)
# Remove actual EOS concentration for next step
n = n_list.pop()
StepResult.save(self.operator, n_list, res_list, [t, t + dt],
source_rate, self._i_res + i, proc_time, path)
StepResult.save(
self.operator,
n,
res,
[t, t + dt],
source_rate,
self._i_res + i,
proc_time,
write_rates=write_rates,
path=path
)
# Update for next step
n = n_end
t += dt
# Final simulation -- in the case that final_step is False, a zero
@ -909,9 +906,18 @@ class Integrator(ABC):
# solve)
if output and final_step and comm.rank == 0:
print(f"[openmc.deplete] t={t} (final operator evaluation)")
res_list = [self.operator(n, source_rate if final_step else 0.0)]
StepResult.save(self.operator, [n], res_list, [t, t],
source_rate, self._i_res + len(self), proc_time, path)
res_final = self.operator(n, source_rate if final_step else 0.0)
StepResult.save(
self.operator,
n,
res_final,
[t, t],
source_rate,
self._i_res + len(self),
proc_time,
write_rates=write_rates,
path=path
)
self.operator.write_bos_data(len(self) + self._i_res)
self.operator.finalize()
@ -1012,6 +1018,36 @@ class Integrator(ABC):
material, composition, rate, rate_units, timesteps)
def add_redox(self, material, buffer, oxidation_states, timesteps=None):
"""Add redox control to depletable material.
Parameters
----------
material : openmc.Material or str or int
Depletable material
buffer : dict
Dictionary of buffer nuclides used to maintain redox balance. Keys
are nuclide names (strings) and values are their respective
fractions (float) that collectively sum to 1.
oxidation_states : dict
User-defined oxidation states for elements. Keys are element symbols
(e.g., 'H', 'He'), and values are their corresponding oxidation
states as integers (e.g., +1, 0).
timesteps : list of int, optional
List of timestep indices where to set external source rates.
Defaults to None, which means the external source rate is set for
all timesteps.
"""
if self.transfer_rates is None:
if hasattr(self.operator, 'model'):
materials = self.operator.model.materials
elif hasattr(self.operator, 'materials'):
materials = self.operator.materials
self.transfer_rates = TransferRates(
self.operator, materials, len(self.timesteps))
self.transfer_rates.set_redox(material, buffer, oxidation_states, timesteps)
@add_params
class SIIntegrator(Integrator):
r"""Abstract class for the Stochastic Implicit Euler integrators
@ -1020,7 +1056,7 @@ class SIIntegrator(Integrator):
the number of particles used in initial transport calculation
"""
_params = r"""
_params = dedent(r"""
Parameters
----------
operator : openmc.deplete.abc.TransportOperator
@ -1108,7 +1144,7 @@ class SIIntegrator(Integrator):
.. versionadded:: 0.12
"""
""")
def __init__(
self,
@ -1140,10 +1176,40 @@ class SIIntegrator(Integrator):
self.operator.settings.particles //= self.n_steps
return inherited
@abstractmethod
def __call__(self, n, rates, dt, source_rate, i):
"""Perform the integration across one time step
Parameters
----------
n : list of numpy.ndarray
List of atom number arrays for each material. Each array has
shape ``(n_nucs,)`` where ``n_nucs`` is the number of nuclides
rates : openmc.deplete.ReactionRates
Reaction rates (from transport operator)
dt : float
Time step in [s]
source_rate : float
Power in [W] or source rate in [neutron/sec]
i : int
Current time step index
Returns
-------
proc_time : float
Time spent in transport simulation
n_end : list of numpy.ndarray
Updated atom number densities for each material
op_result : OperatorResult
Eigenvalue and reaction rates resulting from transport simulation
"""
def integrate(
self,
output: bool = True,
path: PathLike = "depletion_results.h5"
path: PathLike = "depletion_results.h5",
write_rates: bool = False
):
"""Perform the entire depletion process across all steps
@ -1155,11 +1221,17 @@ class SIIntegrator(Integrator):
Path to file to write. Defaults to 'depletion_results.h5'.
.. versionadded:: 0.15.0
write_rates : bool, optional
Whether reaction rates should be written to the results file for
each step. Defaults to ``False`` to reduce file size.
.. versionadded:: 0.15.3
"""
with change_directory(self.operator.output_dir):
n = self.operator.initial_condition()
t, self._i_res = self._get_start_data()
res_end = None # Will be set in first iteration
for i, (dt, p) in enumerate(self):
if output:
print(f"[openmc.deplete] t={t} s, dt={dt} s, source={p}")
@ -1169,28 +1241,38 @@ class SIIntegrator(Integrator):
n, res = self._get_bos_data_from_operator(i, p, n)
else:
n, res = self._get_bos_data_from_restart(p, n)
else:
# Pull rates, k from previous iteration w/o
# re-running transport
res = res_list[-1] # defined in previous i iteration
proc_time, n_list, res_list = self(n, res.rates, dt, p, i)
proc_time, n_end, res_end = self(n, res.rates, dt, p, i)
# Insert BOS concentration, transport results
n_list.insert(0, n)
res_list.insert(0, res)
# Remove actual EOS concentration for next step
n = n_list.pop()
StepResult.save(self.operator, n_list, res_list, [t, t + dt],
p, self._i_res + i, proc_time, path)
StepResult.save(
self.operator,
n,
res,
[t, t + dt],
p,
self._i_res + i,
proc_time,
write_rates=write_rates,
path=path
)
# Update for next step
n = n_end
res = res_end
t += dt
# No final simulation for SIE, use last iteration results
StepResult.save(self.operator, [n], [res_list[-1]], [t, t],
p, self._i_res + len(self), proc_time, path)
StepResult.save(
self.operator,
n,
res_end,
[t, t],
p,
self._i_res + len(self),
proc_time,
write_rates=write_rates,
path=path
)
self.operator.write_bos_data(self._i_res + len(self))
self.operator.finalize()

View file

@ -7,6 +7,7 @@ loaded from an .xml file and all the nuclides are linked together.
from io import StringIO
from itertools import chain
import math
import numpy as np
import re
from collections import defaultdict, namedtuple
from collections.abc import Mapping, Iterable
@ -627,9 +628,15 @@ class Chain:
"""
reactions = set()
# Use DOK matrix as intermediate representation for matrix
n = len(self)
matrix = sp.dok_matrix((n, n))
# we accumulate indices and value entries for everything and create the matrix
# in one step at the end to avoid expensive index checks scipy otherwise does.
rows, cols, vals = [], [], []
def setval(i, j, val):
rows.append(i)
cols.append(j)
vals.append(val)
if fission_yields is None:
fission_yields = self.get_default_fission_yields()
@ -639,7 +646,7 @@ class Chain:
if nuc.half_life is not None:
decay_constant = math.log(2) / nuc.half_life
if decay_constant != 0.0:
matrix[i, i] -= decay_constant
setval(i, i, -decay_constant)
# Gain from radioactive decay
if nuc.n_decay_modes != 0:
@ -650,19 +657,19 @@ class Chain:
if branch_val != 0.0:
if target is not None:
k = self.nuclide_dict[target]
matrix[k, i] += branch_val
setval(k, i, branch_val)
# Produce alphas and protons from decay
if 'alpha' in decay_type:
k = self.nuclide_dict.get('He4')
if k is not None:
count = decay_type.count('alpha')
matrix[k, i] += count * branch_val
setval(k, i, count * branch_val)
elif 'p' in decay_type:
k = self.nuclide_dict.get('H1')
if k is not None:
count = decay_type.count('p')
matrix[k, i] += count * branch_val
setval(k, i, count * branch_val)
if nuc.name in rates.index_nuc:
# Extract all reactions for this nuclide in this cell
@ -679,13 +686,13 @@ class Chain:
if r_type not in reactions:
reactions.add(r_type)
if path_rate != 0.0:
matrix[i, i] -= path_rate
setval(i, i, -path_rate)
# Gain term; allow for total annihilation for debug purposes
if r_type != 'fission':
if target is not None and path_rate != 0.0:
k = self.nuclide_dict[target]
matrix[k, i] += path_rate * br
setval(k, i, path_rate * br)
# Determine light nuclide production, e.g., (n,d) should
# produce H2
@ -693,20 +700,76 @@ class Chain:
for light_nuc in light_nucs:
k = self.nuclide_dict.get(light_nuc)
if k is not None:
matrix[k, i] += path_rate * br
setval(k, i, path_rate * br)
else:
for product, y in fission_yields[nuc.name].items():
yield_val = y * path_rate
if yield_val != 0.0:
k = self.nuclide_dict[product]
matrix[k, i] += yield_val
setval(k, i, yield_val)
# Clear set of reactions
reactions.clear()
# Return CSC representation instead of DOK
return matrix.tocsc()
return sp.csc_matrix((vals, (rows, cols)), shape=(n, n))
def add_redox_term(self, matrix, buffer, oxidation_states):
r"""Adds a redox term to the depletion matrix from data contained in
the matrix itself and a few user-inputs.
The redox term to add to the buffer nuclide :math:`N_j` can be written
as:
.. math::
\frac{dN_j(t)}{dt} = \cdots - \frac{1}{OS_j}\sum_i N_i a_{ij}
\cdot OS_i
where :math:`OS` is the oxidation states vector and :math:`a_{ij}` the
corresponding term in the Bateman matrix.
Parameters
----------
matrix : scipy.sparse.csc_matrix
Sparse matrix representing depletion
buffer : dict
Dictionary of buffer nuclides used to maintain anoins net balance.
Keys are nuclide names (strings) and values are their respective
fractions (float) that collectively sum to 1.
oxidation_states : dict
User-defined oxidation states for elements. Keys are element symbols
(e.g., 'H', 'He'), and values are their corresponding oxidation
states as integers (e.g., +1, 0).
Returns
-------
matrix : scipy.sparse.csc_matrix
Sparse matrix with redox term added
"""
# Elements list with the same size as self.nuclides
elements = [re.split(r'\d+', nuc.name)[0] for nuc in self.nuclides]
# Match oxidation states with all elements and add 0 if not data
os = np.array([oxidation_states[elm] if elm in oxidation_states else 0
for elm in elements])
# Buffer idx with nuclide index as value
buffer_idx = {nuc: self.nuclide_dict[nuc] for nuc in buffer}
array = matrix.toarray()
redox_change = np.array([])
# calculate the redox array
for i in range(len(self)):
# Net redox impact of reaction: multiply the i-th column of the
# depletion matrix by the oxidation states
redox_change = np.append(redox_change, sum(array[:, i]*os))
# Subtract redox vector to the buffer nuclides in the matrix scaling by
# their respective oxidation states
for nuc, idx in buffer_idx.items():
array[idx] -= redox_change * buffer[nuc] / os[idx]
return sp.csc_matrix(array)
def form_rr_term(self, tr_rates, current_timestep, mats):
"""Function to form the transfer rate term matrices.

View file

@ -5,7 +5,7 @@ shutdown dose rate calculations.
"""
from copy import deepcopy
from copy import copy
from typing import Sequence
from math import log, prod
@ -108,14 +108,15 @@ def time_correction_factors(
# Create a 2D array for the time correction factors
h = np.zeros((n_timesteps, n_nuclides))
for i, (dt, rate) in enumerate(zip(timesteps, source_rates)):
# Precompute the exponential terms. Since (1 - exp(-x)) is susceptible to
# roundoff error, use expm1 instead (which computes exp(x) - 1)
g = np.exp(-decay_rate*dt)
one_minus_g = -np.expm1(-decay_rate*dt)
# Precompute all exponential terms with same shape as h
decay_dt = decay_rate[np.newaxis, :] * timesteps[:, np.newaxis]
g = np.exp(-decay_dt)
one_minus_g = -np.expm1(-decay_dt)
# Apply recurrence relation step by step
for i in range(len(timesteps)):
# Eq. (4) in doi:10.1016/j.fusengdes.2019.111399
h[i + 1] = rate*one_minus_g + h[i]*g
h[i + 1] = source_rates[i] * one_minus_g[i] + h[i] * g[i]
return {nuclides[i]: h[:, i] for i in range(n_nuclides)}
@ -163,8 +164,12 @@ def apply_time_correction(
radionuclides = [str(x) for x in tally.filters[i_filter].bins]
tcf = np.array([time_correction_factors[x][index] for x in radionuclides])
# Create copy of tally
new_tally = deepcopy(tally)
# Force tally results to be read and std_dev to be computed
tally.std_dev
# Create shallow copy of tally
new_tally = copy(tally)
new_tally._filters = copy(tally._filters)
# Determine number of bins in other filters
n_bins_before = prod([f.num_bins for f in tally.filters[:i_filter]])
@ -176,32 +181,33 @@ def apply_time_correction(
shape = (n_bins_before, n_radionuclides, n_bins_after, n_nuclides, n_scores)
tally_sum = new_tally.sum.reshape(shape)
tally_sum_sq = new_tally.sum_sq.reshape(shape)
tally_mean = new_tally.mean.reshape(shape)
tally_std_dev = new_tally.std_dev.reshape(shape)
# 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 = None
new_tally._std_dev = None
new_tally._mean = tally_mean * tcf
new_tally._std_dev = tally_std_dev * tcf
shape = (-1, n_nuclides, n_scores)
if sum_nuclides:
# Query the mean and standard deviation
mean = new_tally.mean
std_dev = new_tally.std_dev
# Sum over parent nuclides (note that when combining different bins for
# parent nuclide, we can't work directly on sum_sq)
new_tally._mean = mean.sum(axis=1).reshape(shape)
new_tally._std_dev = np.linalg.norm(std_dev, axis=1).reshape(shape)
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
# 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
new_tally._mean.shape = shape
new_tally._std_dev.shape = shape
return new_tally

View file

@ -22,7 +22,6 @@ class PredictorIntegrator(Integrator):
.. math::
\mathbf{n}_{i+1} = \exp\left(h\mathbf{A}(\mathbf{n}_i) \right) \mathbf{n}_i
"""
_num_stages = 1
@ -47,15 +46,12 @@ class PredictorIntegrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
n_end : list of numpy.ndarray
Concentrations at end of interval
op_results : empty list
Kept for consistency with API. No intermediate calls to operator
with predictor
"""
proc_time, n_end = self._timed_deplete(n, rates, dt, _i)
return proc_time, [n_end], []
return proc_time, n_end
@add_params
@ -99,11 +95,8 @@ class CECMIntegrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from transport simulations
n_end : list of numpy.ndarray
Concentrations at end of interval
"""
# deplete across first half of interval
time0, n_middle = self._timed_deplete(n, rates, dt / 2, _i)
@ -113,7 +106,7 @@ class CECMIntegrator(Integrator):
# MOS reaction rates
time1, n_end = self._timed_deplete(n, res_middle.rates, dt, _i)
return time0 + time1, [n_middle, n_end], [res_middle]
return time0 + time1, n_end
@add_params
@ -163,12 +156,8 @@ class CF4Integrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
n_end : list of numpy.ndarray
Concentrations at end of interval
"""
# Step 1: deplete with matrix 1/2*A(y0)
time1, n_eos1 = self._timed_deplete(
@ -193,9 +182,7 @@ class CF4Integrator(Integrator):
time5, n_eos5 = self._timed_deplete(
n_inter, list_rates, dt, _i, matrix_func=cf4_f4)
return (time1 + time2 + time3 + time4 + time5,
[n_eos1, n_eos2, n_eos3, n_eos5],
[res1, res2, res3])
return time1 + time2 + time3 + time4 + time5, n_eos5
@add_params
@ -241,12 +228,8 @@ class CELIIntegrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulation
n_end : list of numpy.ndarray
Concentrations at end of interval
"""
# deplete to end using BOS rates
proc_time, n_ce = self._timed_deplete(n_bos, rates, dt, _i)
@ -261,7 +244,7 @@ class CELIIntegrator(Integrator):
time_le2, n_end = self._timed_deplete(
n_inter, list_rates, dt, _i, matrix_func=celi_f2)
return proc_time + time_le1 + time_le1, [n_ce, n_end], [res_ce]
return proc_time + time_le1 + time_le2, n_end
@add_params
@ -307,12 +290,8 @@ class EPCRK4Integrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
n_end : list of numpy.ndarray
Concentrations at end of interval
"""
# Step 1: deplete with matrix A(y0) / 2
@ -331,7 +310,7 @@ class EPCRK4Integrator(Integrator):
list_rates = list(zip(rates, res1.rates, res2.rates, res3.rates))
time4, n4 = self._timed_deplete(n, list_rates, dt, _i, matrix_func=rk4_f4)
return (time1 + time2 + time3 + time4, [n1, n2, n3, n4], [res1, res2, res3])
return time1 + time2 + time3 + time4, n4
@add_params
@ -389,12 +368,8 @@ class LEQIIntegrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulation
n_end : list of numpy.ndarray
Concentrations at end of interval
"""
if i == 0:
if self._i_res < 1: # need at least previous transport solution
@ -403,7 +378,7 @@ class LEQIIntegrator(Integrator):
self, n_bos, bos_rates, dt, source_rate, i)
prev_res = self.operator.prev_res[-2]
prev_dt = self.timesteps[i] - prev_res.time[0]
self._prev_rates = prev_res.rates[0]
self._prev_rates = prev_res.rates
else:
prev_dt = self.timesteps[i - 1]
@ -432,9 +407,7 @@ class LEQIIntegrator(Integrator):
# store updated rates
self._prev_rates = copy.deepcopy(bos_res.rates)
return (
time1 + time2 + time3 + time4, [n_eos0, n_eos1],
[bos_res, res_inter])
return time1 + time2 + time3 + time4, n_eos1
@add_params
@ -471,10 +444,9 @@ class SICELIIntegrator(SIIntegrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_bos_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
n_end : list of numpy.ndarray
Concentrations at end of interval
op_result : openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
"""
@ -500,7 +472,7 @@ class SICELIIntegrator(SIIntegrator):
proc_time += time1 + time2
# end iteration
return proc_time, [n_eos, n_inter], [res_bar]
return proc_time, n_inter, res_bar
@add_params
@ -537,10 +509,9 @@ class SILEQIIntegrator(SIIntegrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
n_end : list of numpy.ndarray
Concentrations at end of interval
op_result : openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulation
"""
@ -552,7 +523,7 @@ class SILEQIIntegrator(SIIntegrator):
self, n_bos, bos_rates, dt, source_rate, i)
prev_res = self.operator.prev_res[-2]
prev_dt = self.timesteps[i] - prev_res.time[0]
self._prev_rates = prev_res.rates[0]
self._prev_rates = prev_res.rates
else:
prev_dt = self.timesteps[i - 1]
@ -585,7 +556,10 @@ class SILEQIIntegrator(SIIntegrator):
n_inter, inputs, dt, i, matrix_func=leqi_f4)
proc_time += time1 + time2
return proc_time, [n_eos, n_inter], [res_bar]
# Store updated rates for next step
self._prev_rates = copy.deepcopy(bos_rates)
return proc_time, n_inter, res_bar
integrator_by_name = {

View file

@ -8,8 +8,9 @@ from __future__ import annotations
from collections.abc import Sequence
import shutil
from tempfile import TemporaryDirectory
from typing import Union, TypeAlias
from typing import Union, TypeAlias, Self
import h5py
import pandas as pd
import numpy as np
@ -20,6 +21,7 @@ from openmc.data import REACTION_MT
import openmc
from .chain import Chain, REACTIONS, _get_chain
from .coupled_operator import _find_cross_sections, _get_nuclides_with_data
from ..utility_funcs import h5py_file_or_group
import openmc.lib
from openmc.mpi import comm
@ -47,6 +49,7 @@ def get_microxs_and_flux(
reaction_rate_mode: str = 'direct',
chain_file: PathLike | Chain | None = None,
path_statepoint: PathLike | None = None,
path_input: PathLike | None = None,
run_kwargs=None
) -> tuple[list[np.ndarray], list[MicroXS]]:
"""Generate microscopic cross sections and fluxes for multiple domains.
@ -59,7 +62,7 @@ def get_microxs_and_flux(
.. versionadded:: 0.14.0
.. versionchanged:: 0.15.3
Added `reaction_rate_mode` and `path_statepoint` arguments.
Added `reaction_rate_mode`, `path_statepoint`, `path_input` arguments.
Parameters
----------
@ -90,6 +93,10 @@ def get_microxs_and_flux(
Path to write the statepoint file from the neutron transport solve to.
By default, The statepoint file is written to a temporary directory and
is not kept.
path_input : path-like, optional
Path to write the model XML file from the neutron transport solve to.
By default, the model XML file is written to a temporary directory and
not kept.
run_kwargs : dict, optional
Keyword arguments passed to :meth:`openmc.Model.run`
@ -108,7 +115,7 @@ def get_microxs_and_flux(
check_value('reaction_rate_mode', reaction_rate_mode, {'direct', 'flux'})
# Save any original tallies on the model
original_tallies = model.tallies
original_tallies = list(model.tallies)
# Determine what reactions and nuclides are available in chain
chain = _get_chain(chain_file)
@ -163,14 +170,16 @@ def get_microxs_and_flux(
# Reinitialize with tallies
openmc.lib.init(intracomm=comm)
# create temporary run
with TemporaryDirectory() as temp_dir:
if run_kwargs is None:
run_kwargs = {}
else:
run_kwargs = dict(run_kwargs)
run_kwargs.setdefault('cwd', temp_dir)
# Indicate to run in temporary directory unless being executed through
# openmc.lib, in which case we don't need to specify the cwd
run_kwargs = dict(run_kwargs) if run_kwargs else {}
if not openmc.lib.is_initialized:
run_kwargs.setdefault('cwd', temp_dir)
# Run transport simulation and synchronize
statepoint_path = model.run(**run_kwargs)
comm.barrier()
if comm.rank == 0:
# Move the statepoint file if it is being saved to a specific path
@ -178,15 +187,22 @@ def get_microxs_and_flux(
shutil.move(statepoint_path, path_statepoint)
statepoint_path = path_statepoint
with StatePoint(statepoint_path) as sp:
if reaction_rate_mode == 'direct':
rr_tally = sp.tallies[rr_tally.id]
rr_tally._read_results()
flux_tally = sp.tallies[flux_tally.id]
flux_tally._read_results()
# Export the model to path_input if provided
if path_input is not None:
model.export_to_model_xml(path_input)
# Broadcast updated statepoint path to all ranks
statepoint_path = comm.bcast(statepoint_path)
# Read in tally results (on all ranks)
with StatePoint(statepoint_path) as sp:
if reaction_rate_mode == 'direct':
rr_tally = sp.tallies[rr_tally.id]
rr_tally._read_results()
flux_tally = sp.tallies[flux_tally.id]
flux_tally._read_results()
# Get flux values and make energy groups last dimension
flux_tally = comm.bcast(flux_tally)
flux = flux_tally.get_reshaped_data() # (domains, groups, 1, 1)
flux = np.moveaxis(flux, 1, -1) # (domains, 1, 1, groups)
@ -195,7 +211,6 @@ def get_microxs_and_flux(
if reaction_rate_mode == 'direct':
# Get reaction rates
rr_tally = comm.bcast(rr_tally)
reaction_rates = rr_tally.get_reshaped_data() # (domains, groups, nuclides, reactions)
# Make energy groups last dimension
@ -345,13 +360,17 @@ class MicroXS:
reactions = chain.reactions
mts = [REACTION_MT[name] for name in reactions]
# Normalize multigroup flux
multigroup_flux = np.array(multigroup_flux)
multigroup_flux /= multigroup_flux.sum()
# Create 3D array for microscopic cross sections
microxs_arr = np.zeros((len(nuclides), len(mts), 1))
# If flux is zero, safely return zero cross sections
multigroup_flux = np.array(multigroup_flux)
if (flux_sum := multigroup_flux.sum()) == 0.0:
return cls(microxs_arr, nuclides, reactions)
# Normalize multigroup flux
multigroup_flux /= flux_sum
# Compute microscopic cross sections within a temporary session
with openmc.lib.TemporarySession(**init_kwargs):
# For each nuclide and reaction, compute the flux-averaged xs
@ -383,8 +402,7 @@ class MicroXS:
MicroXS
"""
if 'float_precision' not in kwargs:
kwargs['float_precision'] = 'round_trip'
kwargs.setdefault('float_precision', 'round_trip')
df = pd.read_csv(csv_file, **kwargs)
df.set_index(['nuclides', 'reactions', 'groups'], inplace=True)
@ -419,3 +437,96 @@ class MicroXS:
)
df = pd.DataFrame({'xs': self.data.flatten()}, index=multi_index)
df.to_csv(*args, **kwargs)
def to_hdf5(self, group_or_filename: h5py.Group | PathLike, **kwargs):
"""Export microscopic cross section data to HDF5 format
Parameters
----------
group_or_filename : h5py.Group or path-like
HDF5 group or filename to write to
kwargs : dict, optional
Keyword arguments to pass to :meth:`h5py.Group.create_dataset`.
Defaults to {'compression': 'lzf'}.
"""
kwargs.setdefault('compression', 'lzf')
with h5py_file_or_group(group_or_filename, 'w') as group:
# Store cross section data as 3D dataset
group.create_dataset('data', data=self.data, **kwargs)
# Store metadata as datasets using string encoding
group.create_dataset('nuclides', data=np.array(self.nuclides, dtype='S'))
group.create_dataset('reactions', data=np.array(self.reactions, dtype='S'))
@classmethod
def from_hdf5(cls, group_or_filename: h5py.Group | PathLike) -> Self:
"""Load data from an HDF5 file
Parameters
----------
group_or_filename : h5py.Group or str or PathLike
HDF5 group or path to HDF5 file. If given as an h5py.Group, the
data is read from that group. If given as a string, it is assumed
to be the filename for the HDF5 file.
Returns
-------
MicroXS
"""
with h5py_file_or_group(group_or_filename, 'r') as group:
# Read data from HDF5 group
data = group['data'][:]
nuclides = [nuc.decode('utf-8') for nuc in group['nuclides'][:]]
reactions = [rxn.decode('utf-8') for rxn in group['reactions'][:]]
return cls(data, nuclides, reactions)
def write_microxs_hdf5(
micros: Sequence[MicroXS],
filename: PathLike,
names: Sequence[str] | None = None,
**kwargs
):
"""Write multiple MicroXS objects to an HDF5 file
Parameters
----------
micros : list of MicroXS
List of MicroXS objects
filename : PathLike
Output HDF5 filename
names : list of str, optional
Names for each MicroXS object. If None, uses 'domain_0', 'domain_1',
etc.
**kwargs
Additional keyword arguments passed to :meth:`h5py.Group.create_dataset`
"""
if names is None:
names = [f'domain_{i}' for i in range(len(micros))]
# Open file once and write all domains using group interface
with h5py.File(filename, 'w') as f:
for microxs, name in zip(micros, names):
group = f.create_group(name)
microxs.to_hdf5(group, **kwargs)
def read_microxs_hdf5(filename: PathLike) -> dict[str, MicroXS]:
"""Read multiple MicroXS objects from an HDF5 file
Parameters
----------
filename : path-like
HDF5 filename
Returns
-------
dict
Dictionary mapping domain names to MicroXS objects
"""
with h5py.File(filename, 'r') as f:
return {name: MicroXS.from_hdf5(group) for name, group in f.items()}

View file

@ -109,6 +109,13 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
matrices = [matrix - transfer for (matrix, transfer) in zip(matrices,
transfers)]
if transfer_rates.redox:
for mat_idx, mat_id in enumerate(transfer_rates.local_mats):
if mat_id in transfer_rates.redox:
matrices[mat_idx] = chain.add_redox_term(matrices[mat_idx],
transfer_rates.redox[mat_id][0],
transfer_rates.redox[mat_id][1])
if current_timestep in transfer_rates.index_transfer:
# Gather all on comm.rank 0
matrices = comm.gather(matrices)
@ -125,6 +132,12 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
transfer_matrix = chain.form_rr_term(transfer_rates,
current_timestep,
mat_pair)
# check if destination material has a redox control
if mat_pair[0] in transfer_rates.redox:
transfer_matrix = chain.add_redox_term(transfer_matrix,
transfer_rates.redox[mat_pair[0]][0],
transfer_rates.redox[mat_pair[0]][1])
transfer_pair[mat_pair] = transfer_matrix
# Combine all matrices together in a single matrix of matrices

693
openmc/deplete/r2s.py Normal file
View file

@ -0,0 +1,693 @@
from __future__ import annotations
from collections.abc import Sequence
import copy
from datetime import datetime
import json
from pathlib import Path
import numpy as np
import openmc
from . import IndependentOperator, PredictorIntegrator
from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5
from .results import Results
from ..checkvalue import PathLike
from ..mpi import comm
from openmc.lib import TemporarySession
from openmc.utility_funcs import change_directory
def get_activation_materials(
model: openmc.Model, mmv: openmc.MeshMaterialVolumes
) -> openmc.Materials:
"""Get a list of activation materials for each mesh element/material.
When performing a mesh-based R2S calculation, a unique material is needed
for each activation region, which is a combination of a mesh element and a
material within that mesh element. This function generates a list of such
materials, each with a unique name and volume corresponding to the mesh
element and material.
Parameters
----------
model : openmc.Model
The full model containing the geometry and materials.
mmv : openmc.MeshMaterialVolumes
The mesh material volumes object containing the materials and their
volumes for each mesh element.
Returns
-------
openmc.Materials
A list of materials, each corresponding to a unique mesh element and
material combination.
"""
# Get the material ID, volume, and element index for each element-material
# combination
mat_ids = mmv._materials[mmv._materials > -1]
volumes = mmv._volumes[mmv._materials > -1]
elems, _ = np.where(mmv._materials > -1)
# Get all materials in the model
material_dict = model._get_all_materials()
# Create a new activation material for each element-material combination
materials = openmc.Materials()
for elem, mat_id, vol in zip(elems, mat_ids, volumes):
mat = material_dict[mat_id]
new_mat = mat.clone()
new_mat.depletable = True
new_mat.name = f'Element {elem}, Material {mat_id}'
new_mat.volume = vol
materials.append(new_mat)
return materials
class R2SManager:
"""Manager for Rigorous 2-Step (R2S) method calculations.
This class is responsible for managing the materials and sources needed for
mesh-based or cell-based R2S calculations. It provides methods to get
activation materials and decay photon sources based on the mesh/cells and
materials in the OpenMC model.
This class supports the use of a different models for the neutron and photon
transport calculation. However, for cell-based calculations, it assumes that
the only changes in the model are material assignments. For mesh-based
calculations, it checks material assignments in the photon model and any
element--material combinations that don't appear in the photon model are
skipped.
Parameters
----------
neutron_model : openmc.Model
The OpenMC model to use for neutron transport.
domains : openmc.MeshBase or Sequence[openmc.Cell]
The mesh or a sequence of cells that represent the spatial units over
which the R2S calculation will be performed.
photon_model : openmc.Model, optional
The OpenMC model to use for photon transport calculations. If None, a
shallow copy of the neutron_model will be created and used.
Attributes
----------
domains : openmc.MeshBase or Sequence[openmc.Cell]
The mesh or a sequence of cells that represent the spatial units over
which the R2S calculation will be performed.
neutron_model : openmc.Model
The OpenMC model used for neutron transport.
photon_model : openmc.Model
The OpenMC model used for photon transport calculations.
method : {'mesh-based', 'cell-based'}
Indicates whether the R2S calculation uses mesh elements ('mesh-based')
as the spatial discetization or a list of a cells ('cell-based').
results : dict
A dictionary that stores results from the R2S calculation.
"""
def __init__(
self,
neutron_model: openmc.Model,
domains: openmc.MeshBase | Sequence[openmc.Cell],
photon_model: openmc.Model | None = None,
):
self.neutron_model = neutron_model
if photon_model is None:
# Create a shallow copy of the neutron model for photon transport
self.photon_model = openmc.Model(
geometry=copy.copy(neutron_model.geometry),
materials=copy.copy(neutron_model.materials),
settings=copy.copy(neutron_model.settings),
tallies=copy.copy(neutron_model.tallies),
plots=copy.copy(neutron_model.plots),
)
else:
self.photon_model = photon_model
if isinstance(domains, openmc.MeshBase):
self.method = 'mesh-based'
else:
self.method = 'cell-based'
self.domains = domains
self.results = {}
def run(
self,
timesteps: Sequence[float] | Sequence[tuple[float, str]],
source_rates: float | Sequence[float],
timestep_units: str = 's',
photon_time_indices: Sequence[int] | None = None,
output_dir: PathLike | None = None,
bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
chain_file: PathLike | None = None,
micro_kwargs: dict | None = None,
mat_vol_kwargs: dict | None = None,
run_kwargs: dict | None = None,
operator_kwargs: dict | None = None,
):
"""Run the R2S calculation.
Parameters
----------
timesteps : Sequence[float] or Sequence[tuple[float, str]]
Sequence of timesteps. Note that values are not cumulative. The
units are specified by the `timestep_units` argument when
`timesteps` is an iterable of float. Alternatively, units can be
specified for each step by passing an iterable of (value, unit)
tuples.
source_rates : float or Sequence[float]
Source rate in [neutron/sec] for each interval in `timesteps`.
timestep_units : {'s', 'min', 'h', 'd', 'a'}, optional
Units for values specified in the `timesteps` argument when passing
float values. 's' means seconds, 'min' means minutes, 'h' means
hours, 'd' means days, and 'a' means years (Julian).
photon_time_indices : Sequence[int], optional
Sequence of time indices at which photon transport should be run;
represented as indices into the array of times formed by the
timesteps. For example, if two timesteps are specified, the array of
times would contain three entries, and [2] would indicate computing
photon results at the last time. A value of None indicates to run
photon transport for each time.
output_dir : PathLike, optional
Path to directory where R2S calculation outputs will be saved. If
not provided, a timestamped directory 'r2s_YYYY-MM-DDTHH-MM-SS' is
created. Subdirectories will be created for the neutron transport,
activation, and photon transport steps.
bounding_boxes : dict[int, openmc.BoundingBox], optional
Dictionary mapping cell IDs to bounding boxes used for spatial
source sampling in cell-based R2S calculations. Required if method
is 'cell-based'.
chain_file : PathLike, optional
Path to the depletion chain XML file to use during activation. If
not provided, the default configured chain file will be used.
micro_kwargs : dict, optional
Additional keyword arguments passed to
:func:`openmc.deplete.get_microxs_and_flux` during the neutron
transport step.
mat_vol_kwargs : dict, optional
Additional keyword arguments passed to
:meth:`openmc.MeshBase.material_volumes`.
run_kwargs : dict, optional
Additional keyword arguments passed to :meth:`openmc.Model.run`
during the neutron and photon transport step. By default, output is
disabled.
operator_kwargs : dict, optional
Additional keyword arguments passed to
:class:`openmc.deplete.IndependentOperator`.
Returns
-------
Path
Path to the output directory containing all calculation results
"""
if output_dir is None:
# Create timestamped output directory and broadcast to all ranks for
# consistency (different ranks may have slightly different times)
stamp = datetime.now().strftime('%Y-%m-%dT%H-%M-%S')
output_dir = Path(comm.bcast(f'r2s_{stamp}'))
# Set run_kwargs for the neutron transport step
if micro_kwargs is None:
micro_kwargs = {}
if run_kwargs is None:
run_kwargs = {}
if operator_kwargs is None:
operator_kwargs = {}
run_kwargs.setdefault('output', False)
micro_kwargs.setdefault('run_kwargs', run_kwargs)
# If a chain file is provided, prefer it for steps 1 and 2
if chain_file is not None:
micro_kwargs.setdefault('chain_file', chain_file)
operator_kwargs.setdefault('chain_file', chain_file)
self.step1_neutron_transport(
output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs
)
self.step2_activation(
timesteps, source_rates, timestep_units, output_dir / 'activation',
operator_kwargs=operator_kwargs
)
self.step3_photon_transport(
photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
)
return output_dir
def step1_neutron_transport(
self,
output_dir: PathLike = "neutron_transport",
mat_vol_kwargs: dict | None = None,
micro_kwargs: dict | None = None
):
"""Run the neutron transport step.
This step computes the material volume fractions on the mesh, creates a
mesh-material filter, and retrieves the fluxes and microscopic cross
sections for each mesh/material combination. This step will populate the
'fluxes' and 'micros' keys in the results dictionary. For a mesh-based
calculation, it will also populate the 'mesh_material_volumes' key.
Parameters
----------
output_dir : PathLike, optional
The directory where the results will be saved.
mat_vol_kwargs : dict, optional
Additional keyword arguments based to
:meth:`openmc.MeshBase.material_volumes`.
micro_kwargs : dict, optional
Additional keyword arguments passed to
:func:`openmc.deplete.get_microxs_and_flux`.
"""
output_dir = Path(output_dir).resolve()
output_dir.mkdir(parents=True, exist_ok=True)
if self.method == 'mesh-based':
# Compute material volume fractions on the mesh
if mat_vol_kwargs is None:
mat_vol_kwargs = {}
self.results['mesh_material_volumes'] = mmv = comm.bcast(
self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs))
# Save results to file
if comm.rank == 0:
mmv.save(output_dir / 'mesh_material_volumes.npz')
# Create mesh-material filter based on what combos were found
domains = openmc.MeshMaterialFilter.from_volumes(self.domains, mmv)
else:
domains: Sequence[openmc.Cell] = self.domains
# Check to make sure that each cell is filled with a material and
# that the volume has been set
# TODO: If volumes are not set, run volume calculation for cells
for cell in domains:
if cell.fill is None:
raise ValueError(
f"Cell {cell.id} is not filled with a materials. "
"Please set the fill material for each cell before "
"running the R2S calculation."
)
if cell.volume is None:
raise ValueError(
f"Cell {cell.id} does not have a volume set. "
"Please set the volume for each cell before running "
"the R2S calculation."
)
# Set default keyword arguments for microxs and flux calculation
if micro_kwargs is None:
micro_kwargs = {}
micro_kwargs.setdefault('path_statepoint', output_dir / 'statepoint.h5')
micro_kwargs.setdefault('path_input', output_dir / 'model.xml')
# Run neutron transport and get fluxes and micros. Run via openmc.lib to
# maintain a consistent parallelism strategy with the activation step.
with TemporarySession():
self.results['fluxes'], self.results['micros'] = get_microxs_and_flux(
self.neutron_model, domains, **micro_kwargs)
# Save flux and micros to file
if comm.rank == 0:
np.save(output_dir / 'fluxes.npy', self.results['fluxes'])
write_microxs_hdf5(self.results['micros'], output_dir / 'micros.h5')
def step2_activation(
self,
timesteps: Sequence[float] | Sequence[tuple[float, str]],
source_rates: float | Sequence[float],
timestep_units: str = 's',
output_dir: PathLike = 'activation',
operator_kwargs: dict | None = None,
):
"""Run the activation step.
This step creates a unique copy of each activation material based on the
mesh elements or cells, then solves the depletion equations for each
material using the fluxes and microscopic cross sections obtained in the
neutron transport step. This step will populate the 'depletion_results'
and 'activation_materials' keys in the results dictionary.
Parameters
----------
timesteps : Sequence[float] or Sequence[tuple[float, str]]
Sequence of timesteps. Note that values are not cumulative. The
units are specified by the `timestep_units` argument when
`timesteps` is an iterable of float. Alternatively, units can be
specified for each step by passing an iterable of (value, unit)
tuples.
source_rates : float | Sequence[float]
Source rate in [neutron/sec] for each interval in `timesteps`.
timestep_units : {'s', 'min', 'h', 'd', 'a'}, optional
Units for values specified in the `timesteps` argument when passing
float values. 's' means seconds, 'min' means minutes, 'h' means
hours, 'd' means days, and 'a' means years (Julian).
output_dir : PathLike, optional
Path to directory where activation calculation outputs will be
saved.
operator_kwargs : dict, optional
Additional keyword arguments passed to
:class:`openmc.deplete.IndependentOperator`.
"""
if self.method == 'mesh-based':
# Get unique material for each (mesh, material) combination
mmv = self.results['mesh_material_volumes']
self.results['activation_materials'] = get_activation_materials(self.neutron_model, mmv)
else:
# Create unique material for each cell
activation_mats = openmc.Materials()
for cell in self.domains:
mat = cell.fill.clone()
mat.name = f'Cell {cell.id}'
mat.depletable = True
mat.volume = cell.volume
activation_mats.append(mat)
self.results['activation_materials'] = activation_mats
# Save activation materials to file
output_dir = Path(output_dir)
output_dir.mkdir(parents=True, exist_ok=True)
self.results['activation_materials'].export_to_xml(
output_dir / 'materials.xml')
# Create depletion operator for the activation materials
if operator_kwargs is None:
operator_kwargs = {}
operator_kwargs.setdefault('normalization_mode', 'source-rate')
op = IndependentOperator(
self.results['activation_materials'],
self.results['fluxes'],
self.results['micros'],
**operator_kwargs
)
# Create time integrator and solve depletion equations
integrator = PredictorIntegrator(
op, timesteps, source_rates=source_rates, timestep_units=timestep_units
)
output_path = output_dir / 'depletion_results.h5'
integrator.integrate(final_step=False, path=output_path)
comm.barrier()
# Get depletion results
self.results['depletion_results'] = Results(output_path)
def step3_photon_transport(
self,
time_indices: Sequence[int] | None = None,
bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
output_dir: PathLike = 'photon_transport',
mat_vol_kwargs: dict | None = None,
run_kwargs: dict | None = None,
):
"""Run the photon transport step.
This step performs photon transport calculations using decay photon
sources created from the activated materials. For each specified time,
it creates appropriate photon sources and runs a transport calculation.
In mesh-based mode, the sources are created using the mesh material
volumes, while in cell-based mode, they are created using bounding boxes
for each cell. This step will populate the 'photon_tallies' key in the
results dictionary.
Parameters
----------
time_indices : Sequence[int], optional
Sequence of time indices at which photon transport should be run;
represented as indices into the array of times formed by the
timesteps. For example, if two timesteps are specified, the array of
times would contain three entries, and [2] would indicate computing
photon results at the last time. A value of None indicates to run
photon transport for each time.
bounding_boxes : dict[int, openmc.BoundingBox], optional
Dictionary mapping cell IDs to bounding boxes used for spatial
source sampling in cell-based R2S calculations. Required if method
is 'cell-based'.
output_dir : PathLike, optional
Path to directory where photon transport outputs will be saved.
mat_vol_kwargs : dict, optional
Additional keyword arguments passed to
:meth:`openmc.MeshBase.material_volumes`.
run_kwargs : dict, optional
Additional keyword arguments passed to :meth:`openmc.Model.run`
during the photon transport step. By default, output is disabled.
"""
# TODO: Automatically determine bounding box for each cell
if bounding_boxes is None and self.method == 'cell-based':
raise ValueError("bounding_boxes must be provided for cell-based "
"R2S calculations.")
# Set default run arguments if not provided
if run_kwargs is None:
run_kwargs = {}
run_kwargs.setdefault('output', False)
# Write out JSON file with tally IDs that can be used for loading
# results
output_dir = Path(output_dir)
output_dir.mkdir(parents=True, exist_ok=True)
# Get default time indices if not provided
if time_indices is None:
n_steps = len(self.results['depletion_results'])
time_indices = list(range(n_steps))
# Check whether the photon model is different
neutron_univ = self.neutron_model.geometry.root_universe
photon_univ = self.photon_model.geometry.root_universe
different_photon_model = (neutron_univ != photon_univ)
# For mesh-based calculations, compute material volume fractions for the
# photon model if it is different from the neutron model to account for
# potential material changes
if self.method == 'mesh-based' and different_photon_model:
self.results['mesh_material_volumes_photon'] = photon_mmv = comm.bcast(
self.domains.material_volumes(self.photon_model, **mat_vol_kwargs))
# Save photon MMV results to file
if comm.rank == 0:
photon_mmv.save(output_dir / 'mesh_material_volumes.npz')
if comm.rank == 0:
tally_ids = [tally.id for tally in self.photon_model.tallies]
with open(output_dir / 'tally_ids.json', 'w') as f:
json.dump(tally_ids, f)
self.results['photon_tallies'] = {}
# Get dictionary of cells in the photon model
if different_photon_model:
photon_cells = self.photon_model.geometry.get_all_cells()
for time_index in time_indices:
# Create decay photon source
if self.method == 'mesh-based':
self.photon_model.settings.source = \
self.get_decay_photon_source_mesh(time_index)
else:
sources = []
results = self.results['depletion_results']
for cell, original_mat in zip(self.domains, self.results['activation_materials']):
# Skip if the cell is not in the photon model or the
# material has changed
if different_photon_model:
if cell.id not in photon_cells or \
cell.fill.id != photon_cells[cell.id].fill.id:
continue
# Get bounding box for the cell
bounding_box = bounding_boxes[cell.id]
# Get activated material composition
activated_mat = results[time_index].get_material(str(original_mat.id))
# Create decay photon source source
space = openmc.stats.Box(*bounding_box)
energy = activated_mat.get_decay_photon_energy()
strength = energy.integral() if energy is not None else 0.0
source = openmc.IndependentSource(
space=space,
energy=energy,
particle='photon',
strength=strength,
constraints={'domains': [cell]}
)
sources.append(source)
self.photon_model.settings.source = sources
# Convert time_index (which may be negative) to a normal index
if time_index < 0:
time_index = len(self.results['depletion_results']) + time_index
# Run photon transport calculation
photon_dir = Path(output_dir) / f'time_{time_index}'
with TemporarySession(self.photon_model, cwd=photon_dir):
statepoint_path = self.photon_model.run(**run_kwargs)
# Store tally results
with openmc.StatePoint(statepoint_path) as sp:
self.results['photon_tallies'][time_index] = [
sp.tallies[tally.id] for tally in self.photon_model.tallies
]
def get_decay_photon_source_mesh(
self,
time_index: int = -1
) -> list[openmc.MeshSource]:
"""Create decay photon source for a mesh-based calculation.
This function creates N :class:`MeshSource` objects where N is the
maximum number of unique materials that appears in a single mesh
element. For each mesh element-material combination, and
IndependentSource instance is created with a spatial constraint limited
the sampled decay photons to the correct region.
When the photon transport model is different from the neutron model, the
photon MeshMaterialVolumes is used to determine whether an (element,
material) combination exists in the photon model.
Parameters
----------
time_index : int, optional
Time index for the decay photon source. Default is -1 (last time).
Returns
-------
list of openmc.MeshSource
A list of MeshSource objects, each containing IndependentSource
instances for the decay photons in the corresponding mesh element.
"""
mat_dict = self.neutron_model._get_all_materials()
# Some MeshSource objects will have empty positions; create a "null source"
# that is used for this case
null_source = openmc.IndependentSource(particle='photon', strength=0.0)
# List to hold sources for each MeshSource (length = N)
source_lists = []
# Index in the overall list of activated materials
index_mat = 0
# Get various results from previous steps
mat_vols = self.results['mesh_material_volumes']
materials = self.results['activation_materials']
results = self.results['depletion_results']
photon_mat_vols = self.results.get('mesh_material_volumes_photon')
# Total number of mesh elements
n_elements = mat_vols.num_elements
for index_elem in range(n_elements):
# Determine which materials exist in the photon model for this element
if photon_mat_vols is not None:
photon_materials = {
mat_id
for mat_id, _ in photon_mat_vols.by_element(index_elem)
if mat_id is not None
}
for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)):
# Skip void volume
if mat_id is None:
continue
# Skip if this material doesn't exist in photon model
if photon_mat_vols is not None and mat_id not in photon_materials:
index_mat += 1
continue
# Check whether a new MeshSource object is needed
if j >= len(source_lists):
source_lists.append([null_source]*n_elements)
# Get activated material composition
original_mat = materials[index_mat]
activated_mat = results[time_index].get_material(str(original_mat.id))
# Create decay photon source source
energy = activated_mat.get_decay_photon_energy()
if energy is not None:
strength = energy.integral()
source_lists[j][index_elem] = openmc.IndependentSource(
energy=energy,
particle='photon',
strength=strength,
constraints={'domains': [mat_dict[mat_id]]}
)
# Increment index of activated material
index_mat += 1
# Return list of mesh sources
return [openmc.MeshSource(self.domains, sources) for sources in source_lists]
def load_results(self, path: PathLike):
"""Load results from a previous R2S calculation.
Parameters
----------
path : PathLike
Path to the directory containing the R2S calculation results.
"""
path = Path(path)
# Load neutron transport results
neutron_dir = path / 'neutron_transport'
if self.method == 'mesh-based':
mmv_file = neutron_dir / 'mesh_material_volumes.npz'
if mmv_file.exists():
self.results['mesh_material_volumes'] = \
openmc.MeshMaterialVolumes.from_npz(mmv_file)
fluxes_file = neutron_dir / 'fluxes.npy'
if fluxes_file.exists():
self.results['fluxes'] = list(np.load(fluxes_file, allow_pickle=True))
micros_dict = read_microxs_hdf5(neutron_dir / 'micros.h5')
self.results['micros'] = [
micros_dict[f'domain_{i}'] for i in range(len(micros_dict))
]
# Load activation results
activation_dir = path / 'activation'
activation_results = activation_dir / 'depletion_results.h5'
if activation_results.exists():
self.results['depletion_results'] = Results(activation_results)
activation_mats_file = activation_dir / 'materials.xml'
if activation_mats_file.exists():
self.results['activation_materials'] = \
openmc.Materials.from_xml(activation_mats_file)
# Load photon transport results
photon_dir = path / 'photon_transport'
# Load photon mesh material volumes if they exist (for mesh-based calculations)
if self.method == 'mesh-based':
photon_mmv_file = photon_dir / 'mesh_material_volumes.npz'
if photon_mmv_file.exists():
self.results['mesh_material_volumes_photon'] = \
openmc.MeshMaterialVolumes.from_npz(photon_mmv_file)
# Load tally IDs from JSON file
tally_ids_path = photon_dir / 'tally_ids.json'
if tally_ids_path.exists():
with tally_ids_path.open('r') as f:
tally_ids = json.load(f)
self.results['photon_tallies'] = {}
# For each photon transport calc, load the statepoint and get the
# tally results based on tally_ids
for time_dir in photon_dir.glob('time_*'):
time_index = int(time_dir.name.split('_')[1])
for sp_path in time_dir.glob('statepoint.*.h5'):
with openmc.StatePoint(sp_path) as sp:
self.results['photon_tallies'][time_index] = [
sp.tallies[tally_id] for tally_id in tally_ids
]

View file

@ -203,7 +203,7 @@ class Results(list):
# Evaluate value in each region
for i, result in enumerate(self):
times[i] = result.time[0]
concentrations[i] = result[0, mat_id, nuc]
concentrations[i] = result[mat_id, nuc]
# Unit conversions
times = _get_time_as(times, time_units)
@ -363,7 +363,7 @@ class Results(list):
# Evaluate value in each region
for i, result in enumerate(self):
times[i] = result.time[0]
rates[i] = result.rates[0].get(mat_id, nuc, rx) * result[0, mat, nuc]
rates[i] = result.rates.get(mat_id, nuc, rx) * result[mat, nuc]
return times, rates
@ -397,7 +397,7 @@ class Results(list):
# Get time/eigenvalue at each point
for i, result in enumerate(self):
times[i] = result.time[0]
eigenvalues[i] = result.k[0]
eigenvalues[i] = result.k
# Convert time units if necessary
times = _get_time_as(times, time_units)
@ -630,7 +630,7 @@ class Results(list):
for nuc in result.index_nuc:
if nuc not in available_cross_sections:
continue
atoms = result[0, mat_id, nuc]
atoms = result[mat_id, nuc]
if atoms > 0.0:
atoms_per_barn_cm = 1e-24 * atoms / mat.volume
mat.remove_nuclide(nuc) # Replace if it's there

View file

@ -16,7 +16,7 @@ from openmc.mpi import comm, MPI
from openmc.checkvalue import PathLike
from .reaction_rates import ReactionRates
VERSION_RESULTS = (1, 1)
VERSION_RESULTS = (1, 2)
__all__ = ["StepResult"]
@ -30,8 +30,8 @@ class StepResult:
Attributes
----------
k : list of (float, float)
Eigenvalue and uncertainty for each substep.
k : tuple of (float, float)
Eigenvalue and uncertainty at end of step.
time : list of float
Time at beginning, end of step, in seconds.
source_rate : float
@ -40,8 +40,8 @@ class StepResult:
Number of mats.
n_nuc : int
Number of nuclides.
rates : list of ReactionRates
The reaction rates for each substep.
rates : ReactionRates
The reaction rates at end of step.
volume : dict of str to float
Dictionary mapping mat id to volume.
index_mat : dict of str to int
@ -52,10 +52,8 @@ class StepResult:
A dictionary mapping mat ID as string to global index.
n_hdf5_mats : int
Number of materials in entire geometry.
n_stages : int
Number of stages in simulation.
data : numpy.ndarray
Atom quantity, stored by stage, mat, then by nuclide.
Atom quantity, stored by mat, then by nuclide.
proc_time : int
Average time spent depleting a material across all
materials and processes
@ -86,17 +84,17 @@ class StepResult:
Parameters
----------
pos : tuple
A three-length tuple containing a stage index, mat index and a nuc
index. All can be integers or slices. The second two can be
A two-length tuple containing a mat index and a nuc
index. Both can be integers or slices, or can be
strings corresponding to their respective dictionary.
Returns
-------
float
The atoms for stage, mat, nuc
The atoms for mat, nuc
"""
stage, mat, nuc = pos
mat, nuc = pos
if isinstance(mat, openmc.Material):
mat = str(mat.id)
if isinstance(mat, str):
@ -104,7 +102,7 @@ class StepResult:
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
return self.data[stage, mat, nuc]
return self.data[mat, nuc]
def __setitem__(self, pos, val):
"""Sets an item from results.
@ -112,21 +110,21 @@ class StepResult:
Parameters
----------
pos : tuple
A three-length tuple containing a stage index, mat index and a nuc
index. All can be integers or slices. The second two can be
A two-length tuple containing a mat index and a nuc
index. Both can be integers or slices, or can be
strings corresponding to their respective dictionary.
val : float
The value to set data to.
"""
stage, mat, nuc = pos
mat, nuc = pos
if isinstance(mat, str):
mat = self.index_mat[mat]
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
self.data[stage, mat, nuc] = val
self.data[mat, nuc] = val
@property
def n_mat(self):
@ -140,11 +138,7 @@ class StepResult:
def n_hdf5_mats(self):
return len(self.mat_to_hdf5_ind)
@property
def n_stages(self):
return self.data.shape[0]
def allocate(self, volume, nuc_list, burn_list, full_burn_list, stages):
def allocate(self, volume, nuc_list, burn_list, full_burn_list):
"""Allocate memory for depletion step data
Parameters
@ -157,8 +151,6 @@ class StepResult:
A list of all mat IDs to be burned. Used for sorting the simulation.
full_burn_list : list of str
List of all burnable material IDs
stages : int
Number of stages in simulation.
"""
self.volume = copy.deepcopy(volume)
@ -167,7 +159,7 @@ class StepResult:
self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)}
# Create storage array
self.data = np.zeros((stages, self.n_mat, self.n_nuc))
self.data = np.zeros((self.n_mat, self.n_nuc))
def distribute(self, local_materials, ranges):
"""Create a new object containing data for distributed materials
@ -196,8 +188,8 @@ class StepResult:
for attr in direct_attrs:
setattr(new, attr, getattr(self, attr))
# Get applicable slice of data
new.data = self.data[:, ranges]
new.rates = [r[ranges] for r in self.rates]
new.data = self.data[ranges]
new.rates = self.rates[ranges]
return new
def get_material(self, mat_id):
@ -232,7 +224,7 @@ class StepResult:
f'values are {list(self.volume.keys())}'
) from e
for nuc, _ in sorted(self.index_nuc.items(), key=lambda x: x[1]):
atoms = self[0, mat_id, nuc]
atoms = self[mat_id, nuc]
if atoms <= 0.0:
continue
atom_per_bcm = atoms / vol * 1e-24
@ -240,7 +232,7 @@ class StepResult:
material.volume = vol
return material
def export_to_hdf5(self, filename, step):
def export_to_hdf5(self, filename, step, write_rates: bool = False):
"""Export results to an HDF5 file
Parameters
@ -249,6 +241,8 @@ class StepResult:
The filename to write to
step : int
What step is this?
write_rates : bool, optional
Whether to include reaction rate datasets in the results file.
"""
# Write new file if first time step, else add to existing file
@ -259,7 +253,8 @@ class StepResult:
kwargs['driver'] = 'mpio'
kwargs['comm'] = comm
with h5py.File(filename, **kwargs) as handle:
self._to_hdf5(handle, step, parallel=True)
self._to_hdf5(handle, step, parallel=True,
write_rates=write_rates)
else:
# Gather results at root process
all_results = comm.gather(self)
@ -268,15 +263,18 @@ class StepResult:
if comm.rank == 0:
with h5py.File(filename, **kwargs) as handle:
for res in all_results:
res._to_hdf5(handle, step, parallel=False)
res._to_hdf5(handle, step, parallel=False,
write_rates=write_rates)
def _write_hdf5_metadata(self, handle):
def _write_hdf5_metadata(self, handle, write_rates):
"""Writes result metadata in HDF5 file
Parameters
----------
handle : h5py.File or h5py.Group
An hdf5 file or group type to store this in.
write_rates : bool
Whether reaction rate datasets are being written.
"""
# Create and save the 5 dictionaries:
@ -284,8 +282,8 @@ class StepResult:
# self.index_mat -> self.volume (TODO: support for changing volumes)
# self.index_nuc
# reactions
# self.rates[0].index_nuc (can be different from above, above is superset)
# self.rates[0].index_rx
# self.rates.index_nuc (can be different from above, above is superset)
# self.rates.index_rx
# these are shared by every step of the simulation, and should be deduplicated.
# Store concentration mat and nuclide dictionaries (along with volumes)
@ -295,13 +293,19 @@ class StepResult:
mat_list = sorted(self.mat_to_hdf5_ind, key=int)
nuc_list = sorted(self.index_nuc)
rxn_list = sorted(self.rates[0].index_rx)
include_rates = (
write_rates
and self.rates is not None
and bool(self.rates.index_nuc)
and bool(self.rates.index_rx)
)
rxn_list = sorted(self.rates.index_rx) if include_rates else []
n_mats = self.n_hdf5_mats
n_nuc_number = len(nuc_list)
n_nuc_rxn = len(self.rates[0].index_nuc)
n_nuc_rxn = len(self.rates.index_nuc) if include_rates else 0
n_rxn = len(rxn_list)
n_stages = self.n_stages
mat_group = handle.create_group("materials")
@ -315,41 +319,44 @@ class StepResult:
for nuc in nuc_list:
nuc_single_group = nuc_group.create_group(nuc)
nuc_single_group.attrs["atom number index"] = self.index_nuc[nuc]
if nuc in self.rates[0].index_nuc:
nuc_single_group.attrs["reaction rate index"] = self.rates[0].index_nuc[nuc]
if include_rates and nuc in self.rates.index_nuc:
nuc_single_group.attrs["reaction rate index"] = (
self.rates.index_nuc[nuc])
rxn_group = handle.create_group("reactions")
if include_rates:
rxn_group = handle.create_group("reactions")
for rxn in rxn_list:
rxn_single_group = rxn_group.create_group(rxn)
rxn_single_group.attrs["index"] = self.rates[0].index_rx[rxn]
for rxn in rxn_list:
rxn_single_group = rxn_group.create_group(rxn)
rxn_single_group.attrs["index"] = (
self.rates.index_rx[rxn])
# Construct array storage
handle.create_dataset("number", (1, n_stages, n_mats, n_nuc_number),
maxshape=(None, n_stages, n_mats, n_nuc_number),
handle.create_dataset("number", (1, n_mats, n_nuc_number),
maxshape=(None, n_mats, n_nuc_number),
chunks=True,
dtype='float64')
if n_nuc_rxn > 0 and n_rxn > 0:
handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn),
maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn),
chunks=True,
dtype='float64')
if include_rates and n_nuc_rxn > 0 and n_rxn > 0:
handle.create_dataset(
"reaction rates", (1, n_mats, n_nuc_rxn, n_rxn),
maxshape=(None, n_mats, n_nuc_rxn, n_rxn),
chunks=True, dtype='float64')
handle.create_dataset("eigenvalues", (1, n_stages, 2),
maxshape=(None, n_stages, 2), dtype='float64')
handle.create_dataset("eigenvalues", (1, 2),
maxshape=(None, 2), dtype='float64')
handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
handle.create_dataset("source_rate", (1, n_stages), maxshape=(None, n_stages),
handle.create_dataset("source_rate", (1,), maxshape=(None,),
dtype='float64')
handle.create_dataset(
"depletion time", (1,), maxshape=(None,),
dtype="float64")
def _to_hdf5(self, handle, index, parallel=False):
def _to_hdf5(self, handle, index, parallel=False, write_rates: bool = False):
"""Converts results object into an hdf5 object.
Parameters
@ -360,12 +367,14 @@ class StepResult:
What step is this?
parallel : bool
Being called with parallel HDF5?
write_rates : bool, optional
Whether reaction rate datasets are being written.
"""
if "/number" not in handle:
if parallel:
comm.barrier()
self._write_hdf5_metadata(handle)
self._write_hdf5_metadata(handle, write_rates)
if parallel:
comm.barrier()
@ -417,18 +426,14 @@ class StepResult:
return
# Add data
# Note, for the last step, self.n_stages = 1, even if n_stages != 1.
n_stages = self.n_stages
inds = [self.mat_to_hdf5_ind[mat] for mat in self.index_mat]
low = min(inds)
high = max(inds)
for i in range(n_stages):
number_dset[index, i, low:high+1] = self.data[i]
if has_reactions:
rxn_dset[index, i, low:high+1] = self.rates[i]
if comm.rank == 0:
eigenvalues_dset[index, i] = self.k[i]
number_dset[index, low:high+1] = self.data
if has_reactions:
rxn_dset[index, low:high+1] = self.rates
if comm.rank == 0:
eigenvalues_dset[index] = self.k
time_dset[index] = self.time
source_rate_dset[index] = self.source_rate
if self.proc_time is not None:
@ -459,10 +464,24 @@ class StepResult:
# Older versions used "power" instead of "source_rate"
source_rate_dset = handle["/power"]
results.data = number_dset[step, :, :, :]
results.k = eigenvalues_dset[step, :]
# Check if this is an old format file (with stages dimension) or new format
# Old format: number has shape (n_steps, n_stages, n_mats, n_nucs)
# New format: number has shape (n_steps, n_mats, n_nucs)
has_stages = len(number_dset.shape) == 4
if has_stages:
# Old format - extract data from first stage (index 0)
results.data = number_dset[step, 0, :, :]
results.k = eigenvalues_dset[step, 0, :]
# source_rate had shape (n_steps, n_stages) in old format
results.source_rate = source_rate_dset[step, 0]
else:
# New format - no stages dimension
results.data = number_dset[step, :, :]
results.k = eigenvalues_dset[step, :]
results.source_rate = source_rate_dset[step]
results.time = time_dset[step, :]
results.source_rate = source_rate_dset[step, 0]
if "depletion time" in handle:
proc_time_dset = handle["/depletion time"]
@ -493,33 +512,45 @@ class StepResult:
if "reaction rate index" in nuc_handle.attrs:
rxn_nuc_to_ind[nuc] = nuc_handle.attrs["reaction rate index"]
for rxn, rxn_handle in handle["/reactions"].items():
rxn_to_ind[rxn] = rxn_handle.attrs["index"]
if "reactions" in handle:
for rxn, rxn_handle in handle["/reactions"].items():
rxn_to_ind[rxn] = rxn_handle.attrs["index"]
results.rates = []
# Reconstruct reactions
for i in range(results.n_stages):
rate = ReactionRates(results.index_mat, rxn_nuc_to_ind, rxn_to_ind, True)
if "reaction rates" in handle:
rate[:] = handle["/reaction rates"][step, i, :, :, :]
results.rates.append(rate)
# Reconstruct reaction rates
rate = ReactionRates(results.index_mat, rxn_nuc_to_ind, rxn_to_ind, True)
if "reaction rates" in handle:
if has_stages:
# Old format: (n_steps, n_stages, n_mats, n_nucs, n_rxns)
rate[:] = handle["/reaction rates"][step, 0, :, :, :]
else:
# New format: (n_steps, n_mats, n_nucs, n_rxns)
rate[:] = handle["/reaction rates"][step, :, :, :]
results.rates = rate
return results
@staticmethod
def save(op, x, op_results, t, source_rate, step_ind, proc_time=None,
path: PathLike = "depletion_results.h5"):
def save(
op,
x,
op_results,
t,
source_rate,
step_ind,
proc_time=None,
write_rates: bool = False,
path: PathLike = "depletion_results.h5"
):
"""Creates and writes depletion results to disk
Parameters
----------
op : openmc.deplete.abc.TransportOperator
The operator used to generate these results.
x : list of list of numpy.array
The prior x vectors. Indexed [i][cell] using the above equation.
op_results : list of openmc.deplete.OperatorResult
Results of applying transport operator
x : numpy.array
End-of-step concentrations for each material
op_results : openmc.deplete.OperatorResult
Result of applying transport operator at end of step
t : list of float
Time indices.
source_rate : float
@ -530,7 +561,8 @@ class StepResult:
Total process time spent depleting materials. This may
be process-dependent and will be reduced across MPI
processes.
write_rates : bool, optional
Whether reaction rates should be written to the results file.
path : PathLike
Path to file to write. Defaults to 'depletion_results.h5'.
@ -539,26 +571,20 @@ class StepResult:
# Get indexing terms
vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
stages = len(x)
# Create results
results = StepResult()
results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages)
results.allocate(vol_dict, nuc_list, burn_list, full_burn_list)
n_mat = len(burn_list)
for i in range(stages):
for mat_i in range(n_mat):
results[i, mat_i, :] = x[i][mat_i]
for mat_i in range(n_mat):
results[mat_i, :] = x[mat_i]
ks = []
for r in op_results:
if isinstance(r.k, type(None)):
ks += [(None, None)]
else:
ks += [(r.k.nominal_value, r.k.std_dev)]
results.k = ks
results.rates = [r.rates for r in op_results]
if isinstance(op_results.k, type(None)):
results.k = (None, None)
else:
results.k = (op_results.k.nominal_value, op_results.k.std_dev)
results.rates = op_results.rates
results.time = t
results.source_rate = source_rate
results.proc_time = proc_time
@ -567,7 +593,7 @@ class StepResult:
if not Path(path).is_file():
Path(path).parent.mkdir(parents=True, exist_ok=True)
results.export_to_hdf5(path, step_ind)
results.export_to_hdf5(path, step_ind, write_rates)
def transfer_volumes(self, model):
"""Transfers volumes from depletion results to geometry

View file

@ -49,9 +49,10 @@ class ExternalRates:
self.local_mats = operator.local_mats
self.number_of_timesteps = number_of_timesteps
# initialize transfer rates container dict
#initialize transfer rates container dict
self.external_rates = {mat: defaultdict(list) for mat in self.burnable_mats}
self.external_timesteps = []
self.redox = {}
def _get_material_id(self, val):
"""Helper method for getting material id from Material obj or name.
@ -300,6 +301,46 @@ class TransferRates(ExternalRates):
self.external_timesteps = np.unique(np.concatenate(
[self.external_timesteps, timesteps]))
def set_redox(self, material, buffer, oxidation_states, timesteps=None):
"""Add redox control to depletable material.
Parameters
----------
material : openmc.Material or str or int
Depletable material
buffer : dict
Dictionary of buffer nuclides used to maintain redox balance.
Keys are nuclide names (strings) and values are their respective
fractions (float) that collectively sum to 1.
oxidation_states : dict
User-defined oxidation states for elements.
Keys are element symbols (e.g., 'H', 'He'), and values are their
corresponding oxidation states as integers (e.g., +1, 0).
timesteps : list of int, optional
List of timestep indices where to set external source rates.
Defaults to None, which means the external source rate is set for
all timesteps.
"""
material_id = self._get_material_id(material)
if timesteps is not None:
for timestep in timesteps:
check_value('timestep', timestep, range(self.number_of_timesteps))
timesteps = np.array(timesteps)
else:
timesteps = np.arange(self.number_of_timesteps)
#Check nuclides in buffer exist
for nuc in buffer:
if nuc not in self.chain_nuclides:
raise ValueError(f'{nuc} is not a valid nuclide.')
# Checks element in oxidation states exist
for elm in oxidation_states:
if elm not in ELEMENT_SYMBOL.values():
raise ValueError(f'{elm} is not a valid element.')
self.redox[material_id] = (buffer, oxidation_states)
self.external_timesteps = np.unique(np.concatenate(
[self.external_timesteps, timesteps]))
class ExternalSourceRates(ExternalRates):
"""Class for defining external source rates.
@ -366,7 +407,7 @@ class ExternalSourceRates(ExternalRates):
rate : float
External source rate in units of mass per time. A positive or
negative value corresponds to a feed or removal rate, respectively.
units : {'g/s', 'g/min', 'g/h', 'g/d', 'g/a'}
rate_units : {'g/s', 'g/min', 'g/h', 'g/d', 'g/a'}
Units for values specified in the `rate` argument. 's' for seconds,
'min' for minutes, 'h' for hours, 'a' for Julian years.
timesteps : list of int, optional

View file

@ -83,7 +83,7 @@ def pwr_pin_cell() -> openmc.Model:
constraints={'fissionable': True}
)
plot = openmc.Plot.from_geometry(model.geometry)
plot = openmc.SlicePlot.from_geometry(model.geometry)
plot.pixels = (300, 300)
plot.color_by = 'material'
model.plots.append(plot)
@ -429,7 +429,7 @@ def pwr_core() -> openmc.Model:
model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
[-160, -160, -183], [160, 160, 183]))
plot = openmc.Plot()
plot = openmc.SlicePlot()
plot.origin = (125, 125, 0)
plot.width = (250, 250)
plot.pixels = (3000, 3000)
@ -544,7 +544,7 @@ def pwr_assembly() -> openmc.Model:
constraints={'fissionable': True}
)
plot = openmc.Plot()
plot = openmc.SlicePlot()
plot.origin = (0.0, 0.0, 0)
plot.width = (21.42, 21.42)
plot.pixels = (300, 300)

View file

@ -164,7 +164,7 @@ def plot_inline(plots, openmc_exec='openmc', cwd='.', path_input=None):
Parameters
----------
plots : Iterable of openmc.Plot
plots : Iterable of openmc.PlotBase
Plots to display
openmc_exec : str
Path to OpenMC executable

View file

@ -1839,12 +1839,21 @@ class DistribcellFilter(Filter):
@property
def paths(self):
return self._paths
if self._paths is None:
if not hasattr(self, '_geometry'):
raise ValueError(
"Model must be exported before the 'paths' attribute is" \
"available for a DistribcellFilter.")
@paths.setter
def paths(self, paths):
cv.check_iterable_type('paths', paths, str)
self._paths = paths
# Determine paths for cell instances
self._geometry.determine_paths()
# Get paths for the corresponding cell
cell_id = self.bins[0]
cell = self._geometry.get_all_cells()[cell_id]
self._paths = cell.paths
return self._paths
@Filter.bins.setter
def bins(self, bins):

View file

@ -34,6 +34,10 @@ _dll.openmc_cell_get_temperature.argtypes = [
c_int32, POINTER(c_int32), POINTER(c_double)]
_dll.openmc_cell_get_temperature.restype = c_int
_dll.openmc_cell_get_temperature.errcheck = _error_handler
_dll.openmc_cell_get_density.argtypes = [
c_int32, POINTER(c_int32), POINTER(c_double)]
_dll.openmc_cell_get_density.restype = c_int
_dll.openmc_cell_get_density.errcheck = _error_handler
_dll.openmc_cell_get_name.argtypes = [c_int32, POINTER(c_char_p)]
_dll.openmc_cell_get_name.restype = c_int
_dll.openmc_cell_get_name.errcheck = _error_handler
@ -58,6 +62,10 @@ _dll.openmc_cell_set_temperature.argtypes = [
c_int32, c_double, POINTER(c_int32), c_bool]
_dll.openmc_cell_set_temperature.restype = c_int
_dll.openmc_cell_set_temperature.errcheck = _error_handler
_dll.openmc_cell_set_density.argtypes = [
c_int32, c_double, POINTER(c_int32), c_bool]
_dll.openmc_cell_set_density.restype = c_int
_dll.openmc_cell_set_density.errcheck = _error_handler
_dll.openmc_cell_set_translation.argtypes = [c_int32, POINTER(c_double)]
_dll.openmc_cell_set_translation.restype = c_int
_dll.openmc_cell_set_translation.errcheck = _error_handler
@ -236,6 +244,44 @@ class Cell(_FortranObjectWithID):
_dll.openmc_cell_set_temperature(self._index, T, instance, set_contained)
def get_density(self, instance: int | None = None):
"""Get the density of a cell in [g/cm3]
Parameters
----------
instance : int or None
Which instance of the cell
"""
if instance is not None:
instance = c_int32(instance)
rho = c_double()
_dll.openmc_cell_get_density(self._index, instance, rho)
return rho.value
def set_density(self, rho: float, instance: int | None = None,
set_contained: bool = False):
"""Set the density of a cell
Parameters
----------
rho : float
Density of the cell in [g/cm3]
instance : int or None
Which instance of the cell
set_contained : bool
If cell is not filled by a material, whether to set the density
of all filled cells
"""
if instance is not None:
instance = c_int32(instance)
_dll.openmc_cell_set_density(self._index, rho, instance, set_contained)
@property
def translation(self):
translation = np.zeros(3)

View file

@ -7,12 +7,14 @@ import os
from pathlib import Path
from random import getrandbits
from tempfile import TemporaryDirectory
import traceback as tb
import numpy as np
from numpy.ctypeslib import as_array
from . import _dll
from .error import _error_handler
from ..mpi import comm
from openmc.checkvalue import PathLike
import openmc.lib
import openmc
@ -632,6 +634,9 @@ class TemporarySession:
model : openmc.Model, optional
OpenMC model to use for the session. If None, a minimal working model is
created.
cwd : PathLike, optional
Working directory in which to run OpenMC. If None, a temporary directory
is created and deleted automatically.
**init_kwargs
Keyword arguments to pass to :func:`openmc.lib.init`.
@ -639,10 +644,13 @@ class TemporarySession:
----------
model : openmc.Model
The OpenMC model used for the session.
comm : mpi4py.MPI.Intracomm
The MPI intracommunicator used for the session.
"""
def __init__(self, model=None, **init_kwargs):
self.init_kwargs = init_kwargs
def __init__(self, model=None, cwd=None, **init_kwargs):
self.init_kwargs = dict(init_kwargs)
self.cwd = cwd
if model is None:
surf = openmc.Sphere(boundary_type="vacuum")
cell = openmc.Cell(region=-surf)
@ -652,6 +660,10 @@ class TemporarySession:
particles=1, batches=1, output={'summary': False})
self.model = model
# Determine MPI intercommunicator
self.init_kwargs.setdefault('intracomm', comm)
self.comm = self.init_kwargs['intracomm']
def __enter__(self):
"""Initialize the OpenMC library in a temporary directory."""
# If already initialized, the context manager is a no-op
@ -662,14 +674,24 @@ class TemporarySession:
# Store original working directory
self.orig_dir = Path.cwd()
# Set up temporary directory
self.tmp_dir = TemporaryDirectory()
working_dir = Path(self.tmp_dir.name)
working_dir.mkdir(parents=True, exist_ok=True)
os.chdir(working_dir)
if self.cwd is None:
# Set up temporary directory on rank 0
if self.comm.rank == 0:
self._tmp_dir = TemporaryDirectory()
self.cwd = self._tmp_dir.name
# Export model and initialize OpenMC
self.model.export_to_model_xml()
# Broadcast the path so that all ranks use the same directory
self.cwd = self.comm.bcast(self.cwd)
# Create and change to specified directory
self.cwd = Path(self.cwd)
self.cwd.mkdir(parents=True, exist_ok=True)
os.chdir(self.cwd)
# Export model on first rank and initialize OpenMC
if self.comm.rank == 0:
self.model.export_to_model_xml()
self.comm.barrier()
openmc.lib.init(**self.init_kwargs)
return self
@ -679,11 +701,24 @@ class TemporarySession:
if self.already_initialized:
return
# If an exception occurred, abort all ranks immediately
if exc_type is not None:
# Print exception info on the rank that failed
tb.print_exception(exc_type, exc_value, traceback)
sys.stdout.flush()
# Abort all MPI processes
self.comm.Abort(1)
try:
finalize()
finally:
os.chdir(self.orig_dir)
self.tmp_dir.cleanup()
# Make sure all ranks have finalized before deleting temporary dir
self.comm.barrier()
if hasattr(self, '_tmp_dir'):
self._tmp_dir.cleanup()
class _DLLGlobal:

View file

@ -52,7 +52,7 @@ class _PlotBase(Structure):
C-Type Attributes
-----------------
origin : openmc.lib.plot._Position
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
@ -60,6 +60,8 @@ class _PlotBase(Structure):
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
@ -187,14 +189,6 @@ class _PlotBase(Structure):
def color_overlaps(self, color_overlaps):
self.color_overlaps_ = color_overlaps
@property
def color_overlaps(self):
return self.color_overlaps_
@color_overlaps.setter
def color_overlaps(self, val):
self.color_overlaps_ = val
def __repr__(self):
out_str = ["-----",
"Plot:",

View file

@ -60,6 +60,26 @@ class Material(IDManagerMixin):
temperature : float, optional
Temperature of the material in Kelvin. If not specified, the material
inherits the default temperature applied to the model.
density : float, optional
Density of the material (units defined separately)
density_units : str
Units used for `density`. Can be one of 'g/cm3', 'g/cc', 'kg/m3',
'atom/b-cm', 'atom/cm3', 'sum', or 'macro'. The 'macro' unit only
applies in the case of a multi-group calculation. Defaults to 'sum'.
depletable : bool, optional
Indicate whether the material is depletable. Defaults to False.
volume : float, optional
Volume of the material in cm^3. This can either be set manually or
calculated in a stochastic volume calculation and added via the
:meth:`Material.add_volume_information` method.
components : dict of str to float or dict
Dictionary mapping element or nuclide names to their atom or weight
percent. To specify enrichment of an element, the entry of
``components`` for that element must instead be a dictionary containing
the keyword arguments as well as a value for ``'percent'``
percent_type : {'ao', 'wo'}
Whether the values in `components` should be interpreted as atom percent
('ao') or weight percent ('wo').
Attributes
----------
@ -111,17 +131,28 @@ class Material(IDManagerMixin):
next_id = 1
used_ids = set()
def __init__(self, material_id=None, name='', temperature=None):
def __init__(
self,
material_id: int | None = None,
name: str = "",
temperature: float | None = None,
density: float | None = None,
density_units: str = "sum",
depletable: bool | None = False,
volume: float | None = None,
components: dict | None = None,
percent_type: str = "ao",
):
# Initialize class attributes
self.id = material_id
self.name = name
self.temperature = temperature
self._density = None
self._density_units = 'sum'
self._depletable = False
self._density_units = density_units
self._depletable = depletable
self._paths = None
self._num_instances = None
self._volume = None
self._volume = volume
self._atoms = {}
self._isotropic = []
self._ncrystal_cfg = None
@ -136,6 +167,15 @@ class Material(IDManagerMixin):
# If specified, a list of table names
self._sab = []
# Set density if provided
if density is not None:
self.set_density(density_units, density)
# Add components if provided
if components is not None:
self.add_components(components, percent_type=percent_type)
def __repr__(self) -> str:
string = 'Material\n'
string += '{: <16}=\t{}\n'.format('\tID', self._id)
@ -290,11 +330,13 @@ class Material(IDManagerMixin):
return self.get_decay_photon_energy(0.0)
def get_decay_photon_energy(
self,
clip_tolerance: float = 1e-6,
units: str = 'Bq',
volume: float | None = None
) -> Univariate | None:
self,
clip_tolerance: float = 1e-6,
units: str = 'Bq',
volume: float | None = None,
exclude_nuclides: list[str] | None = None,
include_nuclides: list[str] | None = None
) -> Univariate | None:
r"""Return energy distribution of decay photons from unstable nuclides.
.. versionadded:: 0.14.0
@ -302,22 +344,31 @@ class Material(IDManagerMixin):
Parameters
----------
clip_tolerance : float
Maximum fraction of :math:`\sum_i x_i p_i` for discrete
distributions that will be discarded.
Maximum fraction of :math:`\sum_i x_i p_i` for discrete distributions
that will be discarded.
units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3'}
Specifies the units on the integral of the distribution.
volume : float, optional
Volume of the material. If not passed, defaults to using the
:attr:`Material.volume` attribute.
exclude_nuclides : list of str, optional
Nuclides to exclude from the photon source calculation.
include_nuclides : list of str, optional
Nuclides to include in the photon source calculation. If specified,
only these nuclides are used.
Returns
-------
Univariate or None
Decay photon energy distribution. The integral of this distribution
is the total intensity of the photon source in the requested units.
Decay photon energy distribution. The integral of this distribution is
the total intensity of the photon source in the requested units.
"""
cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3'})
if exclude_nuclides is not None and include_nuclides is not None:
raise ValueError("Cannot specify both exclude_nuclides and include_nuclides")
if units == 'Bq':
multiplier = volume if volume is not None else self.volume
if multiplier is None:
@ -332,6 +383,11 @@ class Material(IDManagerMixin):
dists = []
probs = []
for nuc, atoms_per_bcm in self.get_nuclide_atom_densities().items():
if exclude_nuclides is not None and nuc in exclude_nuclides:
continue
if include_nuclides is not None and nuc not in include_nuclides:
continue
source_per_atom = openmc.data.decay_photon_energy(nuc)
if source_per_atom is not None and atoms_per_bcm > 0.0:
dists.append(source_per_atom)

View file

@ -5,6 +5,7 @@ from collections.abc import Iterable, Sequence, Mapping
from functools import wraps
from math import pi, sqrt, atan2
from numbers import Integral, Real
from pathlib import Path
from typing import Protocol
import h5py
@ -287,6 +288,7 @@ class MeshBase(IDManagerMixin, ABC):
model: openmc.Model,
n_samples: int | tuple[int, int, int] = 10_000,
include_void: bool = True,
material_volumes: MeshMaterialVolumes | None = None,
**kwargs
) -> list[openmc.Material]:
"""Generate homogenized materials over each element in a mesh.
@ -305,8 +307,12 @@ class MeshBase(IDManagerMixin, ABC):
the x, y, and z dimensions.
include_void : bool, optional
Whether homogenization should include voids.
material_volumes : MeshMaterialVolumes, optional
Previously computed mesh material volumes to use for homogenization.
If not provided, they will be computed by calling
:meth:`material_volumes`.
**kwargs
Keyword-arguments passed to :meth:`MeshBase.material_volumes`.
Keyword-arguments passed to :meth:`material_volumes`.
Returns
-------
@ -314,23 +320,16 @@ class MeshBase(IDManagerMixin, ABC):
Homogenized material in each mesh element
"""
vols = self.material_volumes(model, n_samples, **kwargs)
if material_volumes is None:
vols = self.material_volumes(model, n_samples, **kwargs)
else:
vols = material_volumes
mat_volume_by_element = [vols.by_element(i) for i in range(vols.num_elements)]
# Get dictionary of all materials
materials = model._get_all_materials()
# Create homogenized material for each element
materials = model.geometry.get_all_materials()
# Account for materials in DAGMC universes
# TODO: This should really get incorporated in lower-level calls to
# get_all_materials, but right now it requires information from the
# Model object
for cell in model.geometry.get_all_cells().values():
if isinstance(cell.fill, openmc.DAGMCUniverse):
names = cell.fill.material_names
materials.update({
mat.id: mat for mat in model.materials if mat.name in names
})
homogenized_materials = []
for mat_volume_list in mat_volume_by_element:
material_ids, volumes = [list(x) for x in zip(*mat_volume_list)]
@ -402,7 +401,7 @@ class MeshBase(IDManagerMixin, ABC):
# In order to get mesh into model, we temporarily replace the
# tallies with a single mesh tally using the current mesh
original_tallies = model.tallies
original_tallies = list(model.tallies)
new_tally = openmc.Tally()
new_tally.filters = [openmc.MeshFilter(self)]
new_tally.scores = ['flux']
@ -424,7 +423,6 @@ class MeshBase(IDManagerMixin, ABC):
# Restore original tallies
model.tallies = original_tallies
return volumes
@ -2446,6 +2444,7 @@ class UnstructuredMesh(MeshBase):
_UNSUPPORTED_ELEM = -1
_LINEAR_TET = 0
_LINEAR_HEX = 1
_VTK_TETRA = 10
def __init__(self, filename: PathLike, library: str, mesh_id: int | None = None,
name: str = '', length_multiplier: float = 1.0,
@ -2655,7 +2654,8 @@ class UnstructuredMesh(MeshBase):
warnings.warn(
"The 'UnstructuredMesh.write_vtk_mesh' method has been renamed "
"to 'write_data_to_vtk' and will be removed in a future version "
" of OpenMC.", FutureWarning
" of OpenMC.",
FutureWarning,
)
self.write_data_to_vtk(**kwargs)
@ -2673,9 +2673,10 @@ class UnstructuredMesh(MeshBase):
Parameters
----------
filename : str or pathlib.Path
Name of the VTK file to write. If the filename ends in '.vtu' then a
binary VTU format file will be written, if the filename ends in
'.vtk' then a legacy VTK file will be written.
Name of the VTK file to write. If the filename ends in '.vtkhdf'
then a VTKHDF format file will be written. If the filename ends in
'.vtu' then a binary VTU format file will be written. If the
filename ends in '.vtk' then a legacy VTK file will be written.
datasets : dict
Dictionary whose keys are the data labels and values are numpy
appropriately sized arrays of the data
@ -2683,6 +2684,35 @@ class UnstructuredMesh(MeshBase):
Whether or not to normalize the data by the volume of the mesh
elements
"""
if Path(filename).suffix == ".vtkhdf":
self._write_data_to_vtk_hdf5_format(
filename=filename,
datasets=datasets,
volume_normalization=volume_normalization,
)
elif Path(filename).suffix == ".vtk" or Path(filename).suffix == ".vtu":
self._write_data_to_vtk_ascii_format(
filename=filename,
datasets=datasets,
volume_normalization=volume_normalization,
)
else:
raise ValueError(
"Unsupported file extension, The filename must end with "
"'.vtkhdf', '.vtu' or '.vtk'"
)
def _write_data_to_vtk_ascii_format(
self,
filename: PathLike | None = None,
datasets: dict | None = None,
volume_normalization: bool = True,
):
from vtkmodules.util import numpy_support
from vtkmodules import vtkCommonCore
from vtkmodules import vtkCommonDataModel
@ -2690,9 +2720,7 @@ class UnstructuredMesh(MeshBase):
from vtkmodules import vtkIOXML
if self.connectivity is None or self.vertices is None:
raise RuntimeError(
"This mesh has not been loaded from a statepoint file."
)
raise RuntimeError("This mesh has not been loaded from a statepoint file.")
if filename is None:
filename = f"mesh_{self.id}.vtk"
@ -2774,29 +2802,128 @@ class UnstructuredMesh(MeshBase):
writer.Write()
def _write_data_to_vtk_hdf5_format(
self,
filename: PathLike | None = None,
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
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]])
else:
# No -1 values, append the whole cell
trimmed_connectivity.append(cell)
trimmed_connectivity = np.array(trimmed_connectivity, dtype="int32").flatten()
# MOAB meshes supports tet elements only so we know it has 4 points per cell
points_per_cell = 4
# 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)
for name, data in datasets.items():
if data.shape != self.dimension:
raise ValueError(
f'Cannot apply dataset "{name}" with '
f"shape {data.shape} to mesh {self.id} "
f"with dimensions {self.dimension}"
)
with h5py.File(filename, "w") as f:
root = f.create_group("VTKHDF")
vtk_file_format_version = (2, 1)
root.attrs["Version"] = vtk_file_format_version
ascii_type = "UnstructuredGrid".encode("ascii")
root.attrs.create(
"Type",
ascii_type,
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")
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")
)
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)
@classmethod
def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
filename = group['filename'][()].decode()
library = group['library'][()].decode()
if 'options' in group.attrs:
filename = group["filename"][()].decode()
library = group["library"][()].decode()
if "options" in group.attrs:
options = group.attrs['options'].decode()
else:
options = None
mesh = cls(filename=filename, library=library, mesh_id=mesh_id, name=name, options=options)
mesh = cls(
filename=filename,
library=library,
mesh_id=mesh_id,
name=name,
options=options,
)
mesh._has_statepoint_data = True
vol_data = group['volumes'][()]
vol_data = group["volumes"][()]
mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],))
mesh.n_elements = mesh.volumes.size
vertices = group['vertices'][()]
vertices = group["vertices"][()]
mesh._vertices = vertices.reshape((-1, 3))
connectivity = group['connectivity'][()]
connectivity = group["connectivity"][()]
mesh._connectivity = connectivity.reshape((-1, 8))
mesh._element_types = group['element_types'][()]
mesh._element_types = group["element_types"][()]
if 'length_multiplier' in group:
mesh.length_multiplier = group['length_multiplier'][()]
if "length_multiplier" in group:
mesh.length_multiplier = group["length_multiplier"][()]
return mesh
@ -2815,7 +2942,7 @@ class UnstructuredMesh(MeshBase):
element.set("library", self._library)
if self.options is not None:
element.set('options', self.options)
element.set("options", self.options)
subelement = ET.SubElement(element, "filename")
subelement.text = str(self.filename)

View file

@ -13,8 +13,8 @@ GROUP_STRUCTURES = {}
- "XMAS-172_" designed for LWR analysis ([SAR1990]_, [SAN2004]_)
- "SHEM-361_" designed for LWR analysis to eliminate self-shielding calculations
of thermal resonances ([HFA2005]_, [SAN2007]_, [HEB2008]_)
- "SCALE-X" (where X is 44 which is designed for criticality analysis
and 252 is designed for thermal reactors) for the SCALE code suite
- "SCALE-X" (where X is 44 which is designed for criticality analysis, 252 is designed
for thermal reactors and 999 for multipurpose activation) for the SCALE code suite
([ZAL1999]_ and [REARDEN2013]_)
- "MPACT-X" (where X is 51 (PWR), 60 (BWR), 69 (Magnox)) from the MPACT_ reactor
physics code ([KIM2019]_ and [KIM2020]_)
@ -29,6 +29,7 @@ GROUP_STRUCTURES = {}
.. _SCALE44: https://www-nds.iaea.org/publications/indc/indc-czr-0001.pdf
.. _ECCO-33: https://serpent.vtt.fi/mediawiki/index.php/ECCO_33-group_structure
.. _SCALE252: https://oecd-nea.org/science/wpncs/amct/workingarea/meeting2013/EGAMCT2013_08.pdf
.. _SCALE999: https://info.ornl.gov/sites/publications/Files/Pub67728.pdf, https://www.nrc.gov/docs/ML1218/ML12184A002.pdf
.. _MPACT: https://vera.ornl.gov/mpact/
.. _XMAS-172: https://www-nds.iaea.org/wimsd/energy.htm
.. _SHEM-361: http://merlin.polymtl.ca/downloads/FP214.pdf
@ -593,7 +594,208 @@ GROUP_STRUCTURES['CCFE-709'] = np.array([
2.4000e8, 2.8000e8, 3.2000e8, 3.6000e8, 4.0000e8,
4.4000e8, 4.8000e8, 5.2000e8, 5.6000e8, 6.0000e8,
6.4000e8, 6.8000e8, 7.2000e8, 7.6000e8, 8.0000e8,
8.4000e8, 8.8000e8, 9.2000e8, 9.6000e8, 1.0000e9,])
8.4000e8, 8.8000e8, 9.2000e8, 9.6000e8, 1.0000e9])
GROUP_STRUCTURES['SCALE-999'] = np.array([
1.000e-5, 1.000e-4, 5.000e-4, 7.500e-4, 1.000e-3,
1.200e-3, 1.500e-3, 2.000e-3, 2.500e-3, 3.000e-3,
4.000e-3, 5.000e-3, 7.500e-3, 1.000e-2, 1.450e-2,
1.850e-2, 2.100e-2, 2.530e-2, 3.000e-2, 4.000e-2,
5.000e-2, 6.000e-2, 7.000e-2, 8.000e-2, 9.000e-2,
1.000e-1, 1.250e-1, 1.500e-1, 1.750e-1, 1.840e-1,
2.000e-1, 2.250e-1, 2.500e-1, 2.750e-1, 3.000e-1,
3.250e-1, 3.500e-1, 3.668e-1, 3.750e-1, 4.000e-1,
4.140e-1, 4.500e-1, 5.000e-1, 5.316e-1, 5.500e-1,
6.000e-1, 6.250e-1, 6.500e-1, 6.826e-1, 7.000e-1,
7.500e-1, 8.000e-1, 8.500e-1, 8.764e-1, 9.000e-1,
9.250e-1, 9.500e-1, 9.750e-1, 1.000e+0, 1.010e+0,
1.020e+0, 1.030e+0, 1.040e+0, 1.050e+0, 1.060e+0,
1.070e+0, 1.080e+0, 1.090e+0, 1.100e+0, 1.110e+0,
1.120e+0, 1.130e+0, 1.140e+0, 1.150e+0, 1.175e+0,
1.200e+0, 1.225e+0, 1.250e+0, 1.300e+0, 1.350e+0,
1.400e+0, 1.450e+0, 1.500e+0, 1.545e+0, 1.590e+0,
1.635e+0, 1.680e+0, 1.725e+0, 1.770e+0, 1.815e+0,
1.860e+0, 1.900e+0, 1.940e+0, 1.970e+0, 2.000e+0,
2.060e+0, 2.120e+0, 2.165e+0, 2.210e+0, 2.255e+0,
2.300e+0, 2.340e+0, 2.380e+0, 2.425e+0, 2.470e+0,
2.520e+0, 2.570e+0, 2.620e+0, 2.670e+0, 2.720e+0,
2.770e+0, 2.820e+0, 2.870e+0, 2.920e+0, 2.970e+0,
3.000e+0, 3.100e+0, 3.200e+0, 3.300e+0, 3.500e+0,
3.620e+0, 3.730e+0, 3.830e+0, 3.928e+0, 4.000e+0,
4.100e+0, 4.300e+0, 4.500e+0, 4.750e+0, 4.875e+0,
5.000e+0, 5.044e+0, 5.250e+0, 5.400e+0, 5.550e+0,
5.700e+0, 5.850e+0, 6.000e+0, 6.250e+0, 6.375e+0,
6.500e+0, 6.625e+0, 6.750e+0, 6.875e+0, 7.000e+0,
7.075e+0, 7.150e+0, 7.625e+0, 8.100e+0, 8.208e+0,
8.315e+0, 8.708e+0, 9.100e+0, 9.550e+0, 1.000e+1,
1.034e+1, 1.068e+1, 1.109e+1, 1.150e+1, 1.170e+1,
1.190e+1, 1.240e+1, 1.290e+1, 1.333e+1, 1.375e+1,
1.408e+1, 1.440e+1, 1.475e+1, 1.510e+1, 1.555e+1,
1.600e+1, 1.650e+1, 1.700e+1, 1.730e+1, 1.760e+1,
1.805e+1, 1.850e+1, 1.875e+1, 1.900e+1, 1.940e+1,
2.000e+1, 2.050e+1, 2.100e+1, 2.175e+1, 2.250e+1,
2.375e+1, 2.500e+1, 2.625e+1, 2.750e+1, 2.826e+1,
2.902e+1, 2.951e+1, 3.000e+1, 3.063e+1, 3.125e+1,
3.150e+1, 3.175e+1, 3.250e+1, 3.325e+1, 3.350e+1,
3.375e+1, 3.418e+1, 3.460e+1, 3.500e+1, 3.550e+1,
3.600e+1, 3.700e+1, 3.713e+1, 3.727e+1, 3.763e+1,
3.800e+1, 3.855e+1, 3.910e+1, 3.935e+1, 3.960e+1,
4.030e+1, 4.100e+1, 4.170e+1, 4.240e+1, 4.320e+1,
4.400e+1, 4.460e+1, 4.520e+1, 4.610e+1, 4.700e+1,
4.743e+1, 4.785e+1, 4.808e+1, 4.830e+1, 4.875e+1,
4.920e+1, 4.990e+1, 5.060e+1, 5.130e+1, 5.200e+1,
5.270e+1, 5.340e+1, 5.620e+1, 5.800e+1, 6.000e+1,
6.100e+1, 6.122e+1, 6.144e+1, 6.300e+1, 6.500e+1,
6.625e+1, 6.750e+1, 6.975e+1, 7.200e+1, 7.400e+1,
7.600e+1, 7.745e+1, 7.889e+1, 7.945e+1, 8.000e+1,
8.170e+1, 8.200e+1, 8.400e+1, 8.600e+1, 8.800e+1,
9.000e+1, 9.250e+1, 9.500e+1, 9.700e+1, 1.000e+2,
1.012e+2, 1.038e+2, 1.050e+2, 1.080e+2, 1.105e+2,
1.130e+2, 1.160e+2, 1.175e+2, 1.190e+2, 1.205e+2,
1.220e+2, 1.240e+2, 1.260e+2, 1.280e+2, 1.301e+2,
1.325e+2, 1.350e+2, 1.375e+2, 1.400e+2, 1.430e+2,
1.450e+2, 1.475e+2, 1.500e+2, 1.525e+2, 1.550e+2,
1.575e+2, 1.600e+2, 1.625e+2, 1.650e+2, 1.670e+2,
1.700e+2, 1.716e+2, 1.739e+2, 1.763e+2, 1.786e+2,
1.800e+2, 1.877e+2, 1.885e+2, 1.915e+2, 1.930e+2,
1.962e+2, 1.999e+2, 2.020e+2, 2.074e+2, 2.088e+2,
2.095e+2, 2.122e+2, 2.145e+2, 2.175e+2, 2.200e+2,
2.237e+2, 2.269e+2, 2.301e+2, 2.333e+2, 2.367e+2,
2.400e+2, 2.442e+2, 2.484e+2, 2.527e+2, 2.571e+2,
2.615e+2, 2.661e+2, 2.707e+2, 2.754e+2, 2.801e+2,
2.850e+2, 2.899e+2, 2.948e+2, 2.999e+2, 3.050e+2,
3.107e+2, 3.165e+2, 3.224e+2, 3.284e+2, 3.345e+2,
3.408e+2, 3.471e+2, 3.536e+2, 3.591e+2, 3.648e+2,
3.705e+2, 3.764e+2, 3.823e+2, 3.883e+2, 3.944e+2,
4.007e+2, 4.070e+2, 4.134e+2, 4.199e+2, 4.265e+2,
4.332e+2, 4.400e+2, 4.470e+2, 4.540e+2, 4.595e+2,
4.650e+2, 4.706e+2, 4.763e+2, 4.821e+2, 4.879e+2,
4.937e+2, 4.997e+2, 5.057e+2, 5.118e+2, 5.180e+2,
5.243e+2, 5.306e+2, 5.370e+2, 5.435e+2, 5.500e+2,
5.581e+2, 5.662e+2, 5.745e+2, 5.830e+2, 5.932e+2,
6.036e+2, 6.142e+2, 6.250e+2, 6.359e+2, 6.471e+2,
6.585e+2, 6.700e+2, 6.765e+2, 6.830e+2, 6.909e+2,
6.988e+2, 7.069e+2, 7.150e+2, 7.232e+2, 7.316e+2,
7.400e+2, 7.485e+2, 7.598e+2, 7.712e+2, 7.827e+2,
7.945e+2, 8.064e+2, 8.185e+2, 8.308e+2, 8.433e+2,
8.559e+2, 8.688e+2, 8.818e+2, 8.950e+2, 9.085e+2,
9.221e+2, 9.360e+2, 9.500e+2, 9.555e+2, 9.611e+2,
9.720e+2, 9.829e+2, 9.940e+2, 1.005e+3, 1.017e+3,
1.028e+3, 1.040e+3, 1.051e+3, 1.063e+3, 1.075e+3,
1.087e+3, 1.100e+3, 1.112e+3, 1.125e+3, 1.137e+3,
1.150e+3, 1.171e+3, 1.191e+3, 1.213e+3, 1.234e+3,
1.249e+3, 1.265e+3, 1.280e+3, 1.296e+3, 1.312e+3,
1.328e+3, 1.344e+3, 1.361e+3, 1.377e+3, 1.394e+3,
1.411e+3, 1.429e+3, 1.446e+3, 1.464e+3, 1.482e+3,
1.500e+3, 1.525e+3, 1.550e+3, 1.585e+3, 1.610e+3,
1.636e+3, 1.662e+3, 1.689e+3, 1.716e+3, 1.744e+3,
1.772e+3, 1.800e+3, 1.828e+3, 1.856e+3, 1.885e+3,
1.914e+3, 1.943e+3, 1.973e+3, 2.004e+3, 2.035e+3,
2.075e+3, 2.116e+3, 2.158e+3, 2.200e+3, 2.250e+3,
2.290e+3, 2.337e+3, 2.386e+3, 2.435e+3, 2.500e+3,
2.532e+3, 2.580e+3, 2.613e+3, 2.679e+3, 2.747e+3,
2.808e+3, 2.871e+3, 2.935e+3, 3.000e+3, 3.035e+3,
3.112e+3, 3.191e+3, 3.272e+3, 3.355e+3, 3.440e+3,
3.527e+3, 3.616e+3, 3.707e+3, 3.740e+3, 3.819e+3,
3.900e+3, 3.998e+3, 4.099e+3, 4.202e+3, 4.307e+3,
4.400e+3, 4.500e+3, 4.620e+3, 4.740e+3, 4.850e+3,
4.960e+3, 5.100e+3, 5.250e+3, 5.400e+3, 5.531e+3,
5.645e+3, 5.700e+3, 5.879e+3, 6.000e+3, 6.128e+3,
6.258e+3, 6.392e+3, 6.528e+3, 6.667e+3, 6.809e+3,
6.954e+3, 7.102e+3, 7.212e+3, 7.323e+3, 7.437e+3,
7.552e+3, 7.669e+3, 7.787e+3, 7.908e+3, 8.030e+3,
8.159e+3, 8.289e+3, 8.422e+3, 8.557e+3, 8.694e+3,
8.834e+3, 8.975e+3, 9.119e+3, 9.307e+3, 9.500e+3,
9.763e+3, 1.003e+4, 1.031e+4, 1.060e+4, 1.086e+4,
1.114e+4, 1.142e+4, 1.171e+4, 1.202e+4, 1.234e+4,
1.266e+4, 1.300e+4, 1.324e+4, 1.348e+4, 1.373e+4,
1.398e+4, 1.424e+4, 1.450e+4, 1.476e+4, 1.503e+4,
1.527e+4, 1.550e+4, 1.574e+4, 1.599e+4, 1.623e+4,
1.649e+4, 1.674e+4, 1.700e+4, 1.727e+4, 1.755e+4,
1.783e+4, 1.812e+4, 1.841e+4, 1.870e+4, 1.900e+4,
1.931e+4, 1.965e+4, 2.000e+4, 2.045e+4, 2.092e+4,
2.139e+4, 2.188e+4, 2.229e+4, 2.271e+4, 2.314e+4,
2.358e+4, 2.388e+4, 2.418e+4, 2.450e+4, 2.479e+4,
2.500e+4, 2.520e+4, 2.552e+4, 2.580e+4, 2.606e+4,
2.653e+4, 2.700e+4, 2.737e+4, 2.774e+4, 2.812e+4,
2.850e+4, 2.887e+4, 2.924e+4, 2.962e+4, 3.000e+4,
3.045e+4, 3.090e+4, 3.136e+4, 3.183e+4, 3.243e+4,
3.304e+4, 3.367e+4, 3.431e+4, 3.468e+4, 3.507e+4,
3.545e+4, 3.584e+4, 3.624e+4, 3.663e+4, 3.704e+4,
3.744e+4, 3.786e+4, 3.827e+4, 3.869e+4, 3.912e+4,
3.955e+4, 3.998e+4, 4.042e+4, 4.087e+4, 4.136e+4,
4.186e+4, 4.237e+4, 4.288e+4, 4.340e+4, 4.393e+4,
4.446e+4, 4.500e+4, 4.565e+4, 4.631e+4, 4.721e+4,
4.812e+4, 4.905e+4, 5.000e+4, 5.099e+4, 5.200e+4,
5.248e+4, 5.347e+4, 5.448e+4, 5.551e+4, 5.656e+4,
5.740e+4, 5.826e+4, 5.912e+4, 6.000e+4, 6.088e+4,
6.177e+4, 6.267e+4, 6.358e+4, 6.451e+4, 6.545e+4,
6.641e+4, 6.738e+4, 6.851e+4, 6.965e+4, 7.081e+4,
7.200e+4, 7.300e+4, 7.399e+4, 7.500e+4, 7.610e+4,
7.722e+4, 7.835e+4, 7.950e+4, 8.074e+4, 8.200e+4,
8.250e+4, 8.374e+4, 8.500e+4, 8.652e+4, 8.788e+4,
8.926e+4, 9.067e+4, 9.210e+4, 9.355e+4, 9.502e+4,
9.652e+4, 9.804e+4, 1.000e+5, 1.013e+5, 1.027e+5,
1.040e+5, 1.054e+5, 1.068e+5, 1.082e+5, 1.096e+5,
1.111e+5, 1.125e+5, 1.139e+5, 1.153e+5, 1.168e+5,
1.183e+5, 1.197e+5, 1.213e+5, 1.228e+5, 1.241e+5,
1.255e+5, 1.269e+5, 1.283e+5, 1.291e+5, 1.307e+5,
1.323e+5, 1.340e+5, 1.357e+5, 1.374e+5, 1.391e+5,
1.409e+5, 1.426e+5, 1.445e+5, 1.463e+5, 1.481e+5,
1.490e+5, 1.519e+5, 1.538e+5, 1.557e+5, 1.576e+5,
1.596e+5, 1.616e+5, 1.637e+5, 1.657e+5, 1.678e+5,
1.699e+5, 1.721e+5, 1.742e+5, 1.764e+5, 1.786e+5,
1.809e+5, 1.832e+5, 1.855e+5, 1.878e+5, 1.902e+5,
1.926e+5, 1.962e+5, 2.000e+5, 2.024e+5, 2.050e+5,
2.076e+5, 2.102e+5, 2.128e+5, 2.155e+5, 2.182e+5,
2.209e+5, 2.237e+5, 2.265e+5, 2.294e+5, 2.323e+5,
2.352e+5, 2.381e+5, 2.411e+5, 2.442e+5, 2.472e+5,
2.500e+5, 2.527e+5, 2.555e+5, 2.584e+5, 2.612e+5,
2.641e+5, 2.670e+5, 2.700e+5, 2.732e+5, 2.767e+5,
2.802e+5, 2.837e+5, 2.873e+5, 2.909e+5, 2.945e+5,
2.972e+5, 2.985e+5, 3.020e+5, 3.053e+5, 3.088e+5,
3.122e+5, 3.157e+5, 3.192e+5, 3.228e+5, 3.264e+5,
3.300e+5, 3.337e+5, 3.379e+5, 3.422e+5, 3.465e+5,
3.508e+5, 3.553e+5, 3.597e+5, 3.643e+5, 3.688e+5,
3.735e+5, 3.782e+5, 3.829e+5, 3.877e+5, 3.938e+5,
4.000e+5, 4.076e+5, 4.138e+5, 4.200e+5, 4.249e+5,
4.299e+5, 4.349e+5, 4.400e+5, 4.452e+5, 4.505e+5,
4.553e+5, 4.601e+5, 4.650e+5, 4.700e+5, 4.772e+5,
4.845e+5, 4.920e+5, 4.995e+5, 5.054e+5, 5.113e+5,
5.173e+5, 5.234e+5, 5.299e+5, 5.365e+5, 5.432e+5,
5.500e+5, 5.557e+5, 5.614e+5, 5.672e+5, 5.730e+5,
5.784e+5, 5.891e+5, 6.000e+5, 6.081e+5, 6.158e+5,
6.235e+5, 6.313e+5, 6.393e+5, 6.468e+5, 6.545e+5,
6.622e+5, 6.700e+5, 6.790e+5, 6.926e+5, 7.065e+5,
7.154e+5, 7.244e+5, 7.335e+5, 7.427e+5, 7.500e+5,
7.576e+5, 7.653e+5, 7.730e+5, 7.808e+5, 7.904e+5,
8.002e+5, 8.100e+5, 8.200e+5, 8.301e+5, 8.403e+5,
8.506e+5, 8.611e+5, 8.750e+5, 8.874e+5, 9.000e+5,
9.072e+5, 9.200e+5, 9.400e+5, 9.616e+5, 9.800e+5,
1.003e+6, 1.010e+6, 1.040e+6, 1.070e+6, 1.100e+6,
1.108e+6, 1.136e+6, 1.165e+6, 1.200e+6, 1.225e+6,
1.250e+6, 1.287e+6, 1.317e+6, 1.356e+6, 1.400e+6,
1.423e+6, 1.461e+6, 1.500e+6, 1.536e+6, 1.572e+6,
1.612e+6, 1.653e+6, 1.695e+6, 1.738e+6, 1.782e+6,
1.827e+6, 1.850e+6, 1.921e+6, 1.969e+6, 2.019e+6,
2.070e+6, 2.123e+6, 2.176e+6, 2.231e+6, 2.307e+6,
2.354e+6, 2.365e+6, 2.385e+6, 2.466e+6, 2.479e+6,
2.535e+6, 2.592e+6, 2.658e+6, 2.725e+6, 2.794e+6,
2.865e+6, 2.932e+6, 3.000e+6, 3.080e+6, 3.166e+6,
3.247e+6, 3.329e+6, 3.413e+6, 3.499e+6, 3.588e+6,
3.679e+6, 3.772e+6, 3.867e+6, 3.965e+6, 4.066e+6,
4.183e+6, 4.304e+6, 4.398e+6, 4.493e+6, 4.607e+6,
4.724e+6, 4.800e+6, 4.882e+6, 4.966e+6, 5.092e+6,
5.221e+6, 5.353e+6, 5.488e+6, 5.627e+6, 5.770e+6,
5.916e+6, 6.065e+6, 6.219e+6, 6.376e+6, 6.434e+6,
6.592e+6, 6.703e+6, 6.873e+6, 7.047e+6, 7.225e+6,
7.408e+6, 7.596e+6, 7.788e+6, 7.985e+6, 8.187e+6,
8.395e+6, 8.607e+6, 8.825e+6, 9.048e+6, 9.278e+6,
9.512e+6, 9.753e+6, 1.000e+7, 1.025e+7, 1.051e+7,
1.078e+7, 1.105e+7, 1.133e+7, 1.162e+7, 1.191e+7,
1.221e+7, 1.252e+7, 1.284e+7, 1.317e+7, 1.350e+7,
1.384e+7, 1.419e+7, 1.455e+7, 1.492e+7, 1.530e+7,
1.568e+7, 1.608e+7, 1.649e+7, 1.691e+7, 1.733e+7,
1.790e+7, 1.845e+7, 1.900e+7, 1.964e+7, 2.000e+7])
GROUP_STRUCTURES['UKAEA-1102'] = np.array([
1.0000e-5, 1.0471e-5, 1.0965e-5, 1.1482e-5, 1.2023e-5,
1.2589e-5, 1.3183e-5, 1.3804e-5, 1.4454e-5, 1.5136e-5,

View file

@ -10,6 +10,7 @@ import numpy as np
import openmc
import openmc.mgxs
import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from ..tallies import ESTIMATOR_TYPES
@ -555,14 +556,14 @@ class Library:
self.all_mgxs[domain.id][mgxs_type] = mgxs
def add_to_tallies_file(self, tallies_file, merge=True):
"""Add all tallies from all MGXS objects to a tallies file.
def add_to_tallies(self, tallies, merge=True):
"""Add tallies from all MGXS objects to a tallies object.
NOTE: This assumes that :meth:`Library.build_library` has been called
Parameters
----------
tallies_file : openmc.Tallies
tallies : openmc.Tallies
A Tallies collection to add each MGXS' tallies to generate a
'tallies.xml' input file for OpenMC
merge : bool
@ -571,7 +572,7 @@ class Library:
"""
cv.check_type('tallies_file', tallies_file, openmc.Tallies)
cv.check_type('tallies', tallies, openmc.Tallies)
# Add tallies from each MGXS for each domain and mgxs type
for domain in self.domains:
@ -586,7 +587,15 @@ class Library:
= list(range(1, self.num_delayed_groups + 1))
for tally in mgxs.tallies.values():
tallies_file.append(tally, merge=merge)
tallies.append(tally, merge=merge)
def add_to_tallies_file(self, tallies_file, merge=True):
warn(
"The Library.add_to_tallies_file(...) method has been renamed to"
"add_to_tallies(...) and will be removed in a future version of "
"OpenMC.", FutureWarning
)
self.add_to_tallies(tallies_file, merge=merge)
def load_from_statepoint(self, statepoint):
"""Extracts tallies in an OpenMC StatePoint with the data needed to
@ -851,7 +860,7 @@ class Library:
'since a statepoint has not yet been loaded'
raise ValueError(msg)
cv.check_type('filename', filename, str)
cv.check_type('filename', filename, (str, PathLike))
cv.check_type('directory', directory, str)
import h5py
@ -894,7 +903,7 @@ class Library:
"""
cv.check_type('filename', filename, str)
cv.check_type('filename', filename, (str, PathLike))
cv.check_type('directory', directory, str)
# Make directory if it does not exist
@ -930,7 +939,7 @@ class Library:
"""
cv.check_type('filename', filename, str)
cv.check_type('filename', filename, (str, PathLike))
cv.check_type('directory', directory, str)
# Make directory if it does not exist

View file

@ -7,6 +7,7 @@ import numpy as np
import openmc
import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from openmc.mgxs import MGXS
from .mgxs import _DOMAIN_TO_FILTER
@ -722,7 +723,7 @@ class MDGXS(MGXS):
"""
cv.check_type('filename', filename, str)
cv.check_type('filename', filename, (str, PathLike))
cv.check_type('directory', directory, str)
cv.check_value('format', format, ['csv', 'excel', 'pickle', 'latex'])
cv.check_value('xs_type', xs_type, ['macro', 'micro'])

View file

@ -2,6 +2,7 @@ import copy
from numbers import Integral
import os
import warnings
from textwrap import dedent
import h5py
import numpy as np
@ -9,6 +10,7 @@ import numpy as np
import openmc
from openmc.data import REACTION_MT, REACTION_NAME, FISSION_MTS
import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from ..tallies import ESTIMATOR_TYPES
from . import EnergyGroups
@ -164,7 +166,7 @@ class MGXS:
"""
_params = """
_params = dedent("""
Parameters
----------
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh
@ -251,7 +253,7 @@ class MGXS:
.. versionadded:: 0.13.1
"""
""")
# Store whether or not the number density should be removed for microscopic
# values of this data
@ -1981,7 +1983,7 @@ class MGXS:
"""
cv.check_type('filename', filename, str)
cv.check_type('filename', filename, (str, PathLike))
cv.check_type('directory', directory, str)
cv.check_value('format', format, ['csv', 'excel', 'pickle', 'latex'])
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
@ -2125,8 +2127,8 @@ class MGXS:
df['std. dev.'] /= np.tile(densities, tile_factor)
# Replace NaNs by zeros (happens if nuclide density is zero)
df['mean'].replace(np.nan, 0.0, inplace=True)
df['std. dev.'].replace(np.nan, 0.0, inplace=True)
df['mean'] = df['mean'].replace(np.nan, 0.0)
df['std. dev.'] = df['std. dev.'].replace(np.nan, 0.0)
# Sort the dataframe by domain type id (e.g., distribcell id) and
# energy groups such that data is from fast to thermal

View file

@ -11,7 +11,7 @@ import openmc
import openmc.mgxs
from openmc.mgxs import SCATTER_TABULAR, SCATTER_LEGENDRE, SCATTER_HISTOGRAM
from .checkvalue import check_type, check_value, check_greater_than, \
check_iterable_type, check_less_than, check_filetype_version
check_iterable_type, check_less_than, check_filetype_version, PathLike
ROOM_TEMPERATURE_KELVIN = 294.0
@ -2506,7 +2506,7 @@ class MGXSLibrary:
Parameters
----------
filename : str
filename : str or PathLike
Filename of file, default is mgxs.h5.
libver : {'earliest', 'latest'}
Compatibility mode for the HDF5 file. 'latest' will produce files
@ -2514,8 +2514,7 @@ class MGXSLibrary:
"""
check_type('filename', filename, str)
check_type('filename', filename, (str, PathLike))
# Create and write to the HDF5 file
file = h5py.File(filename, "w", libver=libver)
file.attrs['filetype'] = np.bytes_(_FILETYPE_MGXS_LIBRARY)
@ -2554,7 +2553,7 @@ class MGXSLibrary:
raise ValueError("Either path or openmc.config['mg_cross_sections']"
"must be set")
check_type('filename', filename, str)
check_type('filename', filename, (str, PathLike))
file = h5py.File(filename, 'r')
# Check filetype and version

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