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

Merge branch 'develop' of https://github.com/openmc-dev/openmc into openmc-dev-develop

Browse source
This commit is contained in:
church89 2026-01-13 13:30:50 +01:00
commit 8d7ac63151
838 changed files with 43374 additions and 21920 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -75,12 +75,18 @@ 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.
RDMAV_FORK_SAFE: 1
steps:
- name: Setup cmake
uses: jwlawson/actions-setup-cmake@v2
with:
cmake-version: '3.31'
- name: Checkout repository
uses: actions/checkout@v4
with:
@ -148,7 +154,7 @@ jobs:
path: |
~/nndc_hdf5
~/endf-b-vii.1
key: ${{ runner.os }}-build-xs-cache
key: ${{ runner.os }}-build-xs-cache-${{ hashFiles(format('{0}/tools/ci/download-xs.sh', github.workspace)) }}
- name: before
shell: bash
@ -166,11 +172,38 @@ 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 \
--merge-mode-functions=separate \
--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
@ -179,5 +212,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

68
CITATION.cff Normal file
View file

@ -0,0 +1,68 @@
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"
- family-names: Horelik
given-names: Nicholas E.
- family-names: Herman
given-names: Bryan R.
- family-names: Nelson
given-names: Adam G.
orcid: "https://orcid.org/0000-0002-3614-0676"
- family-names: Forget
given-names: Benoit
orcid: "https://orcid.org/0000-0003-1459-7672"
- family-names: Smith
given-names: Kord
contact:
- family-names: Romano
given-names: Paul K.
orcid: "https://orcid.org/0000-0002-1147-045X"
doi: 10.1016/j.anucene.2014.07.048
issn: 0306-4549
volume: 82
journal: Annals of Nuclear Energy
publisher:
name: Elsevier
start: 90
end: 97
year: 2015
month: 8
title: "OpenMC: A state-of-the-art Monte Carlo code for research and development"
type: article
url: "https://doi.org/10.1016/j.anucene.2014.07.048"

64
CMakeLists.txt
View file

@ -36,8 +36,8 @@ option(OPENMC_ENABLE_COVERAGE "Compile with coverage analysis flags"
option(OPENMC_USE_DAGMC "Enable support for DAGMC (CAD) geometry" OFF)
option(OPENMC_USE_LIBMESH "Enable support for libMesh unstructured mesh tallies" OFF)
option(OPENMC_USE_MPI "Enable MPI" OFF)
option(OPENMC_USE_MCPL "Enable MCPL" OFF)
option(OPENMC_USE_UWUW "Enable UWUW" OFF)
option(OPENMC_FORCE_VENDORED_LIBS "Explicitly use submodules defined in 'vendor'" OFF)
message(STATUS "OPENMC_USE_OPENMP ${OPENMC_USE_OPENMP}")
message(STATUS "OPENMC_BUILD_TESTS ${OPENMC_BUILD_TESTS}")
@ -46,8 +46,8 @@ message(STATUS "OPENMC_ENABLE_COVERAGE ${OPENMC_ENABLE_COVERAGE}")
message(STATUS "OPENMC_USE_DAGMC ${OPENMC_USE_DAGMC}")
message(STATUS "OPENMC_USE_LIBMESH ${OPENMC_USE_LIBMESH}")
message(STATUS "OPENMC_USE_MPI ${OPENMC_USE_MPI}")
message(STATUS "OPENMC_USE_MCPL ${OPENMC_USE_MCPL}")
message(STATUS "OPENMC_USE_UWUW ${OPENMC_USE_UWUW}")
message(STATUS "OPENMC_FORCE_VENDORED_LIBS ${OPENMC_FORCE_VENDORED_LIBS}")
# Warnings for deprecated options
foreach(OLD_OPT IN ITEMS "openmp" "profile" "coverage" "dagmc" "libmesh")
@ -189,15 +189,6 @@ if(${HDF5_VERSION} VERSION_GREATER_EQUAL 1.12.0)
list(APPEND cxxflags -DH5Oget_info_by_idx_vers=1 -DH5O_info_t_vers=1)
endif()
#===============================================================================
# MCPL
#===============================================================================
if (OPENMC_USE_MCPL)
find_package(MCPL REQUIRED)
message(STATUS "Found MCPL: ${MCPL_DIR} (found version \"${MCPL_VERSION}\")")
endif()
#===============================================================================
# Set compile/link flags based on which compiler is being used
#===============================================================================
@ -249,31 +240,47 @@ endif()
# pugixml library
#===============================================================================
find_package_write_status(pugixml)
if (NOT pugixml_FOUND)
if(OPENMC_FORCE_VENDORED_LIBS)
add_subdirectory(vendor/pugixml)
set_target_properties(pugixml PROPERTIES CXX_STANDARD 14 CXX_EXTENSIONS OFF)
else()
find_package_write_status(pugixml)
if (NOT pugixml_FOUND)
add_subdirectory(vendor/pugixml)
set_target_properties(pugixml PROPERTIES CXX_STANDARD 14 CXX_EXTENSIONS OFF)
endif()
endif()
#===============================================================================
# {fmt} library
#===============================================================================
find_package_write_status(fmt)
if (NOT fmt_FOUND)
if(OPENMC_FORCE_VENDORED_LIBS)
set(FMT_INSTALL ON CACHE BOOL "Generate the install target.")
add_subdirectory(vendor/fmt)
else()
find_package_write_status(fmt)
if (NOT fmt_FOUND)
set(FMT_INSTALL ON CACHE BOOL "Generate the install target.")
add_subdirectory(vendor/fmt)
endif()
endif()
#===============================================================================
# xtensor header-only library
#===============================================================================
find_package_write_status(xtensor)
if (NOT xtensor_FOUND)
if(OPENMC_FORCE_VENDORED_LIBS)
add_subdirectory(vendor/xtl)
set(xtl_DIR ${CMAKE_CURRENT_BINARY_DIR}/vendor/xtl)
add_subdirectory(vendor/xtensor)
else()
find_package_write_status(xtensor)
if (NOT xtensor_FOUND)
add_subdirectory(vendor/xtl)
set(xtl_DIR ${CMAKE_CURRENT_BINARY_DIR}/vendor/xtl)
add_subdirectory(vendor/xtensor)
endif()
endif()
#===============================================================================
@ -281,9 +288,13 @@ endif()
#===============================================================================
if(OPENMC_BUILD_TESTS)
find_package_write_status(Catch2)
if (NOT Catch2_FOUND)
if (OPENMC_FORCE_VENDORED_LIBS)
add_subdirectory(vendor/Catch2)
else()
find_package_write_status(Catch2)
if (NOT Catch2_FOUND)
add_subdirectory(vendor/Catch2)
endif()
endif()
endif()
@ -327,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
@ -343,6 +355,7 @@ list(APPEND libopenmc_SOURCES
src/geometry.cpp
src/geometry_aux.cpp
src/hdf5_interface.cpp
src/ifp.cpp
src/initialize.cpp
src/lattice.cpp
src/material.cpp
@ -406,6 +419,7 @@ list(APPEND libopenmc_SOURCES
src/tallies/filter_materialfrom.cpp
src/tallies/filter_mesh.cpp
src/tallies/filter_meshborn.cpp
src/tallies/filter_meshmaterial.cpp
src/tallies/filter_meshsurface.cpp
src/tallies/filter_mu.cpp
src/tallies/filter_musurface.cpp
@ -417,6 +431,7 @@ list(APPEND libopenmc_SOURCES
src/tallies/filter_surface.cpp
src/tallies/filter_time.cpp
src/tallies/filter_universe.cpp
src/tallies/filter_weight.cpp
src/tallies/filter_zernike.cpp
src/tallies/tally.cpp
src/tallies/tally_scoring.cpp
@ -489,11 +504,11 @@ else()
endif()
if(OPENMC_USE_DAGMC)
target_compile_definitions(libopenmc PRIVATE DAGMC)
target_compile_definitions(libopenmc PRIVATE OPENMC_DAGMC_ENABLED)
target_link_libraries(libopenmc dagmc-shared)
if(OPENMC_USE_UWUW)
target_compile_definitions(libopenmc PRIVATE OPENMC_UWUW)
target_compile_definitions(libopenmc PRIVATE OPENMC_UWUW_ENABLED)
target_link_libraries(libopenmc uwuw-shared)
endif()
elseif(OPENMC_USE_UWUW)
@ -502,7 +517,7 @@ elseif(OPENMC_USE_UWUW)
endif()
if(OPENMC_USE_LIBMESH)
target_compile_definitions(libopenmc PRIVATE LIBMESH)
target_compile_definitions(libopenmc PRIVATE OPENMC_LIBMESH_ENABLED)
target_link_libraries(libopenmc PkgConfig::LIBMESH)
endif()
@ -525,11 +540,6 @@ if (OPENMC_BUILD_TESTS)
add_subdirectory(tests/cpp_unit_tests)
endif()
if (OPENMC_USE_MCPL)
target_compile_definitions(libopenmc PUBLIC OPENMC_MCPL)
target_link_libraries(libopenmc MCPL::mcpl)
endif()
#===============================================================================
# Log build info that this executable can report later
#===============================================================================

3
Dockerfile
View file

@ -24,7 +24,7 @@ ARG compile_cores=1
ARG build_dagmc=off
ARG build_libmesh=off
FROM debian:bookworm-slim AS dependencies
FROM ubuntu:24.04 AS dependencies
ARG compile_cores
ARG build_dagmc
@ -96,6 +96,7 @@ RUN if [ "$build_dagmc" = "on" ]; then \
# Install addition packages required for DAGMC
apt-get -y install libeigen3-dev libnetcdf-dev libtbb-dev libglfw3-dev \
&& pip install --upgrade numpy \
&& pip install --no-cache-dir setuptools cython \
# Clone and install EMBREE
&& mkdir -p $HOME/EMBREE && cd $HOME/EMBREE \
&& git clone --single-branch -b ${EMBREE_TAG} --depth 1 ${EMBREE_REPO} \

10
cmake/Modules/GetVersionFromGit.cmake
View file

@ -37,11 +37,17 @@ if(EXISTS "${CMAKE_SOURCE_DIR}/.git" AND GIT_FOUND)
WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
OUTPUT_VARIABLE VERSION_STRING
OUTPUT_STRIP_TRAILING_WHITESPACE
ERROR_QUIET
)
# If no tags are found, instruct user to fetch them
# If no tags are found, set version to 0 and show a warning
if(VERSION_STRING STREQUAL "")
message(FATAL_ERROR "No git tags found. Run 'git fetch --tags' and try again.")
set(VERSION_STRING "0.0.0")
message(WARNING
"No git tags found. Version set to 0.0.0.\n"
"Run 'git fetch --tags' to ensure proper versioning.\n"
"For more information, see OpenMC developer documentation."
)
endif()
# Extract the commit hash

15
cmake/OpenMCConfig.cmake.in
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()
@ -29,10 +32,6 @@ if(@OPENMC_USE_OPENMP@)
find_package(OpenMP REQUIRED)
endif()
if(@OPENMC_USE_MCPL@)
find_package(MCPL REQUIRED)
endif()
if(@OPENMC_USE_UWUW@ AND NOT ${DAGMC_BUILD_UWUW})
message(FATAL_ERROR "UWUW is enabled in OpenMC but the DAGMC installation discovered was not configured with UWUW.")
endif()

BIN
docs/source/_images/sphere-mesh-vtk.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 72 KiB

29
docs/source/capi/index.rst
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

2
docs/source/conf.py
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'

8
docs/source/devguide/policies.rst
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.

24
docs/source/devguide/workflow.rst
View file

@ -91,6 +91,30 @@ features and bug fixes. The general steps for contributing are as follows:
6. After the pull request has been thoroughly vetted, it is merged back into the
*develop* branch of openmc-dev/openmc.
Setting Up Upstream Tracking (Required for Versioning)
------------------------------------------------------
By default, your fork **does not** include tags from the upstream OpenMC repository.
OpenMC relies on `git describe --tags` for versioning in source builds, and missing tags can lead
to incorrect version detection (i.e., ``0.0.0``). To ensure proper versioning, follow these steps:
1. **Add the Upstream Repository**
This allows you to fetch updates from the main OpenMC repository.
.. code-block:: sh
git remote add upstream https://github.com/openmc-dev/openmc.git
2. **Fetch and Push Tags**
Retrieve tags from the upstream repository and update your fork:
.. code-block:: sh
git fetch --tags upstream
git push --tags origin
This ensures that both your **local** and **remote** fork have the correct versioning information.
Private Development
-------------------

46
docs/source/io_formats/collision_track.rst Normal file
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.

21
docs/source/io_formats/depletion_chain.rst
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
----------------------

26
docs/source/io_formats/depletion_results.rst
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.

15
docs/source/io_formats/geometry.rst
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"
@ -318,9 +316,10 @@ the following attributes or sub-elements:
*Default*: None
:orientation:
The orientation of the hexagonal lattice. The string "x" indicates that two
sides of the lattice are parallel to the x-axis, whereas the string "y"
indicates that two sides are parallel to the y-axis.
The orientation of the hexagonal lattice. The string "x" indicates that each
lattice element has two faces that are perpendicular to the x-axis, whereas
the string "y" indicates that each lattice element has two faces that are
perpendicular to the y-axis.
*Default*: "y"

1
docs/source/io_formats/index.rst
View file

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

3
docs/source/io_formats/properties.rst
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/**

141
docs/source/io_formats/settings.rst
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
-----------------------------------
@ -488,6 +577,14 @@ found in the :ref:`random ray user guide <random_ray>`.
:type:
The type of the domain. Can be ``material``, ``cell``, or ``universe``.
:diagonal_stabilization_rho:
The rho factor for use with diagonal stabilization. This technique is
applied when negative diagonal (in-group) elements are detected in
the scattering matrix of input MGXS data, which is a common feature
of transport corrected MGXS data.
*Default*: 1.0
----------------------------------
``<resonance_scattering>`` Element
----------------------------------
@ -747,13 +844,18 @@ attributes/sub-elements:
relative source strength of each mesh element or each point in the cloud.
:volume_normalized:
For "mesh" spatial distrubtions, this optional boolean element specifies
For "mesh" spatial distributions, this optional boolean element specifies
whether the vector of relative strengths should be multiplied by the mesh
element volume. This is most common if the strengths represent a source
per unit volume.
*Default*: false
:bias:
For "mesh" and "cloud" spatial distributions, this optional element
specifies floating point values corresponding to alternative probabilities
for each value/component to use for biased sampling.
:angle:
An element specifying the angular distribution of source sites. This element
has the following attributes:
@ -786,6 +888,10 @@ attributes/sub-elements:
are those of a univariate probability distribution (see the description in
:ref:`univariate`).
:bias:
For "isotropic" angular distributions, this optional element specifies a
"mu-phi" angular distribution used for biased sampling.
:energy:
An element specifying the energy distribution of source sites. The necessary
sub-elements/attributes are those of a univariate probability distribution
@ -809,6 +915,10 @@ attributes/sub-elements:
mesh element and follows the format for :ref:`source_element`. The number of
``<source>`` sub-elements should correspond to the number of mesh elements.
.. note:: Biased sampling can be applied to the spatial and energy distributions
of a source by using the ``<bias>`` sub-element (see
:ref:`univariate` for details on how to specify bias distributions).
:constraints:
This sub-element indicates the presence of constraints on sampled source
sites (see :ref:`usersguide_source_constraints` for details). It may have
@ -901,13 +1011,36 @@ variable and whose sub-elements/attributes are as follows:
*Default*: histogram
:pair:
For a "mixture" distribution, this element provides a distribution and its corresponding probability.
For a "mixture" distribution, this element provides a distribution and its
corresponding probability.
:probability:
An attribute or ``pair`` that provides the probability of a univariate distribution within a "mixture" distribution.
An attribute or ``pair`` that provides the probability of a univariate
distribution within a "mixture" distribution.
:dist:
This sub-element of a ``pair`` element provides information on the corresponding univariate distribution.
This sub-element of a ``pair`` element provides information on the
corresponding univariate distribution.
:bias:
This optional element specifies a biased distribution for importance sampling.
For continuous distributions, the ``bias`` element should contain another
univariate distribution with the same support (interval) as the parent
distribution. For discrete distributions, the ``bias`` element should contain
floating point values corresponding to alternative probabilities for each
value/component to be used for biased sampling.
*Default*: None
---------------------------------------
``<source_rejection_fraction>`` Element
---------------------------------------
The ``<source_rejection_fraction>`` element specifies the minimum fraction of
external source sites that must be accepted when applying rejection sampling
based on constraints.
*Default*: 0.05
-------------------------
``<state_point>`` Element

2
docs/source/io_formats/statepoint.rst
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.

3
docs/source/io_formats/summary.rst
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'.

362
docs/source/methods/charged_particles_physics.rst Normal file
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

81
docs/source/methods/energy_deposition.rst
View file

@ -25,19 +25,38 @@ KERMA (Kinetic Energy Release in Materials) [Mack97]_ coefficients for reaction
:math:`\times` cross-section (e.g., eV-barn) and can be used much like a reaction
cross section for the purpose of tallying energy deposition.
KERMA coefficients can be computed using the energy-balance method with
a nuclear data processing code like NJOY, which performs the following
iteration over all reactions :math:`r` for all isotopes :math:`i`
requested
KERMA coefficients can be computed using the energy-balance method with a
nuclear data processing code like NJOY, which estimates the KERMA coefficients
using the following equation:
.. math::
k_{i, r}(E) = \left(E + Q_{i, r} - \bar{E}_{i, r, n}
k_{i, r}(E) = \left(E + Q_{i, r} - \sum\limits_x \bar{E}_{i, r, x}
\right)\sigma_{i, r}(E),
where the summation is over each secondary particle type :math:`x`. This
equation states that the energy deposited is equal to the energy of the incident
particle plus the reaction :math:`Q` value less the energy of secondary
particles that are transported away from the reaction site. For neutron
interactions, the energy-balance KERMA coefficient is
.. math::
k_{i, r}(E) = \left(E + Q_{i, r} - \sum\limits_x \bar{E}_{i, r, n}
- \bar{E}_{i, r, \gamma}\right)\sigma_{i, r}(E),
removing the energy of neutral particles (neutrons and photons) that are
transported away from the reaction site :math:`\bar{E}`, and the reaction
:math:`Q` value.
where :math:`\bar{E}_{i, r, n}` is the average energy of secondary neutrons and
:math:`\bar{E}_{i, r, \gamma}` is the average energy of secondary photons. For
photon and charged particle interactions the KERMA coefficient is
.. math::
:label: energy-balance-photon
k_{i, r}(E) = \left(E + Q_{i, r} - \sum\limits_x \bar{E}_{i, r, x}
\right)\sigma_{i, r}(E).
where the :math:`Q` value is zero for all interactions except for pair
production and positron annihilation.
-------
Fission
@ -120,7 +139,7 @@ run with :math:`N918` reflecting fission heating computed from NJOY.
This modified heating data is stored as the MT=901 reaction and will be scored
if ``heating-local`` is included in :attr:`openmc.Tally.scores`.
Coupled neutron-photon transport
Coupled Neutron-Photon Transport
--------------------------------
Here, OpenMC instructs ``heatr`` to assume that energy from photons is not
@ -138,6 +157,50 @@ Let :math:`N301` represent the total heating number returned from this
This modified heating data is stored as the MT=301 reaction and will be scored
if ``heating`` is included in :attr:`openmc.Tally.scores`.
Photons and Charged Particles
-----------------------------
In OpenMC, energy deposition from photons or charged particles is scored using
the energy balance method based on Equation :eq:`energy-balance-photon`. Special
consideration is given to electrons and positrons as described below.
+++++++++++++++++
Charged Particles
+++++++++++++++++
OpenMC tracks photons interaction by interaction so the energy deposited in each
collision is easily attributed back to the nuclide and reaction for which the
photon interacted with. Charged particles (electrons and photons) aren't tracked
in the same way. For charged particles, OpenMC assumes that all their energy
(less the energy of bremsstrahlung radiation) is deposited in the material in
which they were born. In this way it is harder to trace how much energy should
be attributed in each nuclide.
According to the CSDA approximation (see :ref:`ttb`) the energy deposited by a
charged particle with kinetic energy :math:`T` in the :math:`i`-th element can
be calculated as:
.. math::
E_{i} = \int_{0}^{R(T)} w_{i}S_{\text{col,i}} dx
where :math:`R(T)` is the CSDA range of the charged particle,
:math:`S_{\text{col},i}` is the collision stopping power of the charged particle
in the :math:`i`-th element and :math:`w_i` is the mass fraction of the
:math:`i`-th element. According to the Bethe formula the collision stopping
power of the :math:`i`-th element is proportional to :math:`Z_i/A_i`, so the
fractional collision stopping power from the :math:`i`-th element is:
.. math::
\frac{w_{i}S_{\text{col},i}(T)}{S_{\text{col}}(T)} =
\frac{\frac{w_{i}Z_{i}}{A_{i}}}{\sum_{i}\frac{w_{i}Z_{i}}{A_{i}}} =
\frac{\gamma_i Z_{i}}{\sum_{i}\gamma_i Z_{i}}.
where :math:`\gamma_i` is the atomic fraction of the :math:`i`-th element.
Therefore, the energy deposited by charged particles should be attributed to
a given element according to its fractional charge density.
----------
References
----------

3
docs/source/methods/index.rst
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

19
docs/source/methods/neutron_physics.rst
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

338
docs/source/methods/photon_physics.rst
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

7
docs/source/methods/random_ray.rst
View file

@ -109,9 +109,7 @@ terms on the right hand side.
In Equation :eq:`transport`, :math:`\psi` is the angular neutron flux. This
parameter represents the total distance traveled by all neutrons in a particular
direction inside of a control volume per second, and is often given in units of
:math:`1/(\text{cm}^{2} \text{s})`. As OpenMC does not support time dependence
in the random ray solver mode, we consider the steady state equation, where the
units of flux become :math:`1/\text{cm}^{2}`. The angular direction unit vector,
:math:`1/(\text{cm}^{2} \text{s})`. The angular direction unit vector,
:math:`\mathbf{\Omega}`, represents the direction of travel for the neutron. The
spatial position vector, :math:`\mathbf{r}`, represents the location within the
simulation. The neutron energy, :math:`E`, or speed in continuous space, is
@ -1052,7 +1050,8 @@ random ray and Monte Carlo, however.
regions. Thus, in the OpenMC implementation of random ray, particle sources
are restricted to being volumetric and isotropic, although different energy
spectrums are supported. Fixed sources can be applied to specific materials,
cells, or universes.
cells, or universes. Point sources are "smeared" to fill the volume of the
source region that contains the point source coordinate.
- **Inactive batches:** In Monte Carlo, use of a fixed source implies that all
batches are active batches, as there is no longer a need to develop a fission

130
docs/source/methods/tallies.rst
View file

@ -387,6 +387,130 @@ 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
+++++++++++++++
The figure of merit (FOM) is an indicator that accounts for both the statistical
uncertainty and the execution time and represents how much information is
obtained per unit time in the simulation. The FOM is defined as
.. math::
:label: figure_of_merit
FOM = \frac{1}{r^2 t},
where :math:`t` is the total execution time and :math:`r` is the relative error
defined as
.. math::
:label: relative_error
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). 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
++++++++++++++++++++
@ -494,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
@ -514,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

82
docs/source/methods/variance_reduction.rst
View file

@ -22,12 +22,14 @@ not experience a single scoring event, even after billions of analog histories.
Variance reduction techniques aim to either flatten the global uncertainty
distribution, such that all regions of phase space have a fairly similar
uncertainty, or to reduce the uncertainty in specific locations (such as a
detector). There are two strategies available in OpenMC for variance reduction:
the Monte Carlo MAGIC method and the FW-CADIS method. Both strategies work by
developing a weight window mesh that can be utilized by subsequent Monte Carlo
solves to split particles heading towards areas of lower flux densities while
terminating particles in higher flux regions---all while maintaining a fair
game.
detector). There are three strategies available in OpenMC for variance
reduction: weight windows generated via the MAGIC method or the FW-CADIS method,
and source biasing. Both weight windowing strategies work by developing a mesh
that can be utilized by subsequent Monte Carlo solves to split particles heading
towards areas of lower flux densities while terminating particles in higher flux
regions. In contrast, source biasing modifies source site sampling behavior to
preferentially track particles more likely to reach phase space regions of
interest.
------------
MAGIC Method
@ -132,3 +134,71 @@ aware of this.
:label: variance_fom
\text{FOM} = \frac{1}{\text{Time} \times \sigma^2}
.. _methods_source_biasing:
--------------
Source Biasing
--------------
In contrast to the previous two methods that introduce population controls
during transport, source biasing modifies the sampling of the external source
distribution. The basic premise of the technique is that for each spatial,
angular, energy, or time distribution of a source, an additional distribution
can be specified provided that the two share a common support (set of points
where the distribution is nonzero). Samples are then drawn from this "bias"
distribution, which can be chosen to preferentially direct particles towards
phase space regions of interest. In order to avoid biasing the tally results,
however, a weight adjustment is applied to each sampled site as described below.
Assume that the unbiased probability density function of a random variable
:math:`X:x \rightarrow \mathbb{R}` is given by :math:`f(x)`, but that using the
biased distribution :math:`g(x)` will result in a greater number of particle
trajectories reaching some phase space region of interest. Then a sample
:math:`x_0` may be drawn from :math:`g(x)` while maintaining a fair game,
provided that its weight is adjusted as:
.. math::
:label: source_bias
w = w_0 \times \frac{f(x_0)}{g(x_0)}
where :math:`w_0` is the weight of an unbiased sample from :math:`f(x)`,
typically unity.
Returning now to Equation :eq:`source_bias`, the requirement for common support
becomes evident. If :math:`\mathrm{supp} (g)` fully contains but is not
identical to :math:`\mathrm{supp} (f)`, then some samples from :math:`g(x)` will
correspond to points where :math:`f(x) = 0`. Thus these source sites would be
assigned a starting weight of 0, meaning the particles would be killed
immediately upon transport, effectively wasting computation time. Conversely, if
:math:`\mathrm{supp} (g)` is fully contained by but not identical to
:math:`\mathrm{supp} (f)`, the contributions of some regions outside
:math:`\mathrm{supp} (g)` will not be counted towards the integral, potentially
biasing the tally. The weight assigned to such points would be undefined since
:math:`g(x) = \mathbf{0}` at these points.
When an independent source is sampled in OpenMC, the particle's coordinate in
each variable of phase space :math:`(\mathbf{r},\mathbf{\Omega},E,t)` is
successively drawn from an independent probability distribution. Multiple
variables can be biased, in which case the resultant weight :math:`w` applied to
the particle is the product of the weights assigned from all sampled
distributions: space, angle, energy, and time, as shown in Equation
:eq:`tot_wgt`.
.. math::
:label: tot_wgt
w = w_r \times w_{\Omega} \times w_E \times w_t
Finally, source biasing and weight windows serve different purposes. Source
biasing changes how particles are born, allowing the initial source sites to be
sampled preferentially from important regions of phase space (space, angle,
energy, and time) with an accompanying weight adjustment. Weight windows, by
contrast, apply population control during transport (splitting and Russian
roulette) to help particles reach and contribute in important regions as they
move through the system. Because particle transport proceeds as usual after a
biased source is sampled, particle attenuation in optically thick regions
outside the source volume will not be affected by source biasing; in such
scenarios, transport biasing techniques such as weight windows are often more
effective.

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

@ -37,7 +37,6 @@ Simulation Settings
openmc.read_source_file
openmc.write_source_file
openmc.wwinp_to_wws
Material Specification
----------------------
@ -129,6 +128,7 @@ Constructing Tallies
openmc.SurfaceFilter
openmc.MeshFilter
openmc.MeshBornFilter
openmc.MeshMaterialFilter
openmc.MeshSurfaceFilter
openmc.EnergyFilter
openmc.EnergyoutFilter
@ -143,21 +143,31 @@ Constructing Tallies
openmc.SpatialLegendreFilter
openmc.SphericalHarmonicsFilter
openmc.TimeFilter
openmc.WeightFilter
openmc.ZernikeFilter
openmc.ZernikeRadialFilter
openmc.ParentNuclideFilter
openmc.ParticleFilter
openmc.RegularMesh
openmc.RectilinearMesh
openmc.CylindricalMesh
openmc.SphericalMesh
openmc.UnstructuredMesh
openmc.MeshMaterialVolumes
openmc.Trigger
openmc.TallyDerivative
openmc.Tally
openmc.Tallies
Meshes
------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclassinherit.rst
openmc.RegularMesh
openmc.RectilinearMesh
openmc.CylindricalMesh
openmc.SphericalMesh
openmc.UnstructuredMesh
Geometry Plotting
-----------------
@ -166,7 +176,8 @@ Geometry Plotting
:nosignatures:
:template: myclass.rst
openmc.Plot
openmc.SlicePlot
openmc.VoxelPlot
openmc.WireframeRayTracePlot
openmc.SolidRayTracePlot
openmc.Plots
@ -206,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.
@ -248,8 +262,16 @@ Variance Reduction
:template: myclass
openmc.WeightWindows
openmc.WeightWindowsList
openmc.WeightWindowGenerator
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.hdf5_to_wws
openmc.wwinp_to_wws
Coarse Mesh Finite Difference Acceleration

2
docs/source/pythonapi/capi.rst
View file

@ -89,8 +89,10 @@ Classes
SphericalMesh
SurfaceFilter
Tally
TemporarySession
UniverseFilter
UnstructuredMesh
WeightFilter
WeightWindows
ZernikeFilter
ZernikeRadialFilter

23
docs/source/pythonapi/deplete.rst
View file

@ -78,20 +78,18 @@ A minimal example for performing depletion would be:
>>> import openmc.deplete
>>> geometry = openmc.Geometry.from_xml()
>>> settings = openmc.Settings.from_xml()
>>> model = openmc.model.Model(geometry, settings)
>>> model = openmc.Model(geometry, settings)
# Representation of a depletion chain
>>> chain_file = "chain_casl.xml"
>>> operator = openmc.deplete.CoupledOperator(
... model, chain_file)
>>> operator = openmc.deplete.CoupledOperator(model, chain_file)
# Set up 5 time steps of one day each
>>> dt = [24 * 60 * 60] * 5
>>> power = 1e6 # constant power of 1 MW
# Deplete using mid-point predictor-corrector
>>> cecm = openmc.deplete.CECMIntegrator(
... operator, dt, power)
>>> cecm = openmc.deplete.CECMIntegrator(operator, dt, power)
>>> cecm.integrate()
Internal Classes and Functions
@ -208,14 +206,15 @@ total system energy.
The :class:`openmc.deplete.IndependentOperator` uses inner classes subclassed
from those listed above to perform similar calculations.
The following classes are used to define transfer rates to model continuous
removal or feed of nuclides during depletion.
The following classes are used to define external source rates or transfer rates
to model continuous removal or feed of nuclides during depletion.
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
transfer_rates.ExternalSourceRates
transfer_rates.TransferRates
Intermediate Classes
@ -288,6 +287,16 @@ the following abstract base classes:
abc.SIIntegrator
abc.DepSystemSolver
R2S Automation
--------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
R2SManager
D1S Functions
-------------

10
docs/source/pythonapi/mgxs.rst
View file

@ -11,6 +11,16 @@ Module Variables
.. autodata:: openmc.mgxs.GROUP_STRUCTURES
:annotation:
Functions
+++++++++
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.mgxs.convert_flux_groups
Classes
+++++++

7
docs/source/quickinstall.rst
View file

@ -35,6 +35,13 @@ you wish) with OpenMC installed.
conda create --name openmc-env openmc
conda activate openmc-env
If you are installing on macOS with an Apple silicon ARM-based processor, you
will also need to specify the `--platform` option:
.. code-block:: sh
conda create --name openmc-env --platform osx-64 openmc
You are now in a conda environment called `openmc-env` that has OpenMC
installed.

20
docs/source/releasenotes/0.15.2.rst Normal file
View file

@ -0,0 +1,20 @@
====================
What's New in 0.15.2
====================
.. currentmodule:: openmc
-------
Summary
-------
This is a hotfix release to fix an MPI-related bug that was inadvertently
introduced in the prior release.
---------------------------
Bug Fixes and Small Changes
---------------------------
- Remove errant ``openmc.Settings.random_ray`` check and removed of not useful warning in MG mode (`#3344 <https://github.com/openmc-dev/openmc/pull/3344>`_)
- Throw an error if a spherical harmonics order larger than 10 is provided. (`#3354 <https://github.com/openmc-dev/openmc/pull/3354>`_)
- Correcting the size of the displacement list in the SourceSite MPI interface object (`#3356 <https://github.com/openmc-dev/openmc/pull/3356>`_)

226
docs/source/releasenotes/0.15.3.rst Normal file
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>`_)

2
docs/source/releasenotes/index.rst
View file

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

3
docs/source/usersguide/data.rst
View file

@ -30,7 +30,8 @@ responsible for specifying one or more of the following:
Each of the above files can specified in several ways. In the Python API, a
:ref:`runtime configuration variable <usersguide_data_runtime>`
:data:`openmc.config` can be used to specify any of the above and is initialized
using a set of environment variables.
using a set of environment variables. Data configuration paths set in
:data:`openmc.config` will be expanded to absolute paths.
.. _usersguide_data_runtime:

207
docs/source/usersguide/decay_sources.rst
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
================================
@ -88,5 +235,5 @@ relevant tallies. This can be done with the aid of the
dose_tally = sp.get_tally(name='dose tally')
# Apply time correction factors
tally = d1s.apply_time_correction(tally, factors, time_index)
tally = d1s.apply_time_correction(dose_tally, factors, time_index)

97
docs/source/usersguide/geometry.rst
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
@ -413,11 +415,11 @@ to help figure out how to place universes::
Note that by default, hexagonal lattices are positioned such that each lattice
element has two faces that are parallel to the :math:`y` axis. As one example,
to create a three-ring lattice centered at the origin with a pitch of 10 cm
where all the lattice elements centered along the :math:`y` axis are filled with
universe ``u`` and the remainder are filled with universe ``q``, the following
code would work::
element has two faces that are perpendicular to the :math:`y` axis. As one
example, to create a three-ring lattice centered at the origin with a pitch of
10 cm where all the lattice elements centered along the :math:`y` axis are
filled with universe ``u`` and the remainder are filled with universe ``q``, the
following code would work::
hexlat = openmc.HexLattice()
hexlat.center = (0, 0)
@ -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

1
docs/source/usersguide/index.rst
View file

@ -22,6 +22,7 @@ essential aspects of using OpenMC to perform simulations.
plots
depletion
decay_sources
kinetics
scripts
processing
parallel

29
docs/source/usersguide/install.rst
View file

@ -35,6 +35,13 @@ you wish) with OpenMC installed.
conda create --name openmc-env openmc
conda activate openmc-env
If you are installing on macOS with an Apple silicon ARM-based processor, you
will also need to specify the `--platform` option:
.. code-block:: sh
conda create --name openmc-env --platform osx-64 openmc
You are now in a conda environment called `openmc-env` that has OpenMC
installed.
@ -221,7 +228,7 @@ Prerequisites
OpenMC's built-in plotting capabilities use the libpng library to produce
compressed PNG files. In the absence of this library, OpenMC will fallback
to writing PPM files, which are uncompressed and only supported by select
image viewers. libpng can be installed on Ddebian derivates with::
image viewers. libpng can be installed on Debian derivates with::
sudo apt install libpng-dev
@ -368,10 +375,6 @@ OPENMC_USE_DAGMC
should also be defined as `DAGMC_ROOT` in the CMake configuration command.
(Default: off)
OPENMC_USE_MCPL
Turns on support for reading MCPL_ source files and writing MCPL source points
and surface sources. (Default: off)
OPENMC_USE_LIBMESH
Enables the use of unstructured mesh tallies with libMesh_. (Default: off)
@ -380,6 +383,11 @@ OPENMC_USE_MPI
options, please see the `FindMPI.cmake documentation
<https://cmake.org/cmake/help/latest/module/FindMPI.html>`_.
OPENMC_FORCE_VENDORED_LIBS
Forces OpenMC to use the submodules located in the vendor directory, as
opposed to searching the system for already installed versions of those
modules.
To set any of these options (e.g., turning on profiling), the following form
should be used:
@ -515,10 +523,13 @@ to install the Python package in :ref:`"editable" mode <devguide_editable>`.
Prerequisites
-------------
The Python API works with Python 3.8+. In addition to Python itself, the API
relies on a number of third-party packages. All prerequisites can be installed
using Conda_ (recommended), pip_, or through the package manager in most Linux
distributions.
In addition to Python itself, the OpenMC Python API relies on a number of
third-party packages. All prerequisites can be installed using Conda_
(recommended), pip_, or through the package manager in most Linux distributions.
The current required Python version and up-to-date list of package dependencies
can be found in the `pyproject.toml <https://github.com/openmc-dev/openmc/blob/develop/pyproject.toml>`_
file in the root directory of the OpenMC repository. An overview of these
dependencies is provided below.
.. admonition:: Required
:class: error

133
docs/source/usersguide/kinetics.rst Normal file
View file

@ -0,0 +1,133 @@
.. _kinetics:
===================
Kinetics parameters
===================
OpenMC has the capability to estimate the following adjoint-weighted effective
generation time :math:`\Lambda_{\text{eff}}` and the effective delayed neutron
fraction :math:`\beta_{\text{eff}}`. These parameters are calculated using the
iterated fission probability (IFP) method [Hurwitz_1964]_ based on a similar
approach as in `Serpent 2 <https://doi.org/10.1016/j.anucene.2013.10.032>`_. The
implementation in OpenMC is limited to eigenvalue calculations and is described
in more details in [Dorville_2025]_.
----------------------------------
Iterated Fission Probability (IFP)
----------------------------------
With IFP, additional information needs to be recorded during the simulation
compared to a typical eigenvalue calculation. OpenMC stores an additional
set of values (neutron lifetime or delayed neutron group number for
:math:`\Lambda_{\text{eff}}` or :math:`\beta_{\text{eff}}`, respectively)
for every fission neutron simulated. Each set of values corresponds to
the values that are associated to the :math:`N_{\text{gen}}` direct ancestors
of any given fission neutron.
:math:`N_{\text{gen}}` is referred to as the number of generations in the
IFP method and corresponds to the number of generations between the birth of
a fission neutron and the time its score is added to the IFP tally. By default,
OpenMC considers 10 generations but this value can be modified by the user via
the ``ifp_n_generation`` settings in the Python API::
settings.ifp_n_generation = 5
``ifp_n_generation`` should be greater than 0, but should also be lower than
or equal to the number of inactive batches declared for the calculation.
The respect of these constraints is verified by OpenMC before any calculation.
OpenMC will automatically detect the type of data that needs to be stored based
on the tally scores selected by the user. This guarantees that only information
of interest are stored during a simulation and avoids using extra memory when
only one parameter is needed. The following table shows the tally scores that
are needed to compute kinetics parameters in OpenMC:
.. table:: **OpenMC tally scores needed to calculate adjoint-weighted kinetics parameters**
:align: center
=============================== ============================ ========================== ========
OpenMC tally score \\ Parameter :math:`\Lambda_{\text{eff}}` :math:`\beta_{\text{eff}}` Both
=============================== ============================ ========================== ========
``ifp-time-numerator`` X X
``ifp-beta-numerator`` X X
``ifp-denominator`` X X X
=============================== ============================ ========================== ========
|
.. note:: Because the memory footprint of additional data is generally non-negligible
with IFP, it is recommended to choose the value for ``ifp_n_generation`` carefully.
For example, using one generation for both kinetics parameters corresponds to store
one additional integer (for the delayed neutron group number used with
:math:`\beta_{\text{eff}}`) and one floating point value (for the neutron lifetime
used with :math:`\Lambda_{\text{eff}}`) for every fission neutron simulated once the
asymptotic regime is reached.
-----------------------------
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::
tally = openmc.Tally(name="ifp-scores")
tally.scores = [
"ifp-time-numerator",
"ifp-beta-numerator",
"ifp-denominator"
]
The effective generation time :math:`\Lambda_{\text{eff}}` is calculated
by dividing the result of the ``ifp-time-numerator`` score by the one obtained
for ``ifp-denominator`` and by the :math:`k_{\text{eff}}` of the simulation:
.. math::
:label: lambda_eff
\Lambda_{\text{eff}} = \frac{S_{\text{ifp-time-numerator}}}{S_{\text{ifp-denominator}} \times k_{\text{eff}}}
The effective delayed neutron fraction :math:`\beta_{\text{eff}}` is calculated
by dividing the result of the ``ifp-beta-numerator`` score by the one obtained
for ``ifp-denominator``:
.. math::
:label: beta_eff
\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
.. [Hurwitz_1964] H. Hurwitz Jr., "Naval Reactors Physics Handbook", volume 1, p. 864.
Radkowsky, A. (Ed.), Naval Reactors, Division of Reactor Development, U.S.
Atomic Energy Commission (1964).
.. [Dorville_2025] J. Dorville, L. Labrie-Cleary, and P. K. Romano, "Implementation
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.

46
docs/source/usersguide/plots.rst
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)

95
docs/source/usersguide/processing.rst
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
@ -91,3 +91,82 @@ from a statepoint file, the ``openmc.statepoint`` module can be used. An
`example notebook`_ demontrates how to analyze and plot source information.
.. _example notebook: https://nbviewer.jupyter.org/github/openmc-dev/openmc-notebooks/blob/main/post-processing.ipynb
------------------------
VTK Mesh File Generation
------------------------
VTK files of OpenMC meshes can be created using the
:meth:`openmc.Mesh.write_data_to_vtk` method. Data can be applied to the
elements of the resulting mesh from mesh filter objects. This data can be
provided either as a flat array or, in the case of structured meshes
(:class:`~openmc.RegularMesh`, :class:`~openmc.RectilinearMesh`,
:class:`~openmc.CylindricalMesh`, or :class:`SphericalMesh`), the data can be
shaped with dimensions that match the dimensions of the mesh itself.
.. image:: ../_images/sphere-mesh-vtk.png
:width: 400px
:align: center
:alt: OpenMC spherical mesh exported to VTK
For all mesh types, if a flat data array is provided to the mesh, it is expected
that the data is ordered in the same ordering as the :attr:`openmc.Mesh.indices`
for that mesh object. When providing data directly from a tally, as shown below,
a flat array for a given dataset can be passed directly to this method.
::
# create model above
# create a mesh tally
mesh = openmc.RegularMesh()
mesh.dimension = [10, 20, 30]
mesh.lower_left = [-5, -10, -15]
mesh.upper_right = [5, 10, 15]
mesh_filter = openmc.MeshFilter(mesh)
tally = openmc.Tally()
tally.filters = [mesh_filter]
tally.scores = ['flux']
model.tallies = [tally]
model.run(apply_tally_results=True)
# provide the data as-is to the method
mesh.write_data_to_vtk('flux.vtk', {'flux-mean': tally.mean})
The :class:`~openmc.Tally` object also provides a way to expand the dimensions
of the mesh filter into a meaningful form where indexing the mesh filter
dimensions results in intuitive slicing of structured meshes by setting
``expand_dims=True`` when using :meth:`openmc.Tally.get_reshaped_data`. This
reshaping does cause flat indexing of the data to change, however. As noted
above, provided datasets are allowed to be shaped so long as such datasets have
shapes that match the mesh dimensions. The ability to pass datasets in this way
is useful when additional filters are applied to a tally. The example below
demonstrates such a case for tally with both a :class:`~openmc.MeshFilter` and
:class:`~openmc.EnergyFilter` applied.
::
# create model above
# create a mesh tally with energy filter
mesh = openmc.RegularMesh()
mesh.dimension = [10, 20, 30]
mesh.lower_left = [-5, -10, -15]
mesh.upper_right = [5, 10, 15]
mesh_filter = openmc.MeshFilter(mesh)
energy_filter = openmc.EnergyFilter([0.0, 1.0, 20.0e6])
tally = openmc.Tally()
tally.filters = [mesh_filter, energy_filter]
tally.scores = ['flux']
model.tallies = [tally]
model.run(apply_tally_results=True)
# get the data with mesh dimensions expanded, squeeze out length-one dimensions (nuclides, scores)
flux = tally.get_reshaped_data(expand_dims=True).squeeze() # shape: (10, 20, 30, 2)
# write the lowest energy group to a VTK file
mesh.write_data_to_vtk('flux-group1.vtk', datasets={'flux-mean': flux[..., 0]})

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

@ -11,6 +11,76 @@ active batches <usersguide_batches>`. However, there are a couple of settings
that are unique to the random ray solver and a few areas that the random ray
run strategy differs, both of which will be described in this section.
.. _quick_start:
-----------
Quick Start
-----------
While this page contains a comprehensive guide to the random ray solver and
its various parameters, the process of converting an existing continuous energy
Monte Carlo model to a random ray model can be largely automated via convenience
functions in OpenMC's Python interface::
# Define continuous energy model as normal
model = openmc.Model()
...
# Convert model to multigroup (will auto-generate MGXS library if needed)
model.convert_to_multigroup()
# Convert model to random ray and initialize random ray parameters
# to reasonable defaults based on the specifics of the geometry
model.convert_to_random_ray()
# (Optional) Overlay source region decomposition mesh to improve fidelity of the
# random ray solver. Adjust 'n' for fidelity vs runtime.
n = 100
mesh = openmc.RegularMesh()
mesh.dimension = (n, n, n)
mesh.lower_left = model.geometry.bounding_box.lower_left
mesh.upper_right = model.geometry.bounding_box.upper_right
model.settings.random_ray['source_region_meshes'] = [(mesh, [model.geometry.root_universe])]
# (Optional) Improve fidelity of the random ray solver by enabling linear sources
model.settings.random_ray['source_shape'] = 'linear'
# (Optional) Increase the number of rays/batch, to reduce uncertainty
model.settings.particles = 500
The above strategy first converts the continuous energy model to a multigroup
one using the :meth:`openmc.Model.convert_to_multigroup` method. By default,
this will internally run a coarsely converged continuous energy Monte Carlo
simulation to produce an estimated multigroup macroscopic cross section set for
each material specified in the model, and store this data into a multigroup
cross section library file (``mgxs.h5``) that can be used by the random ray
solver.
The :meth:`openmc.Model.convert_to_random_ray` method enables random ray mode
and performs an analysis of the model geometry to determine reasonable values
for all required parameters. If default behavior is not satisfactory, the user
can manually adjust the settings in the :attr:`~openmc.Settings.random_ray`
dictionary in the :class:`openmc.Settings` as described in the sections below.
Finally a few optional steps are shown. The first (recommended) step overlays a
mesh over the geometry to create smaller source regions so that source
resolution improves and the random ray solver becomes more accurate. Varying the
mesh resolution can be used to trade off between accuracy and runtime.
High-fidelity fission reactor simulation may require source region sizes below 1
cm, while larger fixed source problems with some tolerance for error may be able
to use source regions of 10 or 100 cm.
We also enable linear sources, which can improve the accuracy of the random ray
solver and/or allow for a much coarser mesh resolution to be overlaid. Finally,
the number of rays per batch is adjusted. The goal here is to ensure that the
source region miss rate is below 1%, which is reported by OpenMC at the end of
the simulation (or before via a warning if it is very high).
.. warning::
If using a mesh filter for tallying or weight window generation, ensure that
the same mesh is used for source region decomposition via
``model.settings.random_ray['source_region_meshes']``.
------------------------
Enabling Random Ray Mode
------------------------
@ -557,29 +627,144 @@ variety of problem types (or through a multidimensional parameter sweep of
design variables) with only modest errors and at greatly reduced cost as
compared to using only continuous energy Monte Carlo.
~~~~~~~~~~~~
The Easy Way
~~~~~~~~~~~~
The easiest way to generate a multigroup cross section library is to use the
:meth:`openmc.Model.convert_to_multigroup` method. This method will
automatically output a multigroup cross section library file (``mgxs.h5``) from
a continuous energy Monte Carlo model and alter the material definitions in the
model to use these multigroup cross sections. An example is given below::
# Assume we already have a working continuous energy model
model.convert_to_multigroup(
method="material_wise",
groups="CASMO-2",
nparticles=2000,
overwrite_mgxs_library=False,
mgxs_path="mgxs.h5",
correction=None,
source_energy=None
)
The most important parameter to set is the ``method`` parameter, which can be
either "stochastic_slab", "material_wise", or "infinite_medium". An overview
of these methods is given below:
.. list-table:: Comparison of Automatic MGXS Generation Methods
:header-rows: 1
:widths: 10 30 30 30
* - Method
- Description
- Pros
- Cons
* - ``material_wise`` (default)
- * Higher Fidelity
* Runs a CE simulation with the original geometry and source, tallying
cross sections with a material filter.
- * Typically the most accurate of the three methods
* Accurately captures (averaged over the full problem domain)
both spatial and resonance self shielding effects
- * Potentially slower as the full geometry must be run
* If a material is only present far from the source and doesn't get tallied
to in the CE simulation, the MGXS will be zero for that material.
* - ``stochastic_slab``
- * Medium Fidelity
* Runs a CE simulation with a greatly simplified geometry, where materials
are randomly assigned to layers in a 1D "stochastic slab sandwich" geometry
- * Still captures resonant self shielding and resonance effects between materials
* Fast due to the simplified geometry
* Able to produce cross section data for all materials, regardless of how
far they are from the source in the original geometry
- * Does not capture most spatial self shielding effects, e.g., no lattice physics.
* - ``infinite_medium``
- * Lower Fidelity
* Runs one CE simulation per material independently. Each simulation is just
an infinite medium slowing down problem, with an assumed external source term.
- * Simple
- * Poor accuracy (no spatial information, no lattice physics, no resonance effects
between materials)
* May hang if a material has a k-infinity greater than 1.0
When selecting a non-default energy group structure, you can manually define
group boundaries or specify the name of a known group structure (a list of which
can be found at :data:`openmc.mgxs.GROUP_STRUCTURES`). The ``nparticles``
parameter can be adjusted upward to improve the fidelity of the generated cross
section library. The ``correction`` parameter can be set to ``"P0"`` to enable
P0 transport correction. The ``overwrite_mgxs_library`` parameter can be set to
``True`` to overwrite an existing MGXS library file, or ``False`` to skip
generation and use an existing library file.
.. note::
MGXS transport correction (via setting the ``correction`` parameter in the
:meth:`openmc.Model.convert_to_multigroup` method to ``"P0"``) may
result in negative in-group scattering cross sections, which can cause
numerical instability. To mitigate this, during a random ray solve OpenMC
will automatically apply
`diagonal stabilization <https://doi.org/10.1016/j.anucene.2018.10.036>`_
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
useful in that they can provide a good starting point for a random ray
simulation, and if more fidelity is needed the user may wish to follow the
instructions below or experiment with transport correction techniques to improve
the fidelity of the generated MGXS data.
~~~~~~~~~~~~
The Hard Way
~~~~~~~~~~~~
We give here a quick summary of how to produce a multigroup cross section data
file (``mgxs.h5``) from a starting point of a typical continuous energy Monte
Carlo input file. Notably, continuous energy input files define materials as a
mixture of nuclides with different densities, whereas multigroup materials are
simply defined by which name they correspond to in a ``mgxs.h5`` library file.
Carlo model. Notably, continuous energy models define materials as a mixture of
nuclides with different densities, whereas multigroup materials are simply
defined by which name they correspond to in a ``mgxs.h5`` library file.
To generate the cross section data, we begin with a continuous energy Monte
Carlo input deck and add in the required tallies that will be needed to generate
our library. In this example, we will specify material-wise cross sections and a
two group energy decomposition::
Carlo model and add in the tallies that are needed to generate our library. In
this example, we will specify material-wise cross sections and a two-group
energy decomposition::
# Define geometry
...
...
geometry = openmc.Geometry()
...
...
# Initialize MGXS library with a finished OpenMC geometry object
mgxs_lib = openmc.mgxs.Library(geometry)
# Pick energy group structure
groups = openmc.mgxs.EnergyGroups(openmc.mgxs.GROUP_STRUCTURES['CASMO-2'])
groups = openmc.mgxs.EnergyGroups('CASMO-2')
mgxs_lib.energy_groups = groups
# Disable transport correction
@ -587,7 +772,7 @@ two group energy decomposition::
# Specify needed cross sections for random ray
mgxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission',
'nu-scatter matrix', 'multiplicity matrix', 'chi']
'nu-scatter matrix', 'multiplicity matrix', 'chi']
# Specify a "cell" domain type for the cross section tally filters
mgxs_lib.domain_type = "material"
@ -606,7 +791,7 @@ two group 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()
@ -614,13 +799,13 @@ two group energy decomposition::
...
When selecting an energy decomposition, you can manually define group boundaries
or pick out a group structure already known to OpenMC (a list of which can be
found at :class:`openmc.mgxs.GROUP_STRUCTURES`). Once the above input deck has
been run, the resulting statepoint file will contain the needed flux and
reaction rate tally data so that a MGXS library file can be generated. Below is
the postprocessing script needed to generate the ``mgxs.h5`` library file given
a statepoint file (e.g., ``statepoint.100.h5``) file and summary file (e.g.,
``summary.h5``) that resulted from running our previous example::
or specify the name of known group structure (a list of which can be found at
:data:`openmc.mgxs.GROUP_STRUCTURES`). Once the above model has been run, the
resulting statepoint file will contain the needed flux and reaction rate tally
data so that a MGXS library file can be generated. Below is the postprocessing
script needed to generate the ``mgxs.h5`` library file given a statepoint file
(e.g., ``statepoint.100.h5``) file and summary file (e.g., ``summary.h5``) that
resulted from running our previous example::
import openmc
@ -628,10 +813,7 @@ a statepoint file (e.g., ``statepoint.100.h5``) file and summary file (e.g.,
geom = summary.geometry
mats = summary.materials
statepoint_filename = 'statepoint.100.h5'
sp = openmc.StatePoint(statepoint_filename)
groups = openmc.mgxs.EnergyGroups(openmc.mgxs.GROUP_STRUCTURES['CASMO-2'])
groups = openmc.mgxs.EnergyGroups('CASMO-2')
mgxs_lib = openmc.mgxs.Library(geom)
mgxs_lib.energy_groups = groups
mgxs_lib.correction = None
@ -653,10 +835,10 @@ a statepoint file (e.g., ``statepoint.100.h5``) file and summary file (e.g.,
# Construct all tallies needed for the multi-group cross section library
mgxs_lib.build_library()
mgxs_lib.load_from_statepoint(sp)
with openmc.StatePoint('statepoint.100.h5') as sp:
mgxs_lib.load_from_statepoint(sp)
names = []
for mat in mgxs_lib.domains: names.append(mat.name)
names = [mat.name for mat in mgxs_lib.domains]
# Create a MGXS File which can then be written to disk
mgxs_file = mgxs_lib.create_mg_library(xs_type='macro', xsdata_names=names)
@ -665,8 +847,8 @@ a statepoint file (e.g., ``statepoint.100.h5``) file and summary file (e.g.,
mgxs_file.export_to_hdf5("mgxs.h5")
Notably, the postprocessing script needs to match the same
:class:`openmc.mgxs.Library` settings that were used to generate the tallies,
but otherwise is able to discern the rest of the simulation details from the
:class:`openmc.mgxs.Library` settings that were used to generate the tallies but
is otherwise able to discern the rest of the simulation details from the
statepoint and summary files. Once the postprocessing script is successfully
run, the ``mgxs.h5`` file can be loaded by subsequent runs of OpenMC.
@ -701,11 +883,11 @@ multigroup library instead of defining their isotopic contents, as::
water_data = openmc.Macroscopic('Hot borated water')
# Instantiate some Materials and register the appropriate Macroscopic objects
fuel= openmc.Material(name='UO2 (2.4%)')
fuel = openmc.Material(name='UO2 (2.4%)')
fuel.set_density('macro', 1.0)
fuel.add_macroscopic(fuel_data)
water= openmc.Material(name='Hot borated water')
water = openmc.Material(name='Hot borated water')
water.set_density('macro', 1.0)
water.add_macroscopic(water_data)
@ -763,9 +945,12 @@ Monte Carlo solver.
Currently, all of the following conditions must be met for the particle source
to be valid in random ray mode:
- One or more domain ids must be specified that indicate which cells, universes,
or materials the source applies to. This implicitly limits the source type to
being volumetric. This is specified via the ``domains`` constraint placed on the
- Either a point source must be used, or a domain constraint must be specified
that indicates which cells, universes, or materials the source applies to. In
either case, this implicitly limits the source type to being volumetric, as
even in the point source case the source will be "smeared" throughout the
source region that contains the point source coordinate. A source domain is
specified via the ``domains`` constraint placed on the
:class:`openmc.IndependentSource` Python class.
- The source must be isotropic (default for a source)
- The source must use a discrete (i.e., multigroup) energy distribution. The
@ -946,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])
@ -1030,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])

1
docs/source/usersguide/scripts.rst
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

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

@ -272,6 +272,12 @@ option::
settings.source = [src1, src2]
settings.uniform_source_sampling = True
Additionally, sampling from an :class:`openmc.IndependentSource` may be biased
for local or global variance reduction by modifying the
:attr:`~openmc.IndependentSource.bias` attribute of each of its four main
distributions. Further discussion of source biasing can be found in
:ref:`source_biasing`.
Finally, the :attr:`IndependentSource.particle` attribute can be used to
indicate the source should be composed of particles other than neutrons. For
example, the following would generate a photon source::
@ -756,6 +762,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
-----------------------

18
docs/source/usersguide/tallies.rst
View file

@ -322,6 +322,24 @@ The following tables show all valid scores:
| |particle. Note that this score can only be combined|
| |with a cell filter and an energy filter. |
+----------------------+---------------------------------------------------+
|ifp-time-numerator |Adjoint-weighted lifetime of neutron produced by |
| |fission in units of seconds per source particle. |
| |This score is used to compute kinetics parameters |
| |using the iterated fission probability (IFP) |
| |method. |
+----------------------+---------------------------------------------------+
|ifp-beta-numerator |Adjoint-weighted number of delayed fission events |
| |in units of number of delayed fission event per |
| |source particle. This score is used to compute |
| |kinetics parameters using the iterated fission |
| |probability (IFP) method. |
+----------------------+---------------------------------------------------+
|ifp-denominator |Weights corresponding to the number of fission |
| |events in units of number of fission event per |
| |source particle. This score is used to compute |
| |kinetics parameters using the iterated fission |
| |probability (IFP) method. |
+----------------------+---------------------------------------------------+
.. _usersguide_tally_normalization:

235
docs/source/usersguide/variance_reduction.rst
View file

@ -5,10 +5,12 @@ Variance Reduction
==================
Global variance reduction in OpenMC is accomplished by weight windowing
techniques. OpenMC is capable of generating weight windows using either the
MAGIC or FW-CADIS methods. Both techniques will produce a ``weight_windows.h5``
file that can be loaded and used later on. In this section, we break down the
steps required to both generate and then apply weight windows.
or source biasing techniques, the latter of which additionally provides a
local variance reduction capability. OpenMC is capable of generating weight
windows using either the MAGIC or FW-CADIS methods. Both techniques will
produce a ``weight_windows.h5`` file that can be loaded and used later on. In
this section, we first break down the steps required to generate and apply
weight windows, then describe how source biasing may be applied.
.. _ww_generator:
@ -51,7 +53,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.
@ -88,59 +90,65 @@ random ray mode can be found in the :ref:`Random Ray User Guide <random_ray>`.
ray solver. A high level overview of the current workflow for generation of
weight windows with FW-CADIS using random ray is given below.
1. Produce approximate multigroup cross section data (stored in a ``mgxs.h5``
library). There is more information on generating multigroup cross sections
via OpenMC in the :ref:`multigroup materials <create_mgxs>` user guide, and a
specific example of generating cross section data for use with random ray in
the :ref:`random ray MGXS guide <mgxs_gen>`.
1. Begin by making a deepy copy of your continuous energy Python model and then
convert the copy to be multigroup and use the random ray transport solver.
The conversion process can largely be automated as described in more detail
in the :ref:`random ray quick start guide <quick_start>`, summarized below::
2. Make a copy of your continuous energy Python input file. You'll edit the new
file to work in multigroup mode with random ray for producing weight windows.
# Define continuous energy model
ce_model = openmc.pwr_pin_cell() # example, replace with your model
3. Adjust the material definitions in your new multigroup Python file to utilize
the multigroup cross sections instead of nuclide-wise continuous energy data.
There is a specific example of making this conversion in the :ref:`random ray
MGXS guide <mgxs_gen>`.
# Make a copy to convert to multigroup and random ray
model = copy.deepcopy(ce_model)
4. Configure OpenMC to run in random ray mode (by adding several standard random
ray input flags and settings to the :attr:`openmc.Settings.random_ray`
dictionary). More information can be found in the :ref:`Random Ray User
Guide <random_ray>`.
# Convert model to multigroup (will auto-generate MGXS library if needed)
model.convert_to_multigroup()
5. Add in a :class:`~openmc.WeightWindowGenerator` in a similar manner as for
# Convert model to random ray and initialize random ray parameters
# to reasonable defaults based on the specifics of the geometry
model.convert_to_random_ray()
# (Optional) Overlay source region decomposition mesh to improve fidelity of the
# random ray solver. Adjust 'n' for fidelity vs runtime.
n = 10
mesh = openmc.RegularMesh()
mesh.dimension = (n, n, n)
mesh.lower_left = model.geometry.bounding_box.lower_left
mesh.upper_right = model.geometry.bounding_box.upper_right
model.settings.random_ray['source_region_meshes'] = [(mesh, [model.geometry.root_universe])]
# (Optional) Improve fidelity of the random ray solver by enabling linear sources
model.settings.random_ray['source_shape'] = 'linear'
# (Optional) Increase the number of rays/batch, to reduce uncertainty
model.settings.particles = 500
If you need to improve the fidelity of the MGXS library, there is more
information on generating multigroup cross sections via OpenMC in the
:ref:`random ray MGXS guide <mgxs_gen>`.
2. Add in a :class:`~openmc.WeightWindowGenerator` in a similar manner as for
MAGIC generation with Monte Carlo and set the :attr:`method` attribute set to
``"fw_cadis"``::
# Define weight window spatial mesh
ww_mesh = openmc.RegularMesh()
ww_mesh.dimension = (10, 10, 10)
ww_mesh.lower_left = (0.0, 0.0, 0.0)
ww_mesh.upper_right = (100.0, 100.0, 100.0)
# Create weight window object and adjust parameters
# Create weight window object and adjust parameters, using the same mesh
# we used for source region decomposition
wwg = openmc.WeightWindowGenerator(
method='fw_cadis',
mesh=ww_mesh,
max_realizations=settings.batches
mesh=mesh
)
# Add generator to openmc.settings object
settings.weight_window_generators = wwg
.. warning::
If using FW-CADIS weight window generation, ensure that the selected weight
window mesh does not subdivide any source regions in the problem. This can
be ensured by assigning the weight window tally mesh to the root universe so
as to create source region boundaries that conform to the mesh, as in the
example below.
be ensured by using the same mesh for both source region subdivision (i.e.,
assigning to ``model.settings.random_ray['source_region_meshes']``) and for
weight window generation.
::
root = model.geometry.root_universe
settings.random_ray['source_region_meshes'] = [(ww_mesh, [root])]
6. When running your multigroup random ray input deck, OpenMC will automatically
3. When running your multigroup random ray input deck, OpenMC will automatically
run a forward solve followed by an adjoint solve, with a
``weight_windows.h5`` file generated at the end. The ``weight_windows.h5``
file will contain FW-CADIS generated weight windows. This file can be used in
@ -155,7 +163,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.hdf5_to_wws('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
@ -166,3 +174,148 @@ Weight window mesh information is embedded into the weight window file, so the
mesh does not need to be redefined. Monte Carlo solves that load a weight window
file as above will utilize weight windows to reduce the variance of the
simulation.
.. _source_biasing:
--------------
Source Biasing
--------------
In fixed source problems, source biasing provides a means to reduce the variance
on global or localized responses, depending on the biasing scheme. In either
case, the premise of the method is to sample source sites from a biased
distribution that directs a larger fraction of the simulated histories towards
phase space regions of interest than would be found there under analog sampling.
In order to preserve an unbiased estimate of the tally mean, the weight of these
with analog sampling, divided by the probability assigned by the biased
distribution. While the assignment of statistical weights is outlined in the
:ref:`methods section <methods_source_biasing>`, this section demonstrates the
implementation of source biasing to problems in OpenMC.
Source biasing in OpenMC is accomplished by applying a distribution to the
:attr:`bias` attribute of one or more of the univariate or independent
multivariate distributions which make up an :class:`~openmc.IndependentSource`
instance as follows::
# First create the biased distribution
biased_dist = openmc.stats.PowerLaw(a=0, b=3, n=3)
# Construct a new distribution with the bias applied
dist = openmc.stats.PowerLaw(a=0, b=3, n=2, bias=biased_dist)
# The bias attribute can also be set on an existing "analog" distribution:
sphere_dist = openmc.stats.spherical_uniform(r_outer=3)
sphere_dist.r.bias = biased_dist
Univariate distributions may be sampled via the Python API, returning the
sample(s) along with the associated weight(s)::
sample_vec, wgt_vec = dist.sample(n_samples=100)
Here, if the distribution is unbiased, the weight of each sample will be unity.
Finally, :class:`~openmc.IndependentSource` instances can be constructed with
biased distributions::
# Create a source with a biased spatial distribution
source = openmc.IndependentSource(space=sphere_dist)
During the simulation, source sites are then sampled using the biased
distributions where available and given starting statistical weights
corresponding to the cumulative product of the weights assigned by each
distribution in the source object. Hence multiple source variables (e.g.,
direction and energy) may be biased and the resulting source sites will have
their weights adjusted accordingly.
.. note::
Combining source biasing with weight windows can be a powerful variance
reduction technique if each is constructed appropriately for the response
of interest. For example, if a source biasing scheme is devised for
variance reduction of a specific localized response, the user may be able
to specify their own weight window structure that results in more efficient
transport than if weight windows were generated by either of OpenMC's
automatic weight window generators, which are intended for global variance
reduction.
Biased distributions that could result in degenerate weight mappings are not
recommended; this is most commonly seen when biasing the :math:`\phi`-coordinate
of spherical or cylindrical independent multivariate distributions. In such
cases degenerate behavior will be observed at the pole about which :math:`\phi`
is measured, with all values of :math:`\phi` (hence many possible statistical
weights) mapping to the same point for :math:`r=0` or :math:`\mu=0`, and large
weight gradients in the vicinity. In most cases requiring a spherical
independent source, it would be preferable to reorient the reference vector of
the distribution such that biasing could be applied to the
:math:`\mu`-coordinate instead.
When biasing a distribution, care should also be taken to ensure that both the
unbiased and biased distribution share a common support---that is, every region
of phase space mapped to a nonzero probability density by the unbiased
distribution should likewise map to nonzero probability under the biased
distribution, and vice versa. In OpenMC, this places restrictions on the set of
compatible distributions that may be used to bias sampling of each distribution
type. The following table summarizes the method for each distribution in OpenMC
that permits biased sampling.
.. list-table:: **Distributions that support biased sampling**
:header-rows: 1
:widths: 35 65
* - Discrete Univariate PDFs
- Biasing Method
* - :class:`openmc.stats.Discrete`
- Apply a vector of alternative probabilities to the :attr:`bias`
attribute
.. list-table::
:header-rows: 1
:widths: 35 65
* - Continuous Univariate PDFs
- Biasing Method
* - :class:`openmc.stats.Uniform`,
:class:`openmc.stats.PowerLaw`,
:class:`openmc.stats.Maxwell`,
:class:`openmc.stats.Watt`,
:class:`openmc.stats.Normal`,
:class:`openmc.stats.Tabular`
- Apply a second, unbiased continous univariate PDF to the :attr:`bias`
attribute, ensuring that the :attr:`support` attribute of each
distribution is the same
.. list-table::
:header-rows: 1
:widths: 35 65
* - Mixed Univariate PDFs
- Biasing Method
* - :class:`openmc.stats.Mixture`
- May be constructed from multiple biased univariate distributions, or a
second, unbiased continous univariate PDF may be applied to the
:attr:`bias` attribute
.. list-table::
:header-rows: 1
:widths: 35 65
* - Discrete Multivariate PDFs
- Biasing Method
* - :class:`openmc.stats.PointCloud`,
:class:`openmc.stats.MeshSpatial`
- Apply a vector of the new relative probabilities of each point or mesh
element under biased sampling to the :attr:`bias` attribute
.. list-table::
:header-rows: 1
:widths: 35 65
* - Continuous Multivariate PDFs
- Biasing Method
* - :class:`openmc.stats.CartesianIndependent`,
:class:`openmc.stats.CylindricalIndependent`,
:class:`openmc.stats.SphericalIndependent`,
:class:`openmc.stats.PolarAzimuthal`
- Construct from biased univariate distributions for :attr:`x`, :attr:`y`,
:attr:`z`, etc.
* - :class:`openmc.stats.Isotropic`
- Apply an unbiased :class:`openmc.stats.PolarAzimuthal` to the
:attr:`bias` attribute

4
examples/lattice/hexagonal/build_xml.py
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]

2
examples/lattice/nested/build_xml.py
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]

2
examples/lattice/simple/build_xml.py
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]

101
examples/pincell_pulsed/run_pulse.py Normal file
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()

3
examples/pincell_random_ray/build_xml.py
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])

10
include/openmc/bank.h
View file

@ -20,8 +20,18 @@ 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;
extern vector<vector<double>> ifp_source_lifetime_bank;
extern vector<vector<int>> ifp_fission_delayed_group_bank;
extern vector<vector<double>> ifp_fission_lifetime_bank;
extern vector<int64_t> progeny_per_particle;
} // namespace simulation

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

20
include/openmc/boundary_condition.h
View file

@ -111,6 +111,10 @@ public:
std::string type() const override { return "periodic"; }
int i_surf() const { return i_surf_; }
int j_surf() const { return j_surf_; }
protected:
int i_surf_;
int j_surf_;
@ -134,18 +138,28 @@ 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_;
//! Do we need to flip surfaces senses when applying the transformation?
bool flip_sense_;
//! 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

46
include/openmc/bounding_box.h
View file

@ -13,12 +13,19 @@ namespace openmc {
//==============================================================================
struct BoundingBox {
double xmin = -INFTY;
double xmax = INFTY;
double ymin = -INFTY;
double ymax = INFTY;
double zmin = -INFTY;
double zmax = INFTY;
Position min = {-INFTY, -INFTY, -INFTY};
Position max = {INFTY, INFTY, INFTY};
// Constructors
BoundingBox() = default;
BoundingBox(Position min_, Position max_) : min {min_}, max {max_} {}
// Static factory methods
static BoundingBox infinite() { return {}; }
static BoundingBox inverted()
{
return {{INFTY, INFTY, INFTY}, {-INFTY, -INFTY, -INFTY}};
}
inline BoundingBox operator&(const BoundingBox& other)
{
@ -35,29 +42,26 @@ struct BoundingBox {
// intersect operator
inline BoundingBox& operator&=(const BoundingBox& other)
{
xmin = std::max(xmin, other.xmin);
xmax = std::min(xmax, other.xmax);
ymin = std::max(ymin, other.ymin);
ymax = std::min(ymax, other.ymax);
zmin = std::max(zmin, other.zmin);
zmax = std::min(zmax, other.zmax);
min.x = std::max(min.x, other.min.x);
min.y = std::max(min.y, other.min.y);
min.z = std::max(min.z, other.min.z);
max.x = std::min(max.x, other.max.x);
max.y = std::min(max.y, other.max.y);
max.z = std::min(max.z, other.max.z);
return *this;
}
// union operator
inline BoundingBox& operator|=(const BoundingBox& other)
{
xmin = std::min(xmin, other.xmin);
xmax = std::max(xmax, other.xmax);
ymin = std::min(ymin, other.ymin);
ymax = std::max(ymax, other.ymax);
zmin = std::min(zmin, other.zmin);
zmax = std::max(zmax, other.zmax);
min.x = std::min(min.x, other.min.x);
min.y = std::min(min.y, other.min.y);
min.z = std::min(min.z, other.min.z);
max.x = std::max(max.x, other.max.x);
max.y = std::max(max.y, other.max.y);
max.z = std::max(max.z, other.max.z);
return *this;
}
inline Position min() const { return {xmin, ymin, zmin}; }
inline Position max() const { return {xmax, ymax, zmax}; }
};
} // namespace openmc

4
include/openmc/capi.h
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(

38
include/openmc/cell.h
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,18 @@ 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
//! \param[in] rot The rotation matrix of length 3 or 9
void set_rotation(const vector<double>& rot);
@ -312,12 +336,11 @@ public:
//----------------------------------------------------------------------------
// Data members
int32_t id_; //!< Unique ID
std::string name_; //!< User-defined name
Fill type_; //!< Material, universe, or lattice
int32_t universe_; //!< Universe # this cell is in
int32_t fill_; //!< Universe # filling this cell
int32_t n_instances_ {0}; //!< Number of instances of this cell
int32_t id_; //!< Unique ID
std::string name_; //!< User-defined name
Fill type_; //!< Material, universe, or lattice
int32_t universe_; //!< Universe # this cell is in
int32_t fill_; //!< Universe # filling this cell
//! \brief Index corresponding to this cell in distribcell arrays
int distribcell_index_ {C_NONE};
@ -333,6 +356,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_;

23
include/openmc/collision_track.h Normal file
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

25
include/openmc/constants.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
@ -63,6 +64,16 @@ constexpr int MAX_SAMPLE {100000};
// source region in the random ray solver
constexpr double MIN_HITS_PER_BATCH {1.5};
// The minimum flux value to be considered non-zero when computing adjoint
// sources. Positive values below this cutoff will be treated as zero, so as to
// 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
@ -281,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 };
@ -312,7 +323,10 @@ enum TallyScore {
SCORE_FISS_Q_PROMPT = -14, // prompt fission Q-value
SCORE_FISS_Q_RECOV = -15, // recoverable fission Q-value
SCORE_DECAY_RATE = -16, // delayed neutron precursor decay rate
SCORE_PULSE_HEIGHT = -17 // pulse-height
SCORE_PULSE_HEIGHT = -17, // pulse-height
SCORE_IFP_TIME_NUM = -18, // IFP lifetime numerator
SCORE_IFP_BETA_NUM = -19, // IFP delayed fraction numerator
SCORE_IFP_DENOM = -20 // IFP common denominator
};
// Global tally parameters
@ -322,6 +336,9 @@ enum class GlobalTally { K_COLLISION, K_ABSORPTION, K_TRACKLENGTH, LEAKAGE };
// Miscellaneous
constexpr int C_NONE {-1};
// Default value of generation for IFP
constexpr int DEFAULT_IFP_N_GENERATION {10};
// Interpolation rules
enum class Interpolation {
histogram = 1,

4
include/openmc/dagmc.h
View file

@ -20,7 +20,7 @@ void check_dagmc_root_univ();
} // namespace openmc
#ifdef DAGMC
#ifdef OPENMC_DAGMC_ENABLED
#include "DagMC.hpp"
#include "dagmcmetadata.hpp"
@ -216,6 +216,6 @@ int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ);
} // namespace openmc
#endif // DAGMC
#endif // OPENMC_DAGMC_ENABLED
#endif // OPENMC_DAGMC_H

188
include/openmc/distribution.h
View file

@ -15,6 +15,17 @@
namespace openmc {
//==============================================================================
// Helper function for computing importance weights from biased sampling
//==============================================================================
//! Compute importance weights for biased sampling
//! \param p Unnormalized original probability vector
//! \param b Unnormalized bias probability vector
//! \return Vector of importance weights (p_norm[i] / b_norm[i])
vector<double> compute_importance_weights(
const vector<double>& p, const vector<double>& b);
//==============================================================================
//! Abstract class representing a univariate probability distribution
//==============================================================================
@ -22,11 +33,41 @@ namespace openmc {
class Distribution {
public:
virtual ~Distribution() = default;
virtual double sample(uint64_t* seed) const = 0;
//! Sample a value from the distribution, handling biasing automatically
//! \param seed Pseudorandom number seed pointer
//! \return (sampled value, importance weight)
virtual std::pair<double, double> sample(uint64_t* seed) const;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
virtual double evaluate(double x) const;
//! Return integral of distribution
//! \return Integral of distribution
virtual double integral() const { return 1.0; };
//! Set bias distribution
virtual void set_bias(std::unique_ptr<Distribution> bias)
{
bias_ = std::move(bias);
}
const Distribution* bias() const { return bias_.get(); }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
virtual double sample_unbiased(uint64_t* seed) const = 0;
//! Read bias distribution from XML
//! \param node XML node that may contain a bias child element
void read_bias_from_xml(pugi::xml_node node);
// Biasing distribution
unique_ptr<Distribution> bias_;
};
using UPtrDist = unique_ptr<Distribution>;
@ -50,7 +91,7 @@ public:
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
//! \return Sampled index
size_t sample(uint64_t* seed) const;
// Properties
@ -67,7 +108,7 @@ private:
//! Normalize distribution so that probabilities sum to unity
void normalize();
//! Initialize alias tables for distribution
//! Initialize alias table for sampling
void init_alias();
};
@ -82,20 +123,30 @@ public:
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! \return (sampled value, sample weight)
std::pair<double, double> sample(uint64_t* seed) const override;
double integral() const override { return di_.integral(); };
//! Override set_bias as no-op (bias handled in constructor)
void set_bias(std::unique_ptr<Distribution> bias) override {}
// Properties
const vector<double>& x() const { return x_; }
const vector<double>& prob() const { return di_.prob(); }
const vector<size_t>& alias() const { return di_.alias(); }
const vector<double>& weight() const { return weight_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
vector<double> x_; //!< Possible outcomes
DiscreteIndex di_; //!< discrete probability distribution of
//!< outcome indices
vector<double> x_; //!< Possible outcomes
vector<double> weight_; //!< Importance weights (empty if unbiased)
DiscreteIndex di_; //!< Discrete probability distribution of outcome indices
};
//==============================================================================
@ -107,14 +158,20 @@ public:
explicit Uniform(pugi::xml_node node);
Uniform(double a, double b) : a_ {a}, b_ {b} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
double a() const { return a_; }
double b() const { return b_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
double a_; //!< Lower bound of distribution
double b_; //!< Upper bound of distribution
@ -131,15 +188,21 @@ public:
: offset_ {std::pow(a, n + 1)}, span_ {std::pow(b, n + 1) - offset_},
ninv_ {1 / (n + 1)} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
double a() const { return std::pow(offset_, ninv_); }
double b() const { return std::pow(offset_ + span_, ninv_); }
double n() const { return 1 / ninv_ - 1; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
//! Store processed values in object to allow for faster sampling
double offset_; //!< a^(n+1)
@ -156,13 +219,19 @@ public:
explicit Maxwell(pugi::xml_node node);
Maxwell(double theta) : theta_ {theta} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
double theta() const { return theta_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
double theta_; //!< Factor in exponential [eV]
};
@ -176,14 +245,20 @@ public:
explicit Watt(pugi::xml_node node);
Watt(double a, double b) : a_ {a}, b_ {b} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
double a() const { return a_; }
double b() const { return b_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
double a_; //!< Factor in exponential [eV]
double b_; //!< Factor in square root [1/eV]
@ -200,14 +275,20 @@ public:
Normal(double mean_value, double std_dev)
: mean_value_ {mean_value}, std_dev_ {std_dev} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
double mean_value() const { return mean_value_; }
double std_dev() const { return std_dev_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
double mean_value_; //!< middle of distribution [eV]
double std_dev_; //!< standard deviation [eV]
@ -223,10 +304,10 @@ public:
Tabular(const double* x, const double* p, int n, Interpolation interp,
const double* c = nullptr);
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
// properties
vector<double>& x() { return x_; }
@ -235,6 +316,12 @@ public:
Interpolation interp() const { return interp_; }
double integral() const override { return integral_; };
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
vector<double> x_; //!< tabulated independent variable
vector<double> p_; //!< tabulated probability density
@ -259,13 +346,19 @@ public:
explicit Equiprobable(pugi::xml_node node);
Equiprobable(const double* x, int n) : x_ {x, x + n} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
const vector<double>& x() const { return x_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
vector<double> x_; //! Possible outcomes
};
@ -280,18 +373,25 @@ public:
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! \return (sampled value, sample weight)
std::pair<double, double> sample(uint64_t* seed) const override;
double integral() const override { return integral_; }
private:
// Storrage for probability + distribution
using DistPair = std::pair<double, UPtrDist>;
//! Override set_bias as no-op (bias handled in constructor)
void set_bias(std::unique_ptr<Distribution> bias) override {}
vector<DistPair>
distribution_; //!< sub-distributions + cummulative probabilities
double integral_; //!< integral of distribution
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
vector<UPtrDist> distribution_; //!< Sub-distributions
vector<double> weight_; //!< Importance weights for component selection
DiscreteIndex di_; //!< Discrete probability distribution of indices
double integral_; //!< Integral of distribution
};
} // namespace openmc

45
include/openmc/distribution_multi.h
View file

@ -26,8 +26,8 @@ public:
//! Sample a direction from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Direction sampled
virtual Direction sample(uint64_t* seed) const = 0;
//! \return (sampled Direction, sample weight)
virtual std::pair<Direction, double> sample(uint64_t* seed) const = 0;
Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction
};
@ -43,14 +43,30 @@ public:
//! Sample a direction from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Direction sampled
Direction sample(uint64_t* seed) const override;
//! \return (sampled Direction, sample weight)
std::pair<Direction, double> sample(uint64_t* seed) const override;
//! Sample a direction and return evaluation of the PDF for biased sampling.
//! Note that bias distributions are intended to return unit-weight samples.
//! \param seed Pseudorandom number seed points
//! \return (sampled Direction, value of the PDF at this Direction)
std::pair<Direction, double> sample_as_bias(uint64_t* seed) const;
// Observing pointers
Distribution* mu() const { return mu_.get(); }
Distribution* phi() const { return phi_.get(); }
private:
//! Common sampling implementation
//! \param seed Pseudorandom number seed pointer
//! \param return_pdf If true, return PDF evaluation; if false, return
//! importance weight
//! \return (sampled Direction, weight or PDF value)
std::pair<Direction, double> sample_impl(
uint64_t* seed, bool return_pdf) const;
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
};
@ -64,11 +80,24 @@ Direction isotropic_direction(uint64_t* seed);
class Isotropic : public UnitSphereDistribution {
public:
Isotropic() {};
explicit Isotropic(pugi::xml_node node);
//! Sample a direction from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled direction
Direction sample(uint64_t* seed) const override;
//! \return (sampled direction, sample weight)
std::pair<Direction, double> sample(uint64_t* seed) const override;
// Set or get bias distribution
void set_bias(std::unique_ptr<PolarAzimuthal> bias)
{
bias_ = std::move(bias);
}
const PolarAzimuthal* bias() const { return bias_.get(); }
protected:
// Biasing distribution
unique_ptr<PolarAzimuthal> bias_;
};
//==============================================================================
@ -83,8 +112,8 @@ public:
//! Sample a direction from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled direction
Direction sample(uint64_t* seed) const override;
//! \return (sampled direction, sample weight)
std::pair<Direction, double> sample(uint64_t* seed) const override;
};
using UPtrAngle = unique_ptr<UnitSphereDistribution>;

37
include/openmc/distribution_spatial.h
View file

@ -19,7 +19,9 @@ public:
virtual ~SpatialDistribution() = default;
//! Sample a position from the distribution
virtual Position sample(uint64_t* seed) const = 0;
//! \param seed Pseudorandom number seed pointer
//! \return Sampled (position, importance weight)
virtual std::pair<Position, double> sample(uint64_t* seed) const = 0;
static unique_ptr<SpatialDistribution> create(pugi::xml_node node);
};
@ -34,8 +36,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
// Observer pointers
Distribution* x() const { return x_.get(); }
@ -58,8 +60,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
Distribution* r() const { return r_.get(); }
Distribution* phi() const { return phi_.get(); }
@ -83,8 +85,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
Distribution* r() const { return r_.get(); }
Distribution* cos_theta() const { return cos_theta_.get(); }
@ -109,8 +111,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
//! Sample the mesh for an element and position within that element
//! \param seed Pseudorandom number seed pointer
@ -133,8 +135,8 @@ public:
private:
int32_t mesh_idx_ {C_NONE};
DiscreteIndex elem_idx_dist_; //!< Distribution of
//!< mesh element indices
DiscreteIndex elem_idx_dist_; //!< Distribution of mesh element indices
vector<double> weight_; //!< Importance weights (empty if unbiased)
};
//==============================================================================
@ -149,12 +151,13 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
private:
std::vector<Position> point_cloud_;
DiscreteIndex point_idx_dist_; //!< Distribution of Position indices
vector<double> weight_; //!< Importance weights (empty if unbiased)
};
//==============================================================================
@ -167,8 +170,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
// Properties
bool only_fissionable() const { return only_fissionable_; }
@ -193,8 +196,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
Position r() const { return r_; }

18
include/openmc/geometry_aux.h
View file

@ -15,8 +15,6 @@
namespace openmc {
namespace model {
extern std::unordered_map<int32_t, std::unordered_map<int32_t, int32_t>>
universe_cell_counts;
extern std::unordered_map<int32_t, int32_t> universe_level_counts;
} // namespace model
@ -39,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
@ -80,15 +84,13 @@ void prepare_distribcell(
const std::vector<int32_t>* user_distribcells = nullptr);
//==============================================================================
//! Recursively search through the geometry and count cell instances.
//! Recursively search through the geometry and count universe instances.
//!
//! This function will update the Cell::n_instances value for each cell in the
//! geometry.
//! \param univ_indx The index of the universe to begin searching from (probably
//! the root universe).
//! This function will update Universe.n_instances_ for each
//! universe in the geometry.
//==============================================================================
void count_cell_instances(int32_t univ_indx);
void count_universe_instances();
//==============================================================================
//! Recursively search through universes and count universe instances.

10
include/openmc/hdf5_interface.h
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"
//==============================================================================

187
include/openmc/ifp.h Normal file
View file

@ -0,0 +1,187 @@
#ifndef OPENMC_IFP_H
#define OPENMC_IFP_H
#include "openmc/message_passing.h"
#include "openmc/particle.h"
#include "openmc/particle_data.h"
#include "openmc/settings.h"
namespace openmc {
//! Check the value of the IFP parameter for beta effective or both.
//!
//! \return true if "BetaEffective" or "Both", false otherwise.
bool is_beta_effective_or_both();
//! Check the value of the IFP parameter for generation time or both.
//!
//! \return true if "GenerationTime" or "Both", false otherwise.
bool is_generation_time_or_both();
//! Resize IFP vectors
//!
//! \param[in,out] delayed_groups List of delayed group numbers
//! \param[in,out] lifetimes List of lifetimes
//! \param[in] n Dimension to resize vectors
template<typename T, typename U>
void resize_ifp_data(vector<T>& delayed_groups, vector<U>& lifetimes, int64_t n)
{
if (is_beta_effective_or_both()) {
delayed_groups.resize(n);
}
if (is_generation_time_or_both()) {
lifetimes.resize(n);
}
}
//! Update a list of values by adding a new value if the size
//! of the list can accomodate the new value or by shifting all
//! values to the left (removing the first value of the list
//! and adding the new value at the end of the list).
//!
//! \param[in] value Value to add to the list
//! \param[in] data Initial version of the list
//! \return Updated list
template<typename T>
vector<T> _ifp(const T& value, const vector<T>& data)
{
vector<T> updated;
size_t source_idx = data.size();
if (source_idx < settings::ifp_n_generation) {
updated.resize(source_idx + 1);
for (size_t i = 0; i < source_idx; i++) {
updated[i] = data[i];
}
updated[source_idx] = value;
} else if (source_idx == settings::ifp_n_generation) {
updated.resize(source_idx);
for (size_t i = 0; i < source_idx - 1; i++) {
updated[i] = data[i + 1];
}
updated[source_idx - 1] = value;
}
return updated;
}
//! \brief Iterated Fission Probability (IFP) method.
//!
//! 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.
//!
//! \param[in] p Particle
//! \param[in] idx Bank index from the thread_safe_append call in physics.cpp
void ifp(const Particle& p, int64_t idx);
//! Resize the IFP banks used in the simulation
void resize_simulation_ifp_banks();
//! Retrieve IFP data from the IFP fission banks.
//!
//! \param[in] i_bank Index in the fission banks
//! \param[in,out] delayed_groups Delayed group numbers
//! \param[in,out] lifetimes Lifetimes lists
void copy_ifp_data_from_fission_banks(
int i_bank, vector<int>& delayed_groups, vector<double>& lifetimes);
#ifdef OPENMC_MPI
//! Deserialization information for transfer of IFP data using MPI
struct DeserializationInfo {
int64_t index_local; //!< local index
int64_t n; //!< number of sites sent
};
//! Broadcast the number of generation determined by the size of the first
//! element on the first processor.
//!
//! \param[in] n_generation Number of generations
//! \param[in] delayed_groups List of delayed group numbers lists
//! \param[in] lifetimes List of lifetimes lists
void broadcast_ifp_n_generation(int& n_generation,
const vector<vector<int>>& delayed_groups,
const vector<vector<double>>& lifetimes);
//! Send IFP data using MPI.
//!
//! \param[in] idx Index of the first site
//! \param[in] n Number of sites to send
//! \param[in] n_generation Number of generations
//! \param[in] neighbor Index of the neighboring processor
//! \param[in] requests MPI requests
//! \param[in] delayed_groups List of delayed group numbers lists
//! \param[out] send_delayed_groups Delayed group numbers buffer
//! \param[in] lifetimes List of lifetimes lists
//! \param[out] send_lifetimes Lifetimes buffer
void send_ifp_info(int64_t idx, int64_t n, int n_generation, int neighbor,
vector<MPI_Request>& requests, const vector<vector<int>>& delayed_groups,
vector<int>& send_delayed_groups, const vector<vector<double>>& lifetimes,
vector<double>& send_lifetimes);
//! Receive IFP data using MPI.
//!
//! \param[in] idx Index of the first site
//! \param[in] n Number of sites to receive
//! \param[in] n_generation Number of generations
//! \param[in] neighbor Index of the neighboring processor
//! \param[in] requests MPI requests
//! \param[in] delayed_groups List of delayed group numbers
//! \param[in] lifetimes List of lifetimes
//! \param[out] deserialization Information to deserialize the received data
void receive_ifp_data(int64_t idx, int64_t n, int n_generation, int neighbor,
vector<MPI_Request>& requests, vector<int>& delayed_groups,
vector<double>& lifetimes, vector<DeserializationInfo>& deserialization);
//! Copy partial IFP data from local lists to source banks.
//!
//! \param[in] idx Index of the first site
//! \param[in] n Number of sites to copy
//! \param[in] i_bank Index in the IFP source banks
//! \param[in] delayed_groups List of delayed group numbers lists
//! \param[in] lifetimes List of lifetimes lists
void copy_partial_ifp_data_to_source_banks(int64_t idx, int n, int64_t i_bank,
const vector<vector<int>>& delayed_groups,
const vector<vector<double>>& lifetimes);
//! Deserialize IFP information received using MPI and store it in
//! the IFP source banks.
//!
//! \param[in] n_generation Number of generations
//! \param[out] deserialization Information to deserialize the received data
//! \param[in] delayed_groups List of delayed group numbers
//! \param[in] lifetimes List of lifetimes
void deserialize_ifp_info(int n_generation,
const vector<DeserializationInfo>& deserialization,
const vector<int>& delayed_groups, const vector<double>& lifetimes);
#endif
//! Copy IFP temporary vectors to source banks.
//!
//! \param[in] delayed_groups List of delayed group numbers lists
//! \param[in] lifetimes List of lifetimes lists
void copy_complete_ifp_data_to_source_banks(
const vector<vector<int>>& delayed_groups,
const vector<vector<double>>& lifetimes);
//! Allocate temporary vectors for IFP data.
//!
//! \param[in,out] delayed_groups List of delayed group numbers lists
//! \param[in,out] lifetimes List of delayed group numbers lists
void allocate_temporary_vector_ifp(
vector<vector<int>>& delayed_groups, vector<vector<double>>& lifetimes);
//! Copy local IFP data to IFP fission banks.
//!
//! \param[in] delayed_groups_ptr Pointer to delayed group numbers
//! \param[in] lifetimes_ptr Pointer to lifetimes
void copy_ifp_data_to_fission_banks(
const vector<int>* delayed_groups_ptr, const vector<double>* lifetimes_ptr);
} // namespace openmc
#endif // OPENMC_IFP_H

2
include/openmc/lattice.h
View file

@ -76,7 +76,7 @@ public:
}
//! Populate the distribcell offset tables.
int32_t fill_offset_table(int32_t offset, int32_t target_univ_id, int map,
int32_t fill_offset_table(int32_t target_univ_id, int map,
std::unordered_map<int32_t, int32_t>& univ_count_memo);
//! \brief Check lattice indices.

12
include/openmc/material.h
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_; }
@ -107,6 +114,10 @@ public:
//! \return Density in [g/cm^3]
double density_gpcc() const { return density_gpcc_; }
//! Get charge density in [e/b-cm]
//! \return Charge density in [e/b-cm]
double charge_density() const { return charge_density_; };
//! Get name
//! \return Material name
const std::string& name() const { return name_; }
@ -177,6 +188,7 @@ public:
xt::xtensor<double, 1> atom_density_; //!< Nuclide atom density in [atom/b-cm]
double density_; //!< Total atom density in [atom/b-cm]
double density_gpcc_; //!< Total atom density in [g/cm^3]
double charge_density_; //!< Total charge density in [e/b-cm]
double volume_ {-1.0}; //!< Volume in [cm^3]
vector<bool> p0_; //!< Indicate which nuclides are to be treated with
//!< iso-in-lab scattering

11
include/openmc/math_functions.h
View file

@ -9,6 +9,7 @@
#include <cstdlib>
#include "openmc/position.h"
#include "openmc/search.h"
namespace openmc {
@ -200,5 +201,15 @@ std::complex<double> faddeeva(std::complex<double> z);
//! \return Derivative of Faddeeva function evaluated at z
std::complex<double> w_derivative(std::complex<double> z, int order);
//! Helper function to get index and interpolation function on an incident
//! energy grid
//!
//! \param energies energy grid
//! \param E incident energy
//! \param i grid index
//! \param f interpolation factor
void get_energy_index(
const vector<double>& energies, double E, int& i, double& f);
} // namespace openmc
#endif // OPENMC_MATH_FUNCTIONS_H

50
include/openmc/mcpl_interface.h
View file

@ -9,12 +9,6 @@
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
extern "C" const bool MCPL_ENABLED;
//==============================================================================
// Functions
//==============================================================================
@ -25,18 +19,46 @@ extern "C" const bool MCPL_ENABLED;
//! \return Vector of source sites
vector<SourceSite> mcpl_source_sites(std::string path);
//! Write an MCPL source file
//
//! Write an MCPL source file with stat:sum metadata
//!
//! This function writes particle data to an MCPL file. For MCPL >= 2.1.0,
//! it includes a stat:sum field (key: "openmc_np1") containing the total
//! number of source particles, which is essential for proper file merging
//! and weight normalization when using MCPL files with McStas/McXtrace.
//!
//! The stat:sum field follows the crash-safety pattern:
//! - Initially set to -1 when opening (indicates incomplete file)
//! - Updated with actual particle count before closing
//!
//! \param[in] filename Path to MCPL file
//! \param[in] source_bank Vector of SourceSites to write to file for this
//! MPI rank. Note that this can't be const due to
//! it being used as work space by MPI.
//! \param[in] bank_indx Pointer to vector of site index ranges over all
//! MPI ranks. This can be computed by calling
//! calculate_parallel_index_vector on
//! source_bank.size().
//! MPI rank.
//! \param[in] bank_index Pointer to vector of site index ranges over all
//! MPI ranks.
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();
//! Initialize the MCPL interface
void initialize_mcpl_interface_if_needed();
} // namespace openmc
#endif // OPENMC_MCPL_INTERFACE_H

93
include/openmc/mesh.h
View file

@ -19,14 +19,14 @@
#include "openmc/vector.h"
#include "openmc/xml_interface.h"
#ifdef DAGMC
#ifdef OPENMC_DAGMC_ENABLED
#include "moab/AdaptiveKDTree.hpp"
#include "moab/Core.hpp"
#include "moab/GeomUtil.hpp"
#include "moab/Matrix3.hpp"
#endif
#ifdef LIBMESH
#ifdef OPENMC_LIBMESH_ENABLED
#include "libmesh/bounding_box.h"
#include "libmesh/dof_map.h"
#include "libmesh/elem.h"
@ -61,7 +61,7 @@ extern vector<unique_ptr<Mesh>> meshes;
} // namespace model
#ifdef LIBMESH
#ifdef OPENMC_LIBMESH_ENABLED
namespace settings {
// used when creating new libMesh::MeshBase instances
extern unique_ptr<libMesh::LibMeshInit> libmesh_init;
@ -132,15 +132,18 @@ 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() {};
//! Update a position to the local coordinates of the mesh
virtual void local_coords(Position& r) const {};
//! Return a position in the local coordinates of the mesh
virtual Position local_coords(const Position& r) const { return r; };
@ -241,9 +244,7 @@ public:
//! \return Bounding box of mesh
BoundingBox bounding_box() const
{
auto ll = this->lower_left();
auto ur = this->upper_right();
return {ll.x, ur.x, ll.y, ur.y, ll.z, ur.z};
return {this->lower_left(), this->upper_right()};
}
virtual Position lower_left() const = 0;
@ -261,6 +262,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>;
@ -426,8 +428,7 @@ class PeriodicStructuredMesh : public StructuredMesh {
public:
PeriodicStructuredMesh() = default;
PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {};
void local_coords(Position& r) const override { r -= origin_; };
PeriodicStructuredMesh(hid_t group) : StructuredMesh {group} {};
Position local_coords(const Position& r) const override
{
@ -447,6 +448,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;
@ -486,6 +488,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
@ -497,6 +501,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;
@ -539,6 +544,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;
@ -603,6 +609,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;
@ -671,9 +678,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;
@ -773,13 +780,14 @@ private:
virtual void initialize() = 0;
};
#ifdef DAGMC
#ifdef OPENMC_DAGMC_ENABLED
class MOABMesh : public UnstructuredMesh {
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);
@ -943,12 +951,13 @@ private:
#endif
#ifdef LIBMESH
#ifdef OPENMC_LIBMESH_ENABLED
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);
@ -996,25 +1005,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>
@ -1028,8 +1038,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
@ -1048,6 +1084,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

1
include/openmc/message_passing.h
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

4
include/openmc/nuclide.h
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;

156
include/openmc/particle_data.h
View file

@ -8,7 +8,7 @@
#include "openmc/tallies/filter_match.h"
#include "openmc/vector.h"
#ifdef DAGMC
#ifdef OPENMC_DAGMC_ENABLED
#include "DagMC.hpp"
#endif
@ -47,7 +47,7 @@ struct SourceSite {
double time {0.0};
double wgt {1.0};
int delayed_group {0};
int surf_id {0};
int surf_id {SURFACE_NONE};
ParticleType particle;
// Extra attributes that don't show up in source written to 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]
@ -88,13 +107,37 @@ public:
//! clear data from a single coordinate level
void reset();
Position r; //!< particle position
Direction u; //!< particle direction
int cell {-1};
int universe {-1};
int lattice {-1};
array<int, 3> lattice_i {{-1, -1, -1}};
bool rotated {false}; //!< Is the level rotated?
// accessors
Position& r() { return r_; }
const Position& r() const { return r_; }
Direction& u() { return u_; }
const Direction& u() const { return u_; }
int& cell() { return cell_; }
const int& cell() const { return cell_; }
int& universe() { return universe_; }
const int& universe() const { return universe_; }
int& lattice() { return lattice_; }
int lattice() const { return lattice_; }
array<int, 3>& lattice_index() { return lattice_index_; }
const array<int, 3>& lattice_index() const { return lattice_index_; }
bool& rotated() { return rotated_; }
const bool& rotated() const { return rotated_; }
private:
// Data members
Position r_; //!< particle position
Direction u_; //!< particle direction
int cell_ {-1};
int universe_ {-1};
int lattice_ {-1};
array<int, 3> lattice_index_ {{-1, -1, -1}};
bool rotated_ {false}; //!< Is the level rotated?
};
//==============================================================================
@ -130,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
};
//==============================================================================
@ -186,23 +230,41 @@ struct CacheDataMG {
// Information about nearest boundary crossing
//==============================================================================
struct BoundaryInfo {
double distance {INFINITY}; //!< distance to nearest boundary
int surface {
SURFACE_NONE}; //!< surface token, non-zero if boundary is surface
int coord_level; //!< coordinate level after crossing boundary
array<int, 3>
lattice_translation {}; //!< which way lattice indices will change
class BoundaryInfo {
public:
void reset()
{
distance = INFINITY;
surface = SURFACE_NONE;
coord_level = 0;
lattice_translation = {0, 0, 0};
distance_ = INFINITY;
surface_ = SURFACE_NONE;
coord_level_ = 0;
lattice_translation_ = {0, 0, 0};
}
double& distance() { return distance_; }
const double& distance() const { return distance_; }
int& surface() { return surface_; }
const int& surface() const { return surface_; }
int coord_level() const { return coord_level_; }
int& coord_level() { return coord_level_; }
array<int, 3>& lattice_translation() { return lattice_translation_; }
const array<int, 3>& lattice_translation() const
{
return lattice_translation_;
}
// TODO: off-by-one
int surface_index() const { return std::abs(surface) - 1; }
int surface_index() const { return std::abs(surface()) - 1; }
private:
// Data members
double distance_ {INFINITY}; //!< distance to nearest boundary
int surface_ {
SURFACE_NONE}; //!< surface token, non-zero if boundary is surface
int coord_level_ {0}; //!< coordinate level after crossing boundary
array<int, 3> lattice_translation_ {
0, 0, 0}; //!< which way lattice indices will change
};
/*
@ -301,20 +363,20 @@ public:
const Position& u_last() const { return u_last_; }
// Accessors for position in global coordinates
Position& r() { return coord_[0].r; }
const Position& r() const { return coord_[0].r; }
Position& r() { return coord_[0].r(); }
const Position& r() const { return coord_[0].r(); }
// Accessors for position in local coordinates
Position& r_local() { return coord_[n_coord_ - 1].r; }
const Position& r_local() const { return coord_[n_coord_ - 1].r; }
Position& r_local() { return coord_[n_coord_ - 1].r(); }
const Position& r_local() const { return coord_[n_coord_ - 1].r(); }
// Accessors for direction in global coordinates
Direction& u() { return coord_[0].u; }
const Direction& u() const { return coord_[0].u; }
Direction& u() { return coord_[0].u(); }
const Direction& u() const { return coord_[0].u(); }
// Accessors for direction in local coordinates
Direction& u_local() { return coord_[n_coord_ - 1].u; }
const Direction& u_local() const { return coord_[n_coord_ - 1].u; }
Direction& u_local() { return coord_[n_coord_ - 1].u(); }
const Direction& u_local() const { return coord_[n_coord_ - 1].u(); }
// Surface token for the surface that the particle is currently on
int& surface() { return surface_; }
@ -330,7 +392,7 @@ public:
// Boundary information
BoundaryInfo& boundary() { return boundary_; }
#ifdef DAGMC
#ifdef OPENMC_DAGMC_ENABLED
// DagMC state variables
moab::DagMC::RayHistory& history() { return history_; }
Direction& last_dir() { return last_dir_; }
@ -347,6 +409,11 @@ public:
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
@ -375,7 +442,10 @@ private:
double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV
double sqrtkT_last_ {0.0}; //!< last temperature
#ifdef DAGMC
double density_mult_ {1.0}; //!< density multiplier
double density_mult_last_ {1.0}; //!< last density multiplier
#ifdef OPENMC_DAGMC_ENABLED
moab::DagMC::RayHistory history_;
Direction last_dir_;
#endif
@ -435,6 +505,7 @@ private:
double wgt_ {1.0};
double wgt_born_ {1.0};
double wgt_ww_born_ {-1.0};
double mu_;
double time_ {0.0};
double time_last_ {0.0};
@ -454,6 +525,9 @@ private:
int cell_born_ {-1};
// Iterated Fission Probability
double lifetime_ {0.0}; //!< neutron lifetime [s]
int n_collision_ {0};
bool write_track_ {false};
@ -542,6 +616,10 @@ public:
double& wgt_born() { return wgt_born_; }
double wgt_born() const { return wgt_born_; }
// Weight window value at birth
double& wgt_ww_born() { return wgt_ww_born_; }
const double& wgt_ww_born() const { return wgt_ww_born_; }
// Statistic weight of particle at last collision
double& wgt_last() { return wgt_last_; }
const double& wgt_last() const { return wgt_last_; }
@ -560,14 +638,20 @@ public:
double& time_last() { return time_last_; }
const double& time_last() const { return time_last_; }
// Particle lifetime
double& lifetime() { return lifetime_; }
const double& lifetime() const { return lifetime_; }
// What event took place, described in greater detail below
TallyEvent& event() { return event_; }
const TallyEvent& event() const { return event_; }
bool& fission() { return fission_; } // true if implicit fission
int& event_nuclide() { return event_nuclide_; } // index of collision nuclide
const int& event_nuclide() const { return event_nuclide_; }
int& event_mt() { return event_mt_; } // MT number of collision
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

6
include/openmc/photon.h
View file

@ -33,7 +33,6 @@ public:
int index_subshell; //!< index in SUBSHELLS
int threshold;
double n_electrons;
double binding_energy;
vector<Transition> transitions;
};
@ -90,6 +89,11 @@ public:
xt::xtensor<double, 1> binding_energy_;
xt::xtensor<double, 1> electron_pdf_;
// Map subshells from Compton profile data obtained from Biggs et al,
// "Hartree-Fock Compton profiles for the elements" to ENDF/B atomic
// relaxation data
xt::xtensor<int, 1> subshell_map_;
// Stopping power data
double I_; // mean excitation energy
xt::xtensor<int, 1> n_electrons_;

7
include/openmc/physics.h
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
//==============================================================================

2
include/openmc/plot.h
View file

@ -223,7 +223,7 @@ T SlicePlotBase::get_map() const
GeometryState p;
p.r() = xyz;
p.u() = dir;
p.coord(0).universe = model::root_universe;
p.coord(0).universe() = model::root_universe;
int level = slice_level_;
int j {};

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

@ -27,21 +27,22 @@ 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);
int64_t add_source_to_scalar_flux();
virtual void batch_reset();
void convert_source_regions_to_tallies();
void convert_source_regions_to_tallies(int64_t start_sr_id);
void reset_tally_volumes();
void random_ray_tally();
virtual void accumulate_iteration_flux();
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,10 +55,10 @@ 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
{
return source_regions_.n_source_regions();
@ -66,11 +67,18 @@ 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
static bool volume_normalized_flux_tallies_;
static bool adjoint_; // If the user wants outputs based on the adjoint flux
static double
diagonal_stabilization_rho_; // Adjusts strength of diagonal stabilization
// for transport corrected MGXS data
// Static variables to store source region meshes and domains
static std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>
@ -82,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
@ -106,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.
@ -127,16 +128,36 @@ public:
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
source_region_map_;
// 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, 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
// and random ray if used naively. This flag enables a stabilization
// technique.
bool is_transport_stabilization_needed_ {false};
protected:
//----------------------------------------------------------------------------
// Methods
void apply_external_source_to_source_region(
Discrete* discrete, double strength_factor, int64_t sr);
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);
@ -149,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

2
include/openmc/random_ray/linear_source_domain.h
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;

1
include/openmc/random_ray/random_ray.h
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
//----------------------------------------------------------------------------

9
include/openmc/random_ray/random_ray_simulation.h
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};

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

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

40
include/openmc/settings.h
View file

@ -24,6 +24,32 @@ enum class SSWCellType {
To,
};
// Type of IFP parameters
enum class IFPParameter {
None,
Both,
BetaEffective,
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
//==============================================================================
@ -33,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)?
@ -42,7 +69,8 @@ extern bool
delayed_photon_scaling; //!< Scale fission photon yield to include delayed
extern "C" bool entropy_on; //!< calculate Shannon entropy?
extern "C" bool
event_based; //!< use event-based mode (instead of history-based)
event_based; //!< use event-based mode (instead of history-based)
extern bool ifp_on; //!< Use IFP for kinetics parameters?
extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular?
extern bool material_cell_offsets; //!< create material cells offsets?
extern "C" bool output_summary; //!< write summary.h5?
@ -112,6 +140,10 @@ extern array<double, 4>
energy_cutoff; //!< Energy cutoff in [eV] for each particle type
extern array<double, 4>
time_cutoff; //!< Time cutoff in [s] for each particle type
extern int
ifp_n_generation; //!< Number of generation for Iterated Fission Probability
extern IFPParameter
ifp_parameter; //!< Parameter to calculate for Iterated Fission Probability
extern int
legendre_to_tabular_points; //!< number of points to convert Legendres
extern int max_order; //!< Maximum Legendre order for multigroup data
@ -132,9 +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

1
include/openmc/simulation.h
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?

36
include/openmc/source.h
View file

@ -21,10 +21,9 @@ namespace openmc {
// Constants
//==============================================================================
// Maximum number of external source spatial resamples to encounter before an
// error is thrown.
// Minimum number of external source sites rejected before checking againts the
// source_rejection_fraction
constexpr int EXTSRC_REJECT_THRESHOLD {10000};
constexpr double EXTSRC_REJECT_FRACTION {0.05};
//==============================================================================
// Global variables
@ -140,6 +139,9 @@ public:
DomainType domain_type() const { return domain_type_; }
const std::unordered_set<int32_t>& domain_ids() const { return domain_ids_; }
// Setter for spatial distribution
void set_space(UPtrSpace space) { space_ = std::move(space); }
protected:
// Indicates whether derived class already handles constraints
bool constraints_applied() const override { return true; }
@ -171,7 +173,7 @@ protected:
SourceSite sample(uint64_t* seed) const override;
private:
vector<SourceSite> sites_; //!< Source sites from a file
vector<SourceSite> sites_; //!< Source sites
};
//==============================================================================
@ -206,6 +208,23 @@ typedef unique_ptr<Source> create_compiled_source_t(std::string parameters);
//! Mesh-based source with different distributions for each element
//==============================================================================
// Helper class to sample spatial position on a single mesh element
class MeshElementSpatial : public SpatialDistribution {
public:
MeshElementSpatial(int32_t mesh_index, int elem_index)
: mesh_index_(mesh_index), elem_index_(elem_index)
{}
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return (sampled position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
private:
int32_t mesh_index_ {C_NONE}; //!< Index in global meshes array
int elem_index_; //! Index of mesh element
};
class MeshSource : public Source {
public:
// Constructors
@ -220,18 +239,15 @@ public:
double strength() const override { return space_->total_strength(); }
// Accessors
const std::unique_ptr<Source>& source(int32_t i) const
const unique_ptr<IndependentSource>& source(int32_t i) const
{
return sources_.size() == 1 ? sources_[0] : sources_[i];
}
protected:
bool constraints_applied() const override { return true; }
private:
// Data members
unique_ptr<MeshSpatial> space_; //!< Mesh spatial
vector<std::unique_ptr<Source>> sources_; //!< Source distributions
unique_ptr<MeshSpatial> space_; //!< Mesh spatial
vector<unique_ptr<IndependentSource>> sources_; //!< Source distributions
};
//==============================================================================

35
include/openmc/surface.h
View file

@ -2,6 +2,7 @@
#define OPENMC_SURFACE_H
#include <limits> // For numeric_limits
#include <set>
#include <string>
#include <unordered_map>
@ -378,7 +379,39 @@ public:
// Non-member functions
//==============================================================================
void read_surfaces(pugi::xml_node node);
//! Read surface definitions from XML and populate the global surfaces vector.
//!
//! This function parses surface elements from the XML input, creates the
//! appropriate surface objects, and identifies periodic surfaces along with
//! their albedo values and sense information.
//!
//! \param node XML node containing surface definitions
//! \param[out] periodic_pairs Set of surface ID pairs representing periodic
//! boundary conditions
//! \param[out] albedo_map Map of surface IDs to albedo values for periodic
//! surfaces
//! \param[out] periodic_sense_map Map of surface IDs to their sense values
//! (used to determine orientation for periodic BCs)
void read_surfaces(pugi::xml_node node,
std::set<std::pair<int, int>>& periodic_pairs,
std::unordered_map<int, double>& albedo_map,
std::unordered_map<int, int>& periodic_sense_map);
//! Resolve periodic surface pairs and assign boundary conditions.
//!
//! This function completes the setup of periodic boundary conditions by
//! resolving unpaired periodic surfaces, determining whether each pair
//! represents translational or rotational periodicity based on surface
//! normals, and assigning the appropriate boundary condition objects.
//!
//! \param[inout] periodic_pairs Set of surface ID pairs representing periodic
//! boundary conditions; unpaired entries are resolved
//! \param albedo_map Map of surface IDs to albedo values for periodic surfaces
//! \param periodic_sense_map Map of surface IDs to their sense values (used to
//! determine orientation for periodic BCs)
void prepare_boundary_conditions(std::set<std::pair<int, int>>& periodic_pairs,
std::unordered_map<int, double>& albedo_map,
std::unordered_map<int, int>& periodic_sense_map);
void free_memory_surfaces();

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

@ -33,6 +33,7 @@ enum class FilterType {
MATERIALFROM,
MESH,
MESHBORN,
MESH_MATERIAL,
MESH_SURFACE,
MU,
MUSURFACE,
@ -44,6 +45,7 @@ enum class FilterType {
SURFACE,
TIME,
UNIVERSE,
WEIGHT,
ZERNIKE,
ZERNIKE_RADIAL
};

14
include/openmc/tallies/filter_mesh.h
View file

@ -9,9 +9,9 @@
namespace openmc {
//==============================================================================
//! Indexes the location of particle events to a regular mesh. For tracklength
//! tallies, it will produce multiple valid bins and the bin weight will
//! correspond to the fraction of the track length that lies in that bin.
//! Indexes the location of particle events to a mesh. For tracklength tallies,
//! it will produce multiple valid bins and the bin weight will correspond to
//! the fraction of the track length that lies in that bin.
//==============================================================================
class MeshFilter : public Filter {
@ -51,6 +51,12 @@ public:
virtual bool translated() const { return translated_; }
virtual void set_rotation(const vector<double>& rotation);
virtual const vector<double>& rotation() const { return rotation_; }
virtual bool rotated() const { return rotated_; }
protected:
//----------------------------------------------------------------------------
// Data members
@ -58,6 +64,8 @@ protected:
int32_t mesh_; //!< Index of the mesh
bool translated_ {false}; //!< Whether or not the filter is translated
Position translation_ {0.0, 0.0, 0.0}; //!< Filter translation
bool rotated_ {false}; //!< Whether or not the filter is rotated
vector<double> rotation_; //!< Filter rotation
};
} // namespace openmc

114
include/openmc/tallies/filter_meshmaterial.h Normal file
View file

@ -0,0 +1,114 @@
#ifndef OPENMC_TALLIES_FILTER_MESHMATERIAL_H
#define OPENMC_TALLIES_FILTER_MESHMATERIAL_H
#include <cstdint>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include "openmc/position.h"
#include "openmc/random_ray/source_region.h"
#include "openmc/span.h"
#include "openmc/tallies/filter.h"
#include "openmc/vector.h"
namespace openmc {
//==============================================================================
//! Helper structs that define a combination of a mesh element index and a
//! material index and a functor for hashing to place in an unordered_map
//==============================================================================
struct ElementMat {
//! Check for equality
bool operator==(const ElementMat& other) const
{
return index_element == other.index_element && index_mat == other.index_mat;
}
int32_t index_element;
int32_t index_mat;
};
struct ElementMatHash {
std::size_t operator()(const ElementMat& k) const
{
size_t seed = 0;
hash_combine(seed, k.index_element);
hash_combine(seed, k.index_mat);
return seed;
}
};
//==============================================================================
//! Indexes the location of particle events to combinations of mesh element
//! index and material. For tracklength tallies, it will produce multiple valid
//! bins and the bin weight will correspond to the fraction of the track length
//! that lies in that bin.
//==============================================================================
class MeshMaterialFilter : public Filter {
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~MeshMaterialFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type_str() const override { return "meshmaterial"; }
FilterType type() const override { return FilterType::MESH_MATERIAL; }
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle& p, TallyEstimator estimator,
FilterMatch& match) const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
int32_t mesh() const { return mesh_; }
void set_mesh(int32_t mesh);
//! Set the bins based on a flat vector of alternating element index and
//! material IDs
void set_bins(span<int32_t> bins);
//! Set the bins based on a vector of (element, material index) pairs
void set_bins(vector<ElementMat>&& bins);
virtual void set_translation(const Position& translation);
virtual void set_translation(const double translation[3]);
virtual const Position& translation() const { return translation_; }
virtual bool translated() const { return translated_; }
private:
//----------------------------------------------------------------------------
// Data members
int32_t mesh_; //!< Index of the mesh
bool translated_ {false}; //!< Whether or not the filter is translated
Position translation_ {0.0, 0.0, 0.0}; //!< Filter translation
//! The indices of the mesh element-material combinations binned by this
//! filter.
vector<ElementMat> bins_;
//! The set of materials used in this filter
std::unordered_set<int32_t> materials_;
//! A map from mesh element-material indices to filter bin indices.
std::unordered_map<ElementMat, int32_t, ElementMatHash> map_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_MESHMATERIAL_H

51
include/openmc/tallies/filter_weight.h Normal file
View file

@ -0,0 +1,51 @@
#ifndef OPENMC_TALLIES_FILTER_WEIGHT_H
#define OPENMC_TALLIES_FILTER_WEIGHT_H
#include <string>
#include "openmc/span.h"
#include "openmc/tallies/filter.h"
#include "openmc/vector.h"
namespace openmc {
//==============================================================================
//! Bins the weights of the particles.
//==============================================================================
class WeightFilter : public Filter {
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~WeightFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type_str() const override { return "weight"; }
FilterType type() const override { return FilterType::WEIGHT; }
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle& p, TallyEstimator estimator,
FilterMatch& match) const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
const vector<double>& bins() const { return bins_; }
void set_bins(span<const double> bins);
protected:
//----------------------------------------------------------------------------
// Data members
vector<double> bins_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_WEIGHT_H

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

@ -106,6 +106,8 @@ public:
bool writable() const { return writable_; }
bool higher_moments() const { return higher_moments_; }
//----------------------------------------------------------------------------
// Other methods.
@ -169,10 +171,8 @@ public:
// We need to have quick access to some filters. The following gives indices
// for various filters that could be in the tally or C_NONE if they are not
// present.
int energy_filter_ {C_NONE};
int energyout_filter_ {C_NONE};
int delayedgroup_filter_ {C_NONE};
int cell_filter_ {C_NONE};
vector<Trigger> triggers_;
@ -192,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_;
};
@ -205,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 {
@ -241,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();

10
include/openmc/tallies/tally_scoring.h
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

3
include/openmc/universe.h
View file

@ -6,7 +6,7 @@
namespace openmc {
#ifdef DAGMC
#ifdef OPENMC_DAGMC_ENABLED
class DAGUniverse;
#endif
@ -29,6 +29,7 @@ class Universe {
public:
int32_t id_; //!< Unique ID
vector<int32_t> cells_; //!< Cells within this universe
int32_t n_instances_; //!< Number of instances of this universe
//! \brief Write universe information to an HDF5 group.
//! \param group_id An HDF5 group id.

8
include/openmc/urr.h
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(); }

1
include/openmc/weight_windows.h
View file

@ -71,6 +71,7 @@ struct WeightWindow {
{
lower_weight *= factor;
upper_weight *= factor;
survival_weight *= factor;
}
};

3
man/man1/openmc.1
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

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