Merge branch 'develop' into photonuclear-physics-pt1

AGENTS.md

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

CITATION.cff

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

docs/source/devguide/policies.rst

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

docs/source/io_formats/geometry.rst

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

docs/source/pythonapi/base.rst

@@ -176,7 +176,8 @@ Geometry Plotting
    :nosignatures:
    :template: myclass.rst
 
-   openmc.Plot
+   openmc.SlicePlot
+   openmc.VoxelPlot
    openmc.WireframeRayTracePlot
    openmc.SolidRayTracePlot
    openmc.Plots

docs/source/releasenotes/0.15.3.rst

@@ -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>`_)

docs/source/releasenotes/index.rst

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

docs/source/usersguide/geometry.rst

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

docs/source/usersguide/plots.rst

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

docs/source/usersguide/random_ray.rst

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

docs/source/usersguide/scripts.rst

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

docs/source/usersguide/variance_reduction.rst

@@ -51,7 +51,7 @@ With the :class:`~openmc.WeightWindowGenerator` instance added to the
 :attr:`~openmc.Settings`, the rest of the problem can be defined as normal. When
 running, note that the second iteration and beyond may be several orders of
 magnitude slower than the first. As the weight windows are applied in each
-iteration, particles may be agressively split, resulting in a large number of
+iteration, particles may be aggressively split, resulting in a large number of
 secondary (split) particles being generated per initial source particle. This is
 not necessarily a bad thing, as the split particles are much more efficient at
 exploring low flux regions of phase space as compared to initial particles.
@@ -161,7 +161,7 @@ solver, the Python input just needs to load the h5 file::
 
     settings.weight_window_checkpoints = {'collision': True, 'surface': True}
     settings.survival_biasing = False
-    settings.weight_windows = openmc.WeightWindowsList.from_hdf5('weight_windows.h5')
+    settings.weight_windows_file = "weight_windows.h5"
     settings.weight_windows_on = True
 
 The :class:`~openmc.WeightWindowGenerator` instance is not needed to load an

examples/lattice/hexagonal/build_xml.py

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

examples/lattice/nested/build_xml.py

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

examples/lattice/simple/build_xml.py

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

examples/pincell_random_ray/build_xml.py

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

include/openmc/bank_io.h

@@ -16,8 +16,9 @@
 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 banktype
+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
@@ -30,8 +31,8 @@ void write_bank_dataset(const char* dataset_name, hid_t group_id,
 #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, banktype, dspace, H5P_DEFAULT,
-    H5P_DEFAULT, H5P_DEFAULT);
+  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);
@@ -42,7 +43,7 @@ void write_bank_dataset(const char* dataset_name, hid_t group_id,
   hid_t plist = H5Pcreate(H5P_DATASET_XFER);
   H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
 
-  H5Dwrite(dset, banktype, memspace, dspace, plist, bank.data());
+  H5Dwrite(dset, membanktype, memspace, dspace, plist, bank.data());
 
   H5Sclose(dspace);
   H5Sclose(memspace);
@@ -52,7 +53,7 @@ void write_bank_dataset(const char* dataset_name, hid_t group_id,
   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, banktype, dspace,
+    hid_t dset = H5Dcreate(group_id, dataset_name, filebanktype, dspace,
       H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 
 #ifdef OPENMC_MPI
@@ -75,7 +76,8 @@ void write_bank_dataset(const char* dataset_name, hid_t group_id,
       H5Sselect_hyperslab(
         dspace_rank, H5S_SELECT_SET, start, nullptr, count, nullptr);
 
-      H5Dwrite(dset, banktype, memspace, dspace_rank, H5P_DEFAULT, bank.data());
+      H5Dwrite(
+        dset, membanktype, memspace, dspace_rank, H5P_DEFAULT, bank.data());
 
       H5Sclose(memspace);
       H5Sclose(dspace_rank);

include/openmc/boundary_condition.h

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

include/openmc/bounding_box.h

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

include/openmc/math_functions.h

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

include/openmc/mesh.h

@@ -244,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;

include/openmc/random_ray/flat_source_domain.h

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

include/openmc/tallies/filter_mesh.h

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

man/man1/openmc.1

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

openmc/__init__.py

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

openmc/_sparse_compat.py

@@ -0,0 +1,43 @@
+"""Compatibility module for scipy.sparse arrays
+
+This module provides a compatibility layer for working with scipy.sparse arrays
+across different scipy versions. Sparse arrays were introduced gradually in
+scipy, with full support arriving in scipy 1.15. This module provides a unified
+API that uses sparse arrays when available and falls back to sparse matrices for
+older scipy versions.
+
+For more information on the migration from sparse matrices to sparse arrays,
+see: https://docs.scipy.org/doc/scipy/reference/sparse.migration_to_sparray.html
+"""
+
+import scipy
+from scipy import sparse as sp
+
+# Check scipy version for feature availability
+_SCIPY_VERSION = tuple(map(int, scipy.__version__.split('.')[:2]))
+
+if _SCIPY_VERSION >= (1, 15):
+    # Use sparse arrays
+    csr_array = sp.csr_array
+    csc_array = sp.csc_array
+    dok_array = sp.dok_array
+    lil_array = sp.lil_array
+    eye_array = sp.eye_array
+    block_array = sp.block_array
+else:
+    # Fall back to sparse matrices
+    csr_array = sp.csr_matrix
+    csc_array = sp.csc_matrix
+    dok_array = sp.dok_matrix
+    lil_array = sp.lil_matrix
+    eye_array = sp.eye
+    block_array = sp.bmat
+
+__all__ = [
+    'csr_array',
+    'csc_array',
+    'dok_array',
+    'lil_array',
+    'eye_array',
+    'block_array',
+]

openmc/cmfd.py

@@ -25,6 +25,7 @@ import openmc.lib
 from .checkvalue import (check_type, check_length, check_value,
                          check_greater_than, check_less_than)
 from .exceptions import OpenMCError
+from ._sparse_compat import csr_array
 
 # See if mpi4py module can be imported, define have_mpi global variable
 try:
@@ -980,8 +981,7 @@ class CMFDRun:
         loss_row = self._loss_row
         loss_col = self._loss_col
         temp_data = np.ones(len(loss_row))
-        temp_loss = sparse.csr_matrix((temp_data, (loss_row, loss_col)),
-                                      shape=(n, n))
+        temp_loss = csr_array((temp_data, (loss_row, loss_col)), shape=(n, n))
         temp_loss.sort_indices()
 
         # Pass coremap as 1-d array of 32-bit integers
@@ -1585,7 +1585,7 @@ class CMFDRun:
         # Create csr matrix
         loss_row = self._loss_row
         loss_col = self._loss_col
-        loss = sparse.csr_matrix((data, (loss_row, loss_col)), shape=(n, n))
+        loss = csr_array((data, (loss_row, loss_col)), shape=(n, n))
         loss.sort_indices()
         return loss
 
@@ -1612,7 +1612,7 @@ class CMFDRun:
         # Create csr matrix
         prod_row = self._prod_row
         prod_col = self._prod_col
-        prod = sparse.csr_matrix((data, (prod_row, prod_col)), shape=(n, n))
+        prod = csr_array((data, (prod_row, prod_col)), shape=(n, n))
         prod.sort_indices()
         return prod
 

openmc/dagmc.py

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

openmc/data/decay.py

@@ -13,7 +13,7 @@ import openmc.checkvalue as cv
 from openmc.exceptions import DataError
 from openmc.mixin import EqualityMixin
 from openmc.stats import Discrete, Tabular, Univariate, combine_distributions
-from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
+from .data import ATOMIC_NUMBER, gnds_name
 from .function import INTERPOLATION_SCHEME
 from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
 
@@ -126,9 +126,7 @@ class FissionProductYields(EqualityMixin):
                 for j in range(n_products):
                     Z, A = divmod(int(values[4*j]), 1000)
                     isomeric_state = int(values[4*j + 1])
-                    name = ATOMIC_SYMBOL[Z] + str(A)
-                    if isomeric_state > 0:
-                        name += f'_m{isomeric_state}'
+                    name = gnds_name(Z, A, isomeric_state)
                     yield_j = ufloat(values[4*j + 2], values[4*j + 3])
                     yields[name] = yield_j
 
@@ -256,10 +254,7 @@ class DecayMode(EqualityMixin):
                         A += delta_A
                         Z += delta_Z
 
-        if self._daughter_state > 0:
-            return f'{ATOMIC_SYMBOL[Z]}{A}_m{self._daughter_state}'
-        else:
-            return f'{ATOMIC_SYMBOL[Z]}{A}'
+        return gnds_name(Z, A, self._daughter_state)
 
     @property
     def parent(self):
@@ -348,10 +343,7 @@ class Decay(EqualityMixin):
         self.nuclide['atomic_number'] = Z
         self.nuclide['mass_number'] = A
         self.nuclide['isomeric_state'] = metastable
-        if metastable > 0:
-            self.nuclide['name'] = f'{ATOMIC_SYMBOL[Z]}{A}_m{metastable}'
-        else:
-            self.nuclide['name'] = f'{ATOMIC_SYMBOL[Z]}{A}'
+        self.nuclide['name'] = gnds_name(Z, A, metastable)
         self.nuclide['mass'] = items[1]  # AWR
         self.nuclide['excited_state'] = items[2]  # State of the original nuclide
         self.nuclide['stable'] = (items[4] == 1)  # Nucleus stability flag

openmc/data/neutron.py

@@ -11,7 +11,7 @@ import h5py
 
 from . import HDF5_VERSION, HDF5_VERSION_MAJOR
 from .ace import Library, Table, get_table, get_metadata
-from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV
+from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV, gnds_name
 from .endf import (
     Evaluation, SUM_RULES, get_head_record, get_tab1_record, get_evaluations)
 from .fission_energy import FissionEnergyRelease
@@ -678,11 +678,7 @@ class IncidentNeutron(EqualityMixin):
         temperature = ev.target['temperature']
 
         # Determine name
-        element = ATOMIC_SYMBOL[atomic_number]
-        if metastable > 0:
-            name = f'{element}{mass_number}_m{metastable}'
-        else:
-            name = f'{element}{mass_number}'
+        name = gnds_name(atomic_number, mass_number, metastable)
 
         # Instantiate incident neutron data
         data = cls(name, atomic_number, mass_number, metastable,
@@ -769,6 +765,11 @@ class IncidentNeutron(EqualityMixin):
             for table in lib.tables[1:]:
                 data.add_temperature_from_ace(table)
 
+            # Use name based on ENDF evaluation. The name assigned by from_ace
+            # may be wrong for higher metastable states (e.g., Hf178_m2)
+            ev = evaluation if evaluation is not None else Evaluation(filename)
+            data.name = ev.gnds_name
+
             # Add 0K elastic scattering cross section
             if '0K' not in data.energy:
                 pendf = Evaluation(kwargs['pendf'])
@@ -779,7 +780,6 @@ class IncidentNeutron(EqualityMixin):
                 data[2].xs['0K'] = xs
 
             # Add fission energy release data
-            ev = evaluation if evaluation is not None else Evaluation(filename)
             if (1, 458) in ev.section:
                 data.fission_energy = f = FissionEnergyRelease.from_endf(ev, data)
             else:

openmc/data/njoy.py

@@ -15,19 +15,41 @@ import openmc.data
 # identifiers.
 ThermalTuple = namedtuple('ThermalTuple', ['name', 'zaids', 'nmix'])
 _THERMAL_DATA = {
+    'c_Ag': ThermalTuple('ag', [47107, 47109], 1),
     'c_Al27': ThermalTuple('al27', [13027], 1),
     'c_Al_in_Al2O3': ThermalTuple('asap00', [13027], 1),
+    'c_Al_in_Y3Al5O12': ThermalTuple('alyag', [13027], 1),
+    'c_Au': ThermalTuple('au', [79197], 1),
     'c_Be': ThermalTuple('be', [4009], 1),
     'c_Be_distinct': ThermalTuple('besd', [4009], 1),
     'c_Be_in_BeO': ThermalTuple('bebeo', [4009], 1),
     'c_Be_in_Be2C': ThermalTuple('bebe2c', [4009], 1),
     'c_Be_in_BeF2': ThermalTuple('bebef2', [4009], 1),
     'c_Be_in_FLiBe': ThermalTuple('beflib', [4009], 1),
+    'c_BeO': ThermalTuple('beo', [4009, 8016, 8017, 8018], 2),
+    'c_Bi': ThermalTuple('bi', [83209], 1),
+    'c_Bi_in_Ge3Bi4O12': ThermalTuple('bigbo', [83209], 1),
     'c_C6H6': ThermalTuple('benz', [1001, 6000, 6012], 2),
     'c_C_in_Be2C': ThermalTuple('cbe2c', [6000, 6012, 6013], 1),
     'c_C_in_C5O2H8': ThermalTuple('clucit', [6000, 6012, 6013], 1),
     'c_C_in_C8H8': ThermalTuple('cc8h8', [6000, 6012, 6013], 1),
+    'c_C_in_C19H16_liquid': ThermalTuple('c19liq', [6000, 6012, 6013], 1),
+    'c_C_in_C19H16_solid': ThermalTuple('c19sol', [6000, 6012, 6013], 1),
+    'c_C_in_C2H6O_liquid': ThermalTuple('ethliq', [6000, 6012, 6013], 1),
+    'c_C_in_C2H6O_solid': ThermalTuple('ethsol', [6000, 6012, 6013], 1),
+    'c_C_in_C6H6_liquid': ThermalTuple('benzlq', [6000, 6012, 6013], 1),
+    'c_C_in_C6H6_solid': ThermalTuple('benzsl', [6000, 6012, 6013], 1),
+    'c_C_in_C7H8_liquid': ThermalTuple('tolliq', [6000, 6012, 6013], 1),
+    'c_C_in_C7H8_solid': ThermalTuple('tolsol', [6000, 6012, 6013], 1),
+    'c_C_in_C8H10_liquid': ThermalTuple('xylliq', [6000, 6012, 6013], 1),
+    'c_C_in_C8H10_solid': ThermalTuple('xylsol', [6000, 6012, 6013], 1),
+    'c_C_in_C9H12_liquid': ThermalTuple('mesliq', [6000, 6012, 6013], 1),
+    'c_C_in_C9H12_solid': ThermalTuple('messol', [6000, 6012, 6013], 1),
     'c_C_in_CF2': ThermalTuple('ccf2', [6000, 6012, 6013], 1),
+    'c_C_in_CH2': ThermalTuple('cch2', [6000, 6012, 6013], 1),
+    'c_C_in_CH4_liquid': ThermalTuple('cch4lq', [6000, 6012, 6013], 1),
+    'c_C_in_CH4_solid': ThermalTuple('cch4sl', [6000, 6012, 6013], 1),
+    'c_C_in_Diamond': ThermalTuple('cdiam', [6000, 6012, 6013], 1),
     'c_C_in_SiC': ThermalTuple('csic', [6000, 6012, 6013], 1),
     'c_C_in_UC_100p': ThermalTuple('cuc100', [6000, 6012, 6013], 1),
     'c_C_in_UC_10p': ThermalTuple('cuc10', [6000, 6012, 6013], 1),
@@ -36,16 +58,29 @@ _THERMAL_DATA = {
     'c_C_in_UC_HALEU': ThermalTuple('cuchal', [6000, 6012, 6013], 1),
     'c_C_in_UC_HEU': ThermalTuple('cucheu', [6000, 6012, 6013], 1),
     'c_C_in_ZrC': ThermalTuple('czrc', [6000, 6012, 6013], 1),
+    'c_Ca': ThermalTuple('ca', [20040, 20042, 20043, 20044, 20046, 20048], 1),
     'c_Ca_in_CaH2': ThermalTuple('cacah2', [20040, 20042, 20043, 20044, 20046, 20048], 1),
+    'c_Ca_in_CaO2H2': ThermalTuple('cacaoh', [20040, 20042, 20043, 20044, 20046, 20048], 1),
+    'c_Cr': ThermalTuple('cr', [24050, 24052, 24053, 24054], 1),
+    'c_Cu': ThermalTuple('cu', [29063, 29065], 1),
     'c_D_in_7LiD': ThermalTuple('dlid', [1002], 1),
     'c_D_in_D2O': ThermalTuple('dd2o', [1002], 1),
     'c_D_in_D2O_solid': ThermalTuple('dice', [1002], 1),
+    'c_D_in_MgD2': ThermalTuple('dmgd2', [1002], 1),
     'c_F_in_Be2': ThermalTuple('fbef2', [9019], 1),
     'c_F_in_CF2': ThermalTuple('fcf2', [9019], 1),
     'c_F_in_FLiBe': ThermalTuple('fflibe', [9019], 1),
     'c_F_in_HF': ThermalTuple('f_hf', [9019], 1),
+    'c_F_in_LiF': ThermalTuple('flif', [9019], 1),
     'c_F_in_MgF2': ThermalTuple('fmgf2', [9019], 1),
     'c_Fe56': ThermalTuple('fe56', [26056], 1),
+    'c_Fe_in_Fe_alpha': ThermalTuple('fealph', [26054, 26056, 26057, 26058], 1),
+    'c_Fe_in_Fe_gamma': ThermalTuple('fegamm', [26054, 26056, 26057, 26058], 1),
+    'c_Ga_in_GaN': ThermalTuple('gagan', [31069, 31071], 1),
+    'c_Ga_in_GaSe': ThermalTuple('gagase', [31069, 31071], 1),
+    'c_Ge': ThermalTuple('ge', [32070, 32072, 32073, 32074, 32076], 1),
+    'c_Ge_in_Ge3Bi4O12': ThermalTuple('gegbo', [32070, 32072, 32073, 32074, 32076], 1),
+    'c_Ge_in_GeTe': ThermalTuple('gegete', [32070, 32072, 32073, 32074, 32076], 1),
     'c_Graphite': ThermalTuple('graph', [6000, 6012, 6013], 1),
     'c_Graphite_10p': ThermalTuple('grph10', [6000, 6012, 6013], 1),
     'c_Graphite_20p': ThermalTuple('grph20', [6000, 6012, 6013], 1),
@@ -54,7 +89,20 @@ _THERMAL_DATA = {
     'c_H_in_7LiH': ThermalTuple('hlih', [1001], 1),
     'c_H_in_C5O2H8': ThermalTuple('lucite', [1001], 1),
     'c_H_in_C8H8': ThermalTuple('hc8h8', [1001], 1),
+    'c_H_in_C19H16_liquid': ThermalTuple('h19liq', [1001], 1),
+    'c_H_in_C19H16_solid': ThermalTuple('h19sol', [1001], 1),
+    'c_H_in_C2H6O_liquid': ThermalTuple('hetliq', [1001], 1),
+    'c_H_in_C2H6O_solid': ThermalTuple('hetsol', [1001], 1),
+    'c_H_in_C6H6_liquid': ThermalTuple('hbzliq', [1001], 1),
+    'c_H_in_C6H6_solid': ThermalTuple('hbzsol', [1001], 1),
+    'c_H_in_C7H8_liquid': ThermalTuple('htlliq', [1001], 1),
+    'c_H_in_C7H8_solid': ThermalTuple('htlsol', [1001], 1),
+    'c_H_in_C8H10_liquid': ThermalTuple('hxyliq', [1001], 1),
+    'c_H_in_C8H10_solid': ThermalTuple('hxysol', [1001], 1),
+    'c_H_in_C9H12_liquid': ThermalTuple('hmsliq', [1001], 1),
+    'c_H_in_C9H12_solid': ThermalTuple('hmssol', [1001], 1),
     'c_H_in_CaH2': ThermalTuple('hcah2', [1001], 1),
+    'c_H_in_CaO2H2': ThermalTuple('hcaoh', [1001], 1),
     'c_H1_in_CaH2': ThermalTuple('h1cah2', [1001], 1),
     'c_H2_in_CaH2': ThermalTuple('h2cah2', [1001], 1),
     'c_H_in_CH2': ThermalTuple('hch2', [1001], 1),
@@ -64,32 +112,64 @@ _THERMAL_DATA = {
     'c_H_in_H2O': ThermalTuple('hh2o', [1001], 1),
     'c_H_in_H2O_solid': ThermalTuple('hice', [1001], 1),
     'c_H_in_HF': ThermalTuple('hhf', [1001], 1),
+    'c_H_in_KOH': ThermalTuple('hkoh', [1001], 1),
+    'c_H_in_LiH': ThermalTuple('hlih2', [1001], 1),
     'c_H_in_Mesitylene': ThermalTuple('mesi00', [1001], 1),
     'c_H_in_ParaffinicOil': ThermalTuple('hparaf', [1001], 1),
     'c_H_in_Toluene': ThermalTuple('tol00', [1001], 1),
+    'c_H_in_MgH2': ThermalTuple('hmgh2', [1001], 1),
+    'c_H_in_MgOH2': ThermalTuple('hmgoh', [1001], 1),
+    'c_H_in_NaMgH3': ThermalTuple('hnamg', [1001], 1),
+    'c_H_in_NaOH': ThermalTuple('hnaoh', [1001], 1),
+    'c_H_in_SrH2': ThermalTuple('hsrh2', [1001], 1),
     'c_H_in_UH3': ThermalTuple('huh3', [1001], 1),
     'c_H_in_YH2': ThermalTuple('hyh2', [1001], 1),
     'c_H_in_ZrH': ThermalTuple('hzrh', [1001], 1),
     'c_H_in_ZrH2': ThermalTuple('hzrh2', [1001], 1),
     'c_H_in_ZrHx': ThermalTuple('hzrhx', [1001], 1),
+    'c_I_in_NaI': ThermalTuple('inai', [53127], 1),
+    'c_K': ThermalTuple('k', [19039, 19040, 19041], 1),
+    'c_K_in_KOH': ThermalTuple('kkoh', [19039, 19040, 19041], 1),
     'c_Li_in_FLiBe': ThermalTuple('liflib', [3006, 3007], 1),
     'c_Li_in_7LiD': ThermalTuple('lilid', [3007], 1),
     'c_Li_in_7LiH': ThermalTuple('lilih', [3007], 1),
+    'c_Li_in_LiF': ThermalTuple('lilif', [3006, 3007], 1),
+    'c_Li_in_LiH': ThermalTuple('lilih2', [3006, 3007], 1),
     'c_Mg24': ThermalTuple('mg24', [12024], 1),
     'c_Mg_in_MgF2': ThermalTuple('mgmgf2', [12024, 12025, 12026], 1),
     'c_Mg_in_MgO': ThermalTuple('mgmgo', [12024, 12025, 12026], 1),
+    'c_Mg_in_MgD2': ThermalTuple('mgmgd2', [12024, 12025, 12026], 1),
+    'c_Mg_in_MgH2': ThermalTuple('mgmgh2', [12024, 12025, 12026], 1),
+    'c_Mg_in_MgOH2': ThermalTuple('mgoh2', [12024, 12025, 12026], 1),
+    'c_Mg_in_NaMgH3': ThermalTuple('mgnamg', [12024, 12025, 12026], 1),
+    'c_Mo': ThermalTuple('mo', [42092, 42094, 42095, 42096, 42097, 42098, 42100], 1),
+    'c_N_in_GaN': ThermalTuple('ngan', [7014, 7015], 1),
     'c_N_in_UN_100p': ThermalTuple('nun100', [7014, 7015], 1),
     'c_N_in_UN_10p': ThermalTuple('nun10', [7014, 7015], 1),
     'c_N_in_UN_5p': ThermalTuple('nun5', [7014, 7015], 1),
     'c_N_in_UN': ThermalTuple('n-un', [7014, 7015], 1),
     'c_N_in_UN_HALEU': ThermalTuple('nunhal', [7014, 7015], 1),
     'c_N_in_UN_HEU': ThermalTuple('nunheu', [7014, 7015], 1),
+    'c_Na': ThermalTuple('na', [11023], 1),
+    'c_Na_in_NaI': ThermalTuple('nanai', [11023], 1),
+    'c_Na_in_NaMgH3': ThermalTuple('nanamg', [11023], 1),
+    'c_Na_in_NaOH': ThermalTuple('nanaoh', [11023], 1),
+    'c_Nb': ThermalTuple('nb', [41093], 1),
+    'c_Ni': ThermalTuple('ni', [28058, 28060, 28061, 28062, 28064], 1),
     'c_O_in_Al2O3': ThermalTuple('osap00', [8016, 8017, 8018], 1),
     'c_O_in_BeO': ThermalTuple('obeo', [8016, 8017, 8018], 1),
     'c_O_in_C5O2H8': ThermalTuple('olucit', [8016, 8017, 8018], 1),
+    'c_O_in_C2H6O_liquid': ThermalTuple('oetliq', [8016, 8017, 8018], 1),
+    'c_O_in_C2H6O_solid': ThermalTuple('oetsol', [8016, 8017, 8018], 1),
+    'c_O_in_CaO2H2': ThermalTuple('ocaoh', [8016, 8017, 8018], 1),
     'c_O_in_D2O': ThermalTuple('od2o', [8016, 8017, 8018], 1),
     'c_O_in_H2O_solid': ThermalTuple('oice', [8016, 8017, 8018], 1),
     'c_O_in_MgO': ThermalTuple('omgo', [8016, 8017, 8018], 1),
+    'c_O_in_Ge3Bi4O12': ThermalTuple('ogbo', [8016, 8017, 8018], 1),
+    'c_O_in_H2O': ThermalTuple('oh2o', [8016, 8017, 8018], 1),
+    'c_O_in_KOH': ThermalTuple('okoh', [8016, 8017, 8018], 1),
+    'c_O_in_MgOH2': ThermalTuple('omgoh', [8016, 8017, 8018], 1),
+    'c_O_in_NaOH': ThermalTuple('onaoh', [8016, 8017, 8018], 1),
     'c_O_in_PuO2': ThermalTuple('opuo2', [8016, 8017, 8018], 1),
     'c_O_in_SiO2_alpha': ThermalTuple('osio2a', [8016, 8017, 8018], 1),
     'c_O_in_UO2_100p': ThermalTuple('ouo200', [8016, 8017, 8018], 1),
@@ -98,16 +178,26 @@ _THERMAL_DATA = {
     'c_O_in_UO2': ThermalTuple('ouo2', [8016, 8017, 8018], 1),
     'c_O_in_UO2_HALEU': ThermalTuple('ouo2hl', [8016, 8017, 8018], 1),
     'c_O_in_UO2_HEU': ThermalTuple('ouo2he', [8016, 8017, 8018], 1),
+    'c_O_in_Y3Al5O12': ThermalTuple('oyag', [8016, 8017, 8018], 1),
     'c_ortho_D': ThermalTuple('orthod', [1002], 1),
     'c_ortho_H': ThermalTuple('orthoh', [1001], 1),
     'c_para_D': ThermalTuple('parad', [1002], 1),
     'c_para_H': ThermalTuple('parah', [1001], 1),
+    'c_Pb': ThermalTuple('pb', [82204, 82206, 82207, 82208], 1),
+    'c_Pd': ThermalTuple('pd', [46102, 46104, 46105, 46106, 46108, 46110], 1),
+    'c_Pt': ThermalTuple('pt', [78190, 78192, 78194, 78195, 78196, 78198], 1),
     'c_Pu_in_PuO2': ThermalTuple('puo2', [94239, 94240, 94241, 94242, 94243], 1),
     'c_Si28': ThermalTuple('si00', [14028], 1),
     'c_Si_in_SiC': ThermalTuple('sisic', [14028, 14029, 14030], 1),
     'c_Si_in_SiO2_alpha': ThermalTuple('si_o2a', [14028, 14029, 14030], 1),
     'c_SiO2_alpha': ThermalTuple('sio2-a', [8016, 8017, 8018, 14028, 14029, 14030], 3),
     'c_SiO2_beta': ThermalTuple('sio2-b', [8016, 8017, 8018, 14028, 14029, 14030], 3),
+    'c_S_in_ZnS': ThermalTuple('szns', [16032, 16033, 16034, 16036], 1),
+    'c_Se_in_GaSe': ThermalTuple('segase', [34074, 34076, 34077, 34078, 34080, 34082], 1),
+    'c_Sn': ThermalTuple('sn', [50112, 50114, 50115, 50116, 50117, 50118, 50119, 50120, 50122, 50124], 1),
+    'c_Sr_in_SrH2': ThermalTuple('srsrh2', [38084, 38086, 38087, 38088], 1),
+    'c_Te_in_GeTe': ThermalTuple('tegete', [52120, 52122, 52123, 52124, 52125, 52126, 52128, 52130], 1),
+    'c_Ti': ThermalTuple('ti', [22046, 22047, 22048, 22049, 22050], 1),
     'c_U_metal_100p': ThermalTuple('u-100p', [92233, 92234, 92235, 92236, 92238], 1),
     'c_U_metal_10p': ThermalTuple('u-10p', [92233, 92234, 92235, 92236, 92238], 1),
     'c_U_metal_5p': ThermalTuple('u-5p', [92233, 92234, 92235, 92236, 92238], 1),
@@ -132,7 +222,13 @@ _THERMAL_DATA = {
     'c_U_in_UO2': ThermalTuple('uuo2', [92233, 92234, 92235, 92236, 92238], 1),
     'c_U_in_UO2_HALEU': ThermalTuple('uo2hal', [92233, 92234, 92235, 92236, 92238], 1),
     'c_U_in_UO2_HEU': ThermalTuple('uo2heu', [92233, 92234, 92235, 92236, 92238], 1),
+    'c_V': ThermalTuple('v', [23050, 23051], 1),
+    'c_W': ThermalTuple('w', [74180, 74182, 74183, 74184, 74186], 1),
+    'c_Y_in_Y3Al5O12': ThermalTuple('yyag', [39089], 1),
     'c_Y_in_YH2': ThermalTuple('yyh2', [39089], 1),
+    'c_Zn': ThermalTuple('zn', [30064, 30066, 30067, 30068, 30070], 1),
+    'c_Zn_in_ZnS': ThermalTuple('znzns', [30064, 30066, 30067, 30068, 30070], 1),
+    'c_Zr': ThermalTuple('zr', [40090, 40091, 40092, 40094, 40096], 1),
     'c_Zr_in_ZrC': ThermalTuple('zrzrc', [40000, 40090, 40091, 40092, 40094, 40096], 1),
     'c_Zr_in_ZrH': ThermalTuple('zrzrh', [40000, 40090, 40091, 40092, 40094, 40096], 1),
     'c_Zr_in_ZrH2': ThermalTuple('zrzrh2', [40000, 40090, 40091, 40092, 40094, 40096], 1),
@@ -572,7 +668,7 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
     else:
         with warnings.catch_warnings(record=True) as w:
             proper_name = openmc.data.get_thermal_name(zsymam_thermal)
-            if w:
+            if w or proper_name not in _THERMAL_DATA:
                 raise RuntimeError(
                     f"Thermal scattering material {zsymam_thermal} not "
                     "recognized. Please contact OpenMC developers at "

openmc/data/thermal.py

@@ -27,20 +27,41 @@ from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE,
 
 
 _THERMAL_NAMES = {
+    'c_Ag': ('ag',),
     'c_Al27': ('al', 'al27', 'al-27', '13-al- 27'),
     'c_Al_in_Al2O3': ('asap00', 'asap', 'al(al2o3)'),
-    'c_Be': ('be', 'be-metal', 'be-met', 'be00', 'be-metal', 'be metal', '4-be'),
+    'c_Al_in_Y3Al5O12': ('al(y3al5o1', 'alyag'),
+    'c_Au': ('au',),
+    'c_Be': ('be', 'be-metal', 'be-met', 'be00', 'be-metal', 'be metal', '4-be', '4-be-'),
     'c_BeO': ('beo',),
     'c_Be_distinct': ('besd', 'be+sd'),
     'c_Be_in_BeO': ('bebeo', 'be-beo', 'be-o', 'be/o', 'bbeo00', 'be(beo)', 'be_beo'),
     'c_Be_in_Be2C': ('bebe2c', 'be(be2c)'),
     'c_Be_in_BeF2': ('bebef2', 'be in bef2'),
     'c_Be_in_FLiBe': ('beflib', 'be(flibe)'),
+    'c_Bi': ('83-bi-', 'bi'),
+    'c_Bi_in_Ge3Bi4O12': ('bi(ge3bi4o', 'bigbo'),
     'c_C6H6': ('benz', 'c6h6', 'benzine'),
     'c_C_in_Be2C': ('cbe2c', 'c(be2c)'),
+    'c_C_in_C19H16_liquid': ('c(c19h16)l', 'c19liq'),
+    'c_C_in_C19H16_solid': ('c(c19h16)s', 'c19sol'),
+    'c_C_in_C2H6O_liquid': ('c(c2h6o)l', 'ethliq'),
+    'c_C_in_C2H6O_solid': ('c(c2h6o)s', 'ethsol'),
     'c_C_in_C5O2H8': ('clucit', 'c(lucite)'),
+    'c_C_in_C6H6_liquid': ('c(c6h6)l', 'benzlq'),
+    'c_C_in_C6H6_solid': ('c(c6h6)s', 'benzsl'),
+    'c_C_in_C7H8_liquid': ('c(c7h8)l', 'tolliq'),
+    'c_C_in_C7H8_solid': ('c(c7h8)s', 'tolsol'),
     'c_C_in_C8H8': ('cc8h8', 'c(polystyr'),
+    'c_C_in_C8H10_liquid': ('c(m-c8h10)l', 'xylliq'),
+    'c_C_in_C8H10_solid': ('c(m-c8h10)s', 'xylsol'),
+    'c_C_in_C9H12_liquid': ('c(c9h12)l', 'mesliq'),
+    'c_C_in_C9H12_solid': ('c(c9h12)s', 'messol'),
     'c_C_in_CF2': ('ccf2', 'c(teflon)'),
+    'c_C_in_CH2': ('c(c2h4)n r', 'cch2'),
+    'c_C_in_CH4_liquid': ('c(ch4)l', 'cch4lq'),
+    'c_C_in_CH4_solid': ('c(ch4)s', 'cch4sl'),
+    'c_C_in_Diamond': ('c(c-diamon', 'cdiam'),
     'c_C_in_SiC': ('csic', 'c-sic', 'c(3c-sic)', 'c_sic'),
     'c_C_in_UC_100p': ('cuc100', 'cinuc_100p'),
     'c_C_in_UC_10p': ('cuc10', 'cinuc_10p'),
@@ -49,16 +70,29 @@ _THERMAL_NAMES = {
     'c_C_in_UC_HALEU': ('cuchal', 'cinuc_haleu'),
     'c_C_in_UC_HEU': ('cucheu', 'cinuc_heu'),
     'c_C_in_ZrC': ('czrc', 'c(zrc)'),
+    'c_Ca': ('ca',),
     'c_Ca_in_CaH2': ('cah', 'cah00', 'cacah2', 'ca(cah2)', 'ca_cah2'),
+    'c_Ca_in_CaO2H2': ('ca(caoh2)', 'cacaoh'),
+    'c_Cr': ('cr',),
+    'c_Cu': ('cu',),
     'c_D_in_7LiD': ('dlid', 'd(7lid)'),
     'c_D_in_D2O': ('dd2o', 'd-d2o', 'hwtr', 'hw', 'dhw00', 'd(d2o)'),
     'c_D_in_D2O_solid': ('dice',),
+    'c_D_in_MgD2': ('d(mgd2)', 'dmgd2'),
     'c_F_in_Be2': ('fbef2', 'f in bef2'),
     'c_F_in_CF2': ('fcf2', 'f(teflon)'),
     'c_F_in_FLiBe': ('fflibe', 'f(flibe)'),
     'c_F_in_HF': ('f_hf',),
+    'c_F_in_LiF': ('f(lif)', 'flif'),
     'c_F_in_MgF2': ('fmgf2', 'f in mgf2'),
     'c_Fe56': ('fe', 'fe56', 'fe-56', '26-fe- 56'),
+    'c_Fe_in_Fe_alpha': ('fe(fe-alph', 'fealph'),
+    'c_Fe_in_Fe_gamma': ('fe(fe-gamm', 'fegamm'),
+    'c_Ga_in_GaN': ('ga(gan)', 'gagan'),
+    'c_Ga_in_GaSe': ('ga(gase)', 'gagase'),
+    'c_Ge': ('ge',),
+    'c_Ge_in_Ge3Bi4O12': ('ge(ge3bi4o', 'gegbo'),
+    'c_Ge_in_GeTe': ('ge(gete)', 'gegete'),
     'c_Graphite': ('graph', 'grph', 'gr', 'gr00', 'graphite'),
     'c_Graphite_10p': ('grph10', '10p graphit'),
     'c_Graphite_20p': ('grph20', '20 graphite'),
@@ -68,41 +102,86 @@ _THERMAL_NAMES = {
     'c_H_in_C5O2H8': ('lucite', 'c5o2h8', 'h-luci', 'h(lucite)'),
     'c_H_in_C8H8': ('hc8h8', 'h(polystyr'),
     'c_H_in_CaH2': ('hcah2', 'hca00', 'h(cah2)'),
+    'c_H_in_CaO2H2': ('h(caoh2)', 'hcaoh'),
     'c_H1_in_CaH2': ('h1cah2', 'h1_cah2'),
     'c_H2_in_CaH2': ('h2cah2', 'h2_cah2'),
-    'c_H_in_CH2': ('hch2', 'poly', 'pol', 'h-poly', 'pol00', 'h(ch2)'),
-    'c_H_in_CH4_liquid': ('lch4', 'lmeth', 'l-ch4'),
-    'c_H_in_CH4_solid': ('sch4', 'smeth', 's-ch4'),
+    'c_H_in_C19H16_liquid': ('h(c19h16)l', 'h19liq'),
+    'c_H_in_C19H16_solid': ('h(c19h16)s', 'h19sol'),
+    'c_H_in_C2H6O_liquid': ('h(c2h6o)l', 'hetliq'),
+    'c_H_in_C2H6O_solid': ('h(c2h6o)s', 'hetsol'),
+    'c_H_in_C6H6_liquid': ('h(c6h6)l', 'hbzliq'),
+    'c_H_in_C6H6_solid': ('h(c6h6)s', 'hbzsol'),
+    'c_H_in_C7H8_liquid': ('h(c7h8)l', 'htlliq'),
+    'c_H_in_C7H8_solid': ('h(c7h8)s', 'htlsol'),
+    'c_H_in_C8H10_liquid': ('h(m-c8h10)l', 'hxyliq'),
+    'c_H_in_C8H10_solid': ('h(m-c8h10)s', 'hxysol'),
+    'c_H_in_C9H12_liquid': ('h(c9h12)l', 'hmsliq'),
+    'c_H_in_C9H12_solid': ('h(c9h12)s', 'hmssol'),
+    'c_H_in_CH2': ('hch2', 'poly', 'pol', 'h-poly', 'pol00', 'h(ch2)', 'h(c2h4)n r'),
+    'c_H_in_CH4_liquid': ('lch4', 'lmeth', 'l-ch4', 'h(ch4)l'),
+    'c_H_in_CH4_solid': ('sch4', 'smeth', 's-ch4', 'h(ch4)s'),
     'c_H_in_CH4_solid_phase_II': ('sch4p2',),
     'c_H_in_H2O': ('hh2o', 'h-h2o', 'lwtr', 'lw', 'lw00', 'h(h2o)'),
     'c_H_in_H2O_solid': ('hice', 'h-ice', 'ice00', 'h(ice-ih)', 'h(ice)'),
     'c_H_in_HF': ('hhf', 'h(hf)', 'h_hf'),
+    'c_H_in_KOH': ('h(koh)', 'hkoh'),
+    'c_H_in_LiH': ('h(lih)', 'hlih2'),
     'c_H_in_Mesitylene': ('mesi00', 'mesi', 'mesi-phii'),
+    'c_H_in_MgH2': ('h(mgh2)', 'hmgh2'),
+    'c_H_in_MgOH2': ('h(mgoh2)', 'hmgoh'),
+    'c_H_in_NaMgH3': ('h(namgh3)', 'hnamg'),
+    'c_H_in_NaOH': ('h(naoh)', 'hnaoh'),
     'c_H_in_ParaffinicOil': ('hparaf', 'h(paraffin', 'h(paraffini'),
+    'c_H_in_SrH2': ('h(srh2)', 'hsrh2'),
     'c_H_in_Toluene': ('tol00', 'tol', 'tolue-phii'),
     'c_H_in_UH3': ('huh3', 'h(uh3)'),
     'c_H_in_YH2': ('hyh2', 'h-yh2', 'h(yh2)'),
     'c_H_in_ZrH': ('hzrh', 'h-zrh', 'h-zr', 'h/zr', 'hzr', 'hzr00', 'h(zrh)'),
-    'c_H_in_ZrH2': ('hzrh2', 'h(zrh2)'),
-    'c_H_in_ZrHx': ('hzrhx', 'h(zrhx)'),
+    'c_H_in_ZrH2': ('hzrh2', 'h(zrh2)', 'h(zrh2) in'),
+    'c_H_in_ZrHx': ('hzrhx', 'h(zrhx)', 'h(zrh15) i'),
+    'c_I_in_NaI': ('i(nai)', 'inai'),
+    'c_K': ('k',),
+    'c_K_in_KOH': ('k(koh)', 'kkoh'),
     'c_Li_in_FLiBe': ('liflib', 'li(flibe)'),
     'c_Li_in_7LiD': ('lilid', '7li(7lid)'),
     'c_Li_in_7LiH': ('lilih', '7li(7lih)'),
+    'c_Li_in_LiF': ('li(lif)', 'lilif'),
+    'c_Li_in_LiH': ('li(lih)', 'lilih2'),
     'c_Mg24': ('mg', 'mg24', 'mg00', '24-mg'),
+    'c_Mg_in_MgD2': ('mg(mgd2)', 'mgmgd2'),
     'c_Mg_in_MgF2': ('mgmgf2', 'mg in mgf2'),
+    'c_Mg_in_MgH2': ('mg(mgh2)', 'mgmgh2'),
     'c_Mg_in_MgO': ('mgmgo', 'mg in mgo'),
+    'c_Mg_in_MgOH2': ('mg(mgoh2)', 'mgoh2'),
+    'c_Mg_in_NaMgH3': ('mg(namgh3)', 'mgnamg'),
+    'c_Mo': ('mo',),
+    'c_N_in_GaN': ('n(gan)', 'ngan'),
     'c_N_in_UN_100p': ('nun100', 'n-un-100p'),
     'c_N_in_UN_10p': ('nun10', 'n-un-10p'),
     'c_N_in_UN_5p': ('nun5', 'n-un-5p'),
     'c_N_in_UN': ('n-un', 'n(un)', 'n(un)    l', 'ninun'),
     'c_N_in_UN_HALEU': ('nunhal', 'n-un-haleu'),
     'c_N_in_UN_HEU': ('nunheu', 'n-un-heu'),
+    'c_Na': ('na',),
+    'c_Na_in_NaI': ('na(nai)', 'nanai'),
+    'c_Na_in_NaMgH3': ('na(namgh3)', 'nanamg'),
+    'c_Na_in_NaOH': ('na(naoh)', 'nanaoh'),
+    'c_Nb': ('nb',),
+    'c_Ni': ('ni',),
     'c_O_in_Al2O3': ('osap00', 'osap', 'o(al2o3)'),
     'c_O_in_BeO': ('obeo', 'o-beo', 'o-be', 'o/be', 'obeo00', 'o(beo)', 'o_beo'),
+    'c_O_in_C2H6O_liquid': ('o(c2h6o)l', 'oetliq'),
+    'c_O_in_C2H6O_solid': ('o(c2h6o)s', 'oetsol'),
     'c_O_in_C5O2H8': ('olucit', 'o(lucite)'),
+    'c_O_in_CaO2H2': ('o(caoh2)', 'ocaoh'),
     'c_O_in_D2O': ('od2o', 'o-d2o', 'ohw00', 'o(d2o)'),
+    'c_O_in_H2O': ('o(h2o)', 'oh2o'),
     'c_O_in_H2O_solid': ('oice', 'o-ice', 'o(ice-ih)'),
+    'c_O_in_Ge3Bi4O12': ('o(ge3bi4o1', 'ogbo'),
+    'c_O_in_KOH': ('o(koh)', 'okoh'),
     'c_O_in_MgO': ('omgo', 'o in mgo'),
+    'c_O_in_MgOH2': ('o(mgoh2)', 'omgoh'),
+    'c_O_in_NaOH': ('o(naoh)', 'onaoh'),
     'c_O_in_PuO2': ('opuo2', 'o in puo2'),
     'c_O_in_SiO2_alpha': ('osio2a', 'o_sio2a'),
     'c_O_in_UO2_100p': ('ouo200', 'o-uo2-100p'),
@@ -111,16 +190,26 @@ _THERMAL_NAMES = {
     'c_O_in_UO2': ('ouo2', 'o-uo2', 'o2-u', 'o2/u', 'ouo200', 'o(uo2)'),
     'c_O_in_UO2_HALEU': ('ouo2hl', 'ouo2-haleu'),
     'c_O_in_UO2_HEU': ('ouo2he', 'o_uo2-heu'),
+    'c_O_in_Y3Al5O12': ('o(y3al5o12', 'oyag'),
     'c_ortho_D': ('orthod', 'orthoD', 'dortho', 'od200', 'ortod', 'ortho-d'),
     'c_ortho_H': ('orthoh', 'orthoH', 'hortho', 'oh200', 'ortoh', 'ortho-h'),
     'c_para_D': ('parad', 'paraD', 'dpara', 'pd200', 'para-d'),
     'c_para_H': ('parah', 'paraH', 'hpara', 'ph200', 'para-h'),
+    'c_Pb': ('pb',),
+    'c_Pd': ('pd',),
+    'c_Pt': ('pt',),
     'c_Pu_in_PuO2': ('puo2', 'pu in puo2'),
     'c_Si28': ('si00', 'sili', 'si'),
     'c_Si_in_SiC': ('sisic', 'si-sic', 'si(3c-sic)', 'si_sic'),
     'c_Si_in_SiO2_alpha': ('si_o2a', 'si_sio2a'),
-    'c_SiO2_alpha': ('sio2', 'sio2a', 'sio2alpha'),
-    'c_SiO2_beta': ('sio2b', 'sio2beta'),
+    'c_SiO2_alpha': ('sio2', 'sio2a', 'sio2alpha', 'sio2-a'),
+    'c_SiO2_beta': ('sio2b', 'sio2beta', 'sio2-b'),
+    'c_S_in_ZnS': ('s(zns-spha', 'szns'),
+    'c_Se_in_GaSe': ('se(gase)', 'segase'),
+    'c_Sn': ('sn',),
+    'c_Sr_in_SrH2': ('sr(srh2)', 'srsrh2'),
+    'c_Te_in_GeTe': ('te(gete)', 'tegete'),
+    'c_Ti': ('ti',),
     'c_U_metal_100p': ('u-100p',),
     'c_U_metal_10p': ('u-10p',),
     'c_U_metal_5p': ('u-5p',),
@@ -145,11 +234,17 @@ _THERMAL_NAMES = {
     'c_U_in_UO2': ('uuo2', 'u-uo2', 'u-o2', 'u/o2', 'uuo200', 'u(uo2)'),
     'c_U_in_UO2_HALEU': ('uo2hal', 'uuo2-haleu'),
     'c_U_in_UO2_HEU': ('uo2heu', 'u_uo2-heu'),
+    'c_V': ('v',),
+    'c_W': ('w',),
+    'c_Y_in_Y3Al5O12': ('y(y3al5o12', 'yyag'),
     'c_Y_in_YH2': ('yyh2', 'y-yh2', 'y(yh2)'),
+    'c_Zn': ('zn',),
+    'c_Zn_in_ZnS': ('zn(zns-sph', 'znzns'),
+    'c_Zr': ('zr',),
     'c_Zr_in_ZrC': ('zrzrc', 'zr(zrc)'),
     'c_Zr_in_ZrH': ('zrzrh', 'zr-zrh', 'zr-h', 'zr/h', 'zr(zrh)'),
-    'c_Zr_in_ZrH2': ('zrzrh2', 'zr(zrh2)'),
-    'c_Zr_in_ZrHx': ('zrzrhx', 'zr(zrhx)'),
+    'c_Zr_in_ZrH2': ('zrzrh2', 'zr(zrh2)', 'zr(zrh2) i'),
+    'c_Zr_in_ZrHx': ('zrzrhx', 'zr(zrhx)', 'zr(zrh15)'),
 }
 
 

openmc/deplete/abc.py

@@ -600,7 +600,7 @@ class Integrator(ABC):
         User-supplied functions are expected to have the following signature:
         ``solver(A, n0, t) -> n1`` where
 
-            * ``A`` is a :class:`scipy.sparse.csc_matrix` making up the
+            * ``A`` is a :class:`scipy.sparse.csc_array` making up the
               depletion matrix
             * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
               for a given material in atoms/cm3
@@ -1134,7 +1134,7 @@ class SIIntegrator(Integrator):
         User-supplied functions are expected to have the following signature:
         ``solver(A, n0, t) -> n1`` where
 
-            * ``A`` is a :class:`scipy.sparse.csc_matrix` making up the
+            * ``A`` is a :class:`scipy.sparse.csc_array` making up the
               depletion matrix
             * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
               for a given material in atoms/cm3
@@ -1297,7 +1297,7 @@ class DepSystemSolver(ABC):
 
         Parameters
         ----------
-        A : scipy.sparse.csc_matrix
+        A : scipy.sparse.csc_array
             Sparse transmutation matrix ``A[j, i]`` describing rates at
             which isotope ``i`` transmutes to isotope ``j``
         n0 : numpy.ndarray

openmc/deplete/chain.py

@@ -17,13 +17,13 @@ from warnings import warn
 from typing import List
 
 import lxml.etree as ET
-import scipy.sparse as sp
 
 from openmc.checkvalue import check_type, check_greater_than, PathLike
 from openmc.data import gnds_name, zam
 from openmc.exceptions import DataError
 from .nuclide import FissionYieldDistribution, Nuclide
 from .._xml import get_text
+from .._sparse_compat import csc_array, dok_array
 import openmc.data
 
 
@@ -619,7 +619,7 @@ class Chain:
 
         Returns
         -------
-        scipy.sparse.csc_matrix
+        scipy.sparse.csc_array
             Sparse matrix representing depletion.
 
         See Also
@@ -713,7 +713,7 @@ class Chain:
                 reactions.clear()
 
         # Return CSC representation instead of DOK
-        return sp.csc_matrix((vals, (rows, cols)), shape=(n, n))
+        return csc_array((vals, (rows, cols)), shape=(n, n))
 
     def add_redox_term(self, matrix, buffer, oxidation_states):
         r"""Adds a redox term to the depletion matrix from data contained in
@@ -731,7 +731,7 @@ class Chain:
 
         Parameters
         ----------
-        matrix : scipy.sparse.csc_matrix
+        matrix : scipy.sparse.csc_array
             Sparse matrix representing depletion
         buffer : dict
             Dictionary of buffer nuclides used to maintain anoins net balance.
@@ -743,7 +743,7 @@ class Chain:
             states as integers (e.g., +1, 0).
         Returns
         -------
-        matrix : scipy.sparse.csc_matrix
+        matrix : scipy.sparse.csc_array
             Sparse matrix with redox term added
         """
         # Elements list with the same size as self.nuclides
@@ -769,7 +769,7 @@ class Chain:
         for nuc, idx in buffer_idx.items():
             array[idx] -= redox_change * buffer[nuc] / os[idx]
 
-        return sp.csc_matrix(array)
+        return csc_array(array)
 
     def form_rr_term(self, tr_rates, current_timestep, mats):
         """Function to form the transfer rate term matrices.
@@ -800,13 +800,13 @@ class Chain:
 
         Returns
         -------
-        scipy.sparse.csc_matrix
+        scipy.sparse.csc_array
             Sparse matrix representing transfer term.
 
         """
         # Use DOK as intermediate representation
         n = len(self)
-        matrix = sp.dok_matrix((n, n))
+        matrix = dok_array((n, n))
 
         for i, nuc in enumerate(self.nuclides):
             elm = re.split(r'\d+', nuc.name)[0]
@@ -857,7 +857,7 @@ class Chain:
 
         Returns
         -------
-        scipy.sparse.csc_matrix
+        scipy.sparse.csc_array
             Sparse vector representing external source term.
 
         """
@@ -865,7 +865,7 @@ class Chain:
             return
         # Use DOK as intermediate representation
         n = len(self)
-        vector = sp.dok_matrix((n, 1))
+        vector = dok_array((n, 1))
 
         for i, nuc in enumerate(self.nuclides):
             # Build source term vector

openmc/deplete/cram.py

@@ -6,11 +6,11 @@ Implements two different forms of CRAM for use in openmc.deplete.
 import numbers
 
 import numpy as np
-import scipy.sparse as sp
 import scipy.sparse.linalg as sla
 
 from openmc.checkvalue import check_type, check_length
 from .abc import DepSystemSolver
+from .._sparse_compat import csc_array, eye_array
 
 __all__ = ["CRAM16", "CRAM48", "Cram16Solver", "Cram48Solver", "IPFCramSolver"]
 
@@ -60,7 +60,7 @@ class IPFCramSolver(DepSystemSolver):
 
         Parameters
         ----------
-        A : scipy.sparse.csr_matrix
+        A : scipy.sparse.csc_array
             Sparse transmutation matrix ``A[j, i]`` desribing rates at
             which isotope ``i`` transmutes to isotope ``j``
         n0 : numpy.ndarray
@@ -75,9 +75,9 @@ class IPFCramSolver(DepSystemSolver):
             Final compositions after ``dt``
 
         """
-        A = dt * sp.csc_matrix(A, dtype=np.float64)
+        A = dt * csc_array(A, dtype=np.float64)
         y = n0.copy()
-        ident = sp.eye(A.shape[0], format='csc')
+        ident = eye_array(A.shape[0], format='csc')
         for alpha, theta in zip(self.alpha, self.theta):
             y += 2*np.real(alpha*sla.spsolve(A - theta*ident, y))
         return y * self.alpha0

openmc/deplete/microxs.py

@@ -177,8 +177,9 @@ def get_microxs_and_flux(
         if not openmc.lib.is_initialized:
             run_kwargs.setdefault('cwd', temp_dir)
 
-        # Run transport simulation
+        # Run transport simulation and synchronize
         statepoint_path = model.run(**run_kwargs)
+        comm.barrier()
 
         if comm.rank == 0:
             # Move the statepoint file if it is being saved to a specific path

openmc/deplete/pool.py

@@ -5,10 +5,11 @@ Provided to avoid some circular imports
 from itertools import repeat, starmap
 from multiprocessing import Pool
 
-from scipy.sparse import bmat, hstack, vstack, csc_matrix
 import numpy as np
+from scipy.sparse import hstack
 
 from openmc.mpi import comm
+from .._sparse_compat import block_array
 
 # Configurable switch that enables / disables the use of
 # multiprocessing routines during depletion
@@ -159,7 +160,7 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
                             cols.append(None)
 
                     rows.append(cols)
-                matrix = bmat(rows)
+                matrix = block_array(rows)
 
                 # Concatenate vectors of nuclides in one
                 n_multi = np.concatenate(n)
@@ -194,7 +195,7 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
         # of the nuclide vectors
         for i, matrix in enumerate(matrices):
             if not np.equal(*matrix.shape):
-                matrices[i] = vstack([matrix, csc_matrix([0]*matrix.shape[1])])
+                matrix.resize(matrix.shape[1], matrix.shape[1])
                 n[i] = np.append(n[i], 1.0)
 
     inputs = zip(matrices, n, repeat(dt))

openmc/deplete/r2s.py

@@ -392,6 +392,7 @@ class R2SManager:
         )
         output_path = output_dir / 'depletion_results.h5'
         integrator.integrate(final_step=False, path=output_path)
+        comm.barrier()
 
         # Get depletion results
         self.results['depletion_results'] = Results(output_path)

openmc/examples.py

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

openmc/executor.py

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

openmc/filter.py

@@ -833,6 +833,25 @@ class MeshFilter(Filter):
     translation : Iterable of float
         This array specifies a vector that is used to translate (shift) the mesh
         for this filter
+    rotation : Iterable of float
+        This array specifies the angles in degrees about the x, y, and z axes
+        that the mesh should be rotated. The rotation applied is an intrinsic
+        rotation with specified Tait-Bryan angles. That is to say, if the angles
+        are :math:`(\phi, \theta, \psi)`, then the rotation matrix applied is
+        :math:`R_z(\psi) R_y(\theta) R_x(\phi)` or
+
+        .. math::
+
+           \left [ \begin{array}{ccc} \cos\theta \cos\psi & -\cos\phi \sin\psi
+           + \sin\phi \sin\theta \cos\psi & \sin\phi \sin\psi + \cos\phi
+           \sin\theta \cos\psi \\ \cos\theta \sin\psi & \cos\phi \cos\psi +
+           \sin\phi \sin\theta \sin\psi & -\sin\phi \cos\psi + \cos\phi
+           \sin\theta \sin\psi \\ -\sin\theta & \sin\phi \cos\theta & \cos\phi
+           \cos\theta \end{array} \right ]
+
+        A rotation matrix can also be specified directly by setting this
+        attribute to a nested list (or 2D numpy array) that specifies each
+        element of the matrix.
     bins : list of tuple
         A list of mesh indices for each filter bin, e.g. [(1, 1, 1), (2, 1, 1),
         ...]
@@ -845,6 +864,7 @@ class MeshFilter(Filter):
         self.mesh = mesh
         self.id = filter_id
         self._translation = None
+        self._rotation = None
 
     def __hash__(self):
         string = type(self).__name__ + '\n'
@@ -856,6 +876,7 @@ class MeshFilter(Filter):
         string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id)
         string += '{: <16}=\t{}\n'.format('\tID', self.id)
         string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation)
+        string += '{: <16}=\t{}\n'.format('\tRotation', self.rotation)
         return string
 
     @classmethod
@@ -879,6 +900,10 @@ class MeshFilter(Filter):
         if translation:
             out.translation = translation[()]
 
+        rotation = group.get('rotation')
+        if rotation:
+            out.rotation = rotation[()]
+
         return out
 
     @property
@@ -911,6 +936,15 @@ class MeshFilter(Filter):
         cv.check_length('mesh filter translation', t, 3)
         self._translation = np.asarray(t)
 
+    @property
+    def rotation(self):
+        return self._rotation
+
+    @rotation.setter
+    def rotation(self, rotation):
+        cv.check_length('mesh filter rotation', rotation, 3)
+        self._rotation = np.asarray(rotation)
+
     def can_merge(self, other):
         # Mesh filters cannot have more than one bin
         return False
@@ -996,6 +1030,8 @@ class MeshFilter(Filter):
         subelement.text = str(self.mesh.id)
         if self.translation is not None:
             element.set('translation', ' '.join(map(str, self.translation)))
+        if self.rotation is not None:
+            element.set('rotation', ' '.join(map(str, self.rotation.ravel())))
         return element
 
     @classmethod
@@ -1008,6 +1044,13 @@ class MeshFilter(Filter):
         translation = get_elem_list(elem, "translation", float) or []
         if translation:
             out.translation = translation
+
+        rotation = get_elem_list(elem, 'rotation', float) or []
+        if rotation:
+            if len(rotation) == 3:
+                out.rotation = rotation
+            elif len(rotation) == 9:
+                out.rotation = np.array(rotation).reshape(3, 3)
         return out
 
 
@@ -1839,12 +1882,21 @@ class DistribcellFilter(Filter):
 
     @property
     def paths(self):
-        return self._paths
+        if self._paths is None:
+            if not hasattr(self, '_geometry'):
+                raise ValueError(
+                    "Model must be exported before the 'paths' attribute is" \
+                    "available for a DistribcellFilter.")
 
-    @paths.setter
-    def paths(self, paths):
-        cv.check_iterable_type('paths', paths, str)
-        self._paths = paths
+            # Determine paths for cell instances
+            self._geometry.determine_paths()
+
+            # Get paths for the corresponding cell
+            cell_id = self.bins[0]
+            cell = self._geometry.get_all_cells()[cell_id]
+            self._paths = cell.paths
+
+        return self._paths
 
     @Filter.bins.setter
     def bins(self, bins):

openmc/lib/core.py

@@ -7,6 +7,7 @@ import os
 from pathlib import Path
 from random import getrandbits
 from tempfile import TemporaryDirectory
+import traceback as tb
 
 import numpy as np
 from numpy.ctypeslib import as_array
@@ -700,6 +701,15 @@ class TemporarySession:
         if self.already_initialized:
             return
 
+        # If an exception occurred, abort all ranks immediately
+        if exc_type is not None:
+            # Print exception info on the rank that failed
+            tb.print_exception(exc_type, exc_value, traceback)
+            sys.stdout.flush()
+
+            # Abort all MPI processes
+            self.comm.Abort(1)
+
         try:
             finalize()
         finally:

openmc/lib/filter.py

@@ -97,6 +97,14 @@ _dll.openmc_mesh_filter_get_translation.errcheck = _error_handler
 _dll.openmc_mesh_filter_set_translation.argtypes = [c_int32, POINTER(c_double*3)]
 _dll.openmc_mesh_filter_set_translation.restype = c_int
 _dll.openmc_mesh_filter_set_translation.errcheck = _error_handler
+_dll.openmc_mesh_filter_get_rotation.argtypes = [c_int32, POINTER(c_double),
+    POINTER(c_size_t)]
+_dll.openmc_mesh_filter_get_rotation.restype = c_int
+_dll.openmc_mesh_filter_get_rotation.errcheck = _error_handler
+_dll.openmc_mesh_filter_set_rotation.argtypes = [
+    c_int32, POINTER(c_double), c_size_t]
+_dll.openmc_mesh_filter_set_rotation.restype = c_int
+_dll.openmc_mesh_filter_set_rotation.errcheck = _error_handler
 _dll.openmc_meshborn_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
 _dll.openmc_meshborn_filter_get_mesh.restype = c_int
 _dll.openmc_meshborn_filter_get_mesh.errcheck = _error_handler
@@ -393,6 +401,10 @@ class MeshFilter(Filter):
         Mesh used for the filter
     translation : Iterable of float
         3-D coordinates of the translation vector
+    rotation : Iterable of float
+        The rotation matrix or angles of the filter mesh. This can either be 
+        a fully specified 3 x 3 rotation matrix or an Iterable of length 3 
+        with the angles in degrees about the x, y, and z axes, respectively.
 
     """
     filter_type = 'mesh'
@@ -422,6 +434,34 @@ class MeshFilter(Filter):
     def translation(self, translation):
         _dll.openmc_mesh_filter_set_translation(self._index, (c_double*3)(*translation))
 
+    @property
+    def rotation(self):
+        rotation_data = np.zeros(12)
+        rot_size = c_size_t()
+
+        _dll.openmc_mesh_filter_get_rotation(
+            self._index, rotation_data.ctypes.data_as(POINTER(c_double)),
+            rot_size)
+        rot_size = rot_size.value
+
+        if rot_size == 9:
+            return rotation_data[:rot_size].shape(3, 3)
+        elif rot_size in (0, 12):
+            # If size is 0, rotation_data[9:] will be zeros. This indicates no
+            # rotation and is the most straightforward way to always return
+            # an iterable of floats
+            return rotation_data[9:]
+        else:
+            raise ValueError(
+                f'Invalid size of rotation matrix: {rot_size}')
+        
+    @rotation.setter
+    def rotation(self, rotation_data):
+        flat_rotation = np.asarray(rotation_data, dtype=float).flatten()
+
+        _dll.openmc_mesh_filter_set_rotation(
+            self._index, flat_rotation.ctypes.data_as(POINTER(c_double)),
+            c_size_t(len(flat_rotation)))
 
 class MeshBornFilter(Filter):
     """MeshBorn filter stored internally.

openmc/lib/plot.py

@@ -52,7 +52,7 @@ class _PlotBase(Structure):
 
     C-Type Attributes
     -----------------
-    origin : openmc.lib.plot._Position
+    origin_ : openmc.lib.plot._Position
         A position defining the origin of the plot.
     width_ : openmc.lib.plot._Position
         The width of the plot along the x, y, and z axes, respectively
@@ -60,6 +60,8 @@ class _PlotBase(Structure):
         The axes basis of the plot view.
     pixels_ : c_size_t[3]
         The resolution of the plot in the horizontal and vertical dimensions
+    color_overlaps_ : c_bool
+        Whether to assign unique IDs (-3) to overlapping regions.
     level_ : c_int
         The universe level for the plot view
 
@@ -187,14 +189,6 @@ class _PlotBase(Structure):
     def color_overlaps(self, color_overlaps):
         self.color_overlaps_ = color_overlaps
 
-    @property
-    def color_overlaps(self):
-        return self.color_overlaps_
-
-    @color_overlaps.setter
-    def color_overlaps(self, val):
-        self.color_overlaps_ = val
-
     def __repr__(self):
         out_str = ["-----",
                    "Plot:",

openmc/material.py

@@ -60,6 +60,26 @@ class Material(IDManagerMixin):
     temperature : float, optional
         Temperature of the material in Kelvin. If not specified, the material
         inherits the default temperature applied to the model.
+    density : float, optional
+        Density of the material (units defined separately)
+    density_units : str
+        Units used for `density`. Can be one of 'g/cm3', 'g/cc', 'kg/m3',
+        'atom/b-cm', 'atom/cm3', 'sum', or 'macro'.  The 'macro' unit only
+        applies in the case of a multi-group calculation. Defaults to 'sum'.
+    depletable : bool, optional
+        Indicate whether the material is depletable. Defaults to False.
+    volume : float, optional
+        Volume of the material in cm^3. This can either be set manually or
+        calculated in a stochastic volume calculation and added via the
+        :meth:`Material.add_volume_information` method.
+    components : dict of str to float or dict
+        Dictionary mapping element or nuclide names to their atom or weight
+        percent. To specify enrichment of an element, the entry of
+        ``components`` for that element must instead be a dictionary containing
+        the keyword arguments as well as a value for ``'percent'``
+    percent_type : {'ao', 'wo'}
+        Whether the values in `components` should be interpreted as atom percent
+        ('ao') or weight percent ('wo').
 
     Attributes
     ----------
@@ -111,17 +131,28 @@ class Material(IDManagerMixin):
     next_id = 1
     used_ids = set()
 
-    def __init__(self, material_id=None, name='', temperature=None):
+    def __init__(
+        self,
+        material_id: int | None = None,
+        name: str = "",
+        temperature: float | None = None,
+        density: float | None = None,
+        density_units: str = "sum",
+        depletable: bool | None = False,
+        volume: float | None = None,
+        components: dict | None = None,
+        percent_type: str = "ao",
+    ):
         # Initialize class attributes
         self.id = material_id
         self.name = name
         self.temperature = temperature
         self._density = None
-        self._density_units = 'sum'
-        self._depletable = False
+        self._density_units = density_units
+        self._depletable = depletable
         self._paths = None
         self._num_instances = None
-        self._volume = None
+        self._volume = volume
         self._atoms = {}
         self._isotropic = []
         self._ncrystal_cfg = None
@@ -136,6 +167,15 @@ class Material(IDManagerMixin):
         # If specified, a list of table names
         self._sab = []
 
+        # Set density if provided
+        if density is not None:
+            self.set_density(density_units, density)
+
+        # Add components if provided
+        if components is not None:
+            self.add_components(components, percent_type=percent_type)
+
+
     def __repr__(self) -> str:
         string = 'Material\n'
         string += '{: <16}=\t{}\n'.format('\tID', self._id)

openmc/model/funcs.py

@@ -39,8 +39,8 @@ def borated_water(boron_ppm, temperature=293., pressure=0.1013, temp_unit='K',
     press_unit : {'MPa', 'psi'}
         The units used for the `pressure` argument.
     density : float
-        Water density in [g / cm^3].  If specified, this value overrides the
-        temperature and pressure arguments.
+        Water density in [g / cm^3].  If specified, this value overrides 
+        the value that is computed from the temperature and pressure arguments.
     **kwargs
         All keyword arguments are passed to the created Material object.
 
@@ -95,10 +95,7 @@ def borated_water(boron_ppm, temperature=293., pressure=0.1013, temp_unit='K',
     frac_B = boron_ppm * 1e-6 / M_B
 
     # Build the material.
-    if density is None:
-        out = openmc.Material(temperature=T, **kwargs)
-    else:
-        out = openmc.Material(**kwargs)
+    out = openmc.Material(temperature=T, **kwargs)
     out.add_element('H', frac_H, 'ao')
     out.add_element('O', frac_O, 'ao')
     out.add_element('B', frac_B, 'ao')

openmc/model/model.py

@@ -23,7 +23,7 @@ from openmc.dummy_comm import DummyCommunicator
 from openmc.executor import _process_CLI_arguments
 from openmc.checkvalue import check_type, check_value, PathLike
 from openmc.exceptions import InvalidIDError
-from openmc.plots import add_plot_params, _BASIS_INDICES
+from openmc.plots import add_plot_params, _BASIS_INDICES, id_map_to_rgb
 from openmc.utility_funcs import change_directory
 
 
@@ -546,6 +546,13 @@ class Model:
                 depletion_operator.cleanup_when_done = True
                 depletion_operator.finalize()
 
+    def _link_geometry_to_filters(self):
+        """Establishes a link between distribcell filters and the geometry"""
+        for tally in self.tallies:
+            for f in tally.filters:
+                if isinstance(f, openmc.DistribcellFilter):
+                    f._geometry = self.geometry
+
     def export_to_xml(self, directory: PathLike = '.', remove_surfs: bool = False,
                       nuclides_to_ignore: Iterable[str] | None = None):
         """Export model to separate XML files.
@@ -587,6 +594,8 @@ class Model:
         if self.plots:
             self.plots.export_to_xml(d)
 
+        self._link_geometry_to_filters()
+
     def export_to_model_xml(self, path: PathLike = 'model.xml', remove_surfs: bool = False,
                             nuclides_to_ignore: Iterable[str] | None = None):
         """Export model to a single XML file.
@@ -666,6 +675,8 @@ class Model:
                 fh.write(ET.tostring(plots_element, encoding="unicode"))
             fh.write("</model>\n")
 
+        self._link_geometry_to_filters()
+
     def import_properties(self, filename: PathLike):
         """Import physical properties
 
@@ -1026,6 +1037,7 @@ class Model:
         width: Sequence[float] | None = None,
         pixels: int | Sequence[int] = 40000,
         basis: str = 'xy',
+        color_overlaps: bool = False,
         **init_kwargs
     ) -> np.ndarray:
         """Generate an ID map for domains based on the plot parameters
@@ -1054,6 +1066,10 @@ class Model:
             total and the image aspect ratio based on the width argument.
         basis : {'xy', 'yz', 'xz'}, optional
             Basis of the plot.
+        color_overlaps : bool, optional
+            Whether to assign unique IDs (-3) to overlapping regions. If False,
+            overlapping regions will be assigned the ID of the lowest-numbered
+            cell that occupies that region. Defaults to False.
         **init_kwargs
             Keyword arguments passed to :meth:`Model.init_lib`.
 
@@ -1078,6 +1094,7 @@ class Model:
         plot_obj.h_res = pixels[0]
         plot_obj.v_res = pixels[1]
         plot_obj.basis = basis
+        plot_obj.color_overlaps = color_overlaps
 
         # Silence output by default. Also set arguments to start in volume
         # calculation mode to avoid loading cross sections
@@ -1097,13 +1114,12 @@ class Model:
         color_by: str = 'cell',
         colors: dict | None = None,
         seed: int | None = None,
-        openmc_exec: PathLike = 'openmc',
         axes=None,
         legend: bool = False,
         axis_units: str = 'cm',
         outline: bool | str = False,
         show_overlaps: bool = False,
-        overlap_color: Sequence[int] | str | None = None,
+        overlap_color: Sequence[int] | str = (255, 0, 0),
         n_samples: int | None = None,
         plane_tolerance: float = 1.,
         legend_kwargs: dict | None = None,
@@ -1115,7 +1131,6 @@ class Model:
 
         .. versionadded:: 0.15.1
         """
-        import matplotlib.image as mpimg
         import matplotlib.patches as mpatches
         import matplotlib.pyplot as plt
 
@@ -1145,125 +1160,108 @@ class Model:
         y_min = (origin[y] - 0.5*width[1]) * axis_scaling_factor[axis_units]
         y_max = (origin[y] + 0.5*width[1]) * axis_scaling_factor[axis_units]
 
-        # Determine whether any materials contains macroscopic data and if so,
-        # set energy mode accordingly
-        _energy_mode = self.settings._energy_mode
-        for mat in self.geometry.get_all_materials().values():
-            if mat._macroscopic is not None:
-                self.settings.energy_mode = 'multi-group'
-                break
+        # Get ID map from the C API
+        id_map = self.id_map(
+            origin=origin,
+            width=width,
+            pixels=pixels,
+            basis=basis,
+            color_overlaps=show_overlaps
+        )
 
-        with TemporaryDirectory() as tmpdir:
-            _plot_seed = self.settings.plot_seed
-            if seed is not None:
-                self.settings.plot_seed = seed
-
-            # Create plot object matching passed arguments
-            plot = openmc.Plot()
-            plot.origin = origin
-            plot.width = width
-            plot.pixels = pixels
-            plot.basis = basis
+        # Generate colors if not provided
+        if colors is None and seed is not None:
+            # Use the colorize method to generate random colors
+            plot = openmc.SlicePlot()
             plot.color_by = color_by
-            plot.show_overlaps = show_overlaps
-            if overlap_color is not None:
-                plot.overlap_color = overlap_color
-            if colors is not None:
-                plot.colors = colors
-            self.plots.append(plot)
+            plot.colorize(self.geometry, seed=seed)
+            colors = plot.colors
 
-            # Run OpenMC in geometry plotting mode
-            self.plot_geometry(False, cwd=tmpdir, openmc_exec=openmc_exec)
+        # Convert ID map to RGB image
+        img = id_map_to_rgb(
+            id_map=id_map, 
+            color_by=color_by, 
+            colors=colors,
+            overlap_color=overlap_color
+        )
 
-            # Undo changes to model
-            self.plots.pop()
-            self.settings._plot_seed = _plot_seed
-            self.settings._energy_mode = _energy_mode
+        # Create a figure sized such that the size of the axes within
+        # exactly matches the number of pixels specified
+        if axes is None:
+            px = 1/plt.rcParams['figure.dpi']
+            fig, axes = plt.subplots()
+            axes.set_xlabel(xlabel)
+            axes.set_ylabel(ylabel)
+            params = fig.subplotpars
+            width_px = pixels[0]*px/(params.right - params.left)
+            height_px = pixels[1]*px/(params.top - params.bottom)
+            fig.set_size_inches(width_px, height_px)
 
-            # Read image from file
-            img_path = Path(tmpdir) / f'plot_{plot.id}.png'
-            if not img_path.is_file():
-                img_path = img_path.with_suffix('.ppm')
-            img = mpimg.imread(str(img_path))
+        if outline:
+            # Combine R, G, B values into a single int for contour detection
+            rgb = (img * 256).astype(int)
+            image_value = (rgb[..., 0] << 16) + \
+                (rgb[..., 1] << 8) + (rgb[..., 2])
 
-            # Create a figure sized such that the size of the axes within
-            # exactly matches the number of pixels specified
-            if axes is None:
-                px = 1/plt.rcParams['figure.dpi']
-                fig, axes = plt.subplots()
-                axes.set_xlabel(xlabel)
-                axes.set_ylabel(ylabel)
-                params = fig.subplotpars
-                width = pixels[0]*px/(params.right - params.left)
-                height = pixels[1]*px/(params.top - params.bottom)
-                fig.set_size_inches(width, height)
+            # Set default arguments for contour()
+            if contour_kwargs is None:
+                contour_kwargs = {}
+            contour_kwargs.setdefault('colors', 'k')
+            contour_kwargs.setdefault('linestyles', 'solid')
+            contour_kwargs.setdefault('algorithm', 'serial')
 
-            if outline:
-                # Combine R, G, B values into a single int
-                rgb = (img * 256).astype(int)
-                image_value = (rgb[..., 0] << 16) + \
-                    (rgb[..., 1] << 8) + (rgb[..., 2])
+            axes.contour(
+                image_value,
+                origin="upper",
+                levels=np.unique(image_value),
+                extent=(x_min, x_max, y_min, y_max),
+                **contour_kwargs
+            )
+            
+            # If only showing outline, set the axis limits and aspect explicitly
+            if outline == 'only':
+                axes.set_xlim(x_min, x_max)
+                axes.set_ylim(y_min, y_max)
+                axes.set_aspect('equal')
 
-                # Set default arguments for contour()
-                if contour_kwargs is None:
-                    contour_kwargs = {}
-                contour_kwargs.setdefault('colors', 'k')
-                contour_kwargs.setdefault('linestyles', 'solid')
-                contour_kwargs.setdefault('algorithm', 'serial')
+        # Add legend showing which colors represent which material or cell
+        if legend:
+            if colors is None or len(colors) == 0:
+                raise ValueError("Must pass 'colors' dictionary if you "
+                                 "are adding a legend via legend=True.")
 
-                axes.contour(
-                    image_value,
-                    origin="upper",
-                    levels=np.unique(image_value),
-                    extent=(x_min, x_max, y_min, y_max),
-                    **contour_kwargs
-                )
+            if color_by == "cell":
+                expected_key_type = openmc.Cell
+            else:
+                expected_key_type = openmc.Material
 
-            # add legend showing which colors represent which material
-            # or cell if that was requested
-            if legend:
-                if plot.colors == {}:
-                    raise ValueError("Must pass 'colors' dictionary if you "
-                                     "are adding a legend via legend=True.")
+            patches = []
+            for key, color in colors.items():
+                if isinstance(key, int):
+                    raise TypeError(
+                        "Cannot use IDs in colors dict for auto legend.")
+                elif not isinstance(key, expected_key_type):
+                    raise TypeError(
+                        "Color dict key type does not match color_by")
 
-                if color_by == "cell":
-                    expected_key_type = openmc.Cell
+                # this works whether we're doing cells or materials
+                label = key.name if key.name != '' else key.id
+
+                # matplotlib takes RGB on 0-1 scale rather than 0-255
+                if len(color) == 3 and not isinstance(color, str):
+                    scaled_color = (
+                        color[0]/255, color[1]/255, color[2]/255)
                 else:
-                    expected_key_type = openmc.Material
+                    scaled_color = color
 
-                patches = []
-                for key, color in plot.colors.items():
+                key_patch = mpatches.Patch(color=scaled_color, label=label)
+                patches.append(key_patch)
 
-                    if isinstance(key, int):
-                        raise TypeError(
-                            "Cannot use IDs in colors dict for auto legend.")
-                    elif not isinstance(key, expected_key_type):
-                        raise TypeError(
-                            "Color dict key type does not match color_by")
-
-                    # this works whether we're doing cells or materials
-                    label = key.name if key.name != '' else key.id
-
-                    # matplotlib takes RGB on 0-1 scale rather than 0-255. at
-                    # this point PlotBase has already checked that 3-tuple
-                    # based colors are already valid, so if the length is three
-                    # then we know it just needs to be converted to the 0-1
-                    # format.
-                    if len(color) == 3 and not isinstance(color, str):
-                        scaled_color = (
-                            color[0]/255, color[1]/255, color[2]/255)
-                    else:
-                        scaled_color = color
-
-                    key_patch = mpatches.Patch(color=scaled_color, label=label)
-                    patches.append(key_patch)
-
-                axes.legend(handles=patches, **legend_kwargs)
-
-            # Plot image and return the axes
-            if outline != 'only':
-                axes.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs)
+            axes.legend(handles=patches, **legend_kwargs)
 
+        # Plot image and return the axes
+        if outline != 'only':
+            axes.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs)
 
         if n_samples:
             # Sample external source particles
@@ -1276,8 +1274,8 @@ class Model:
             tol = plane_tolerance
             for particle in particles:
                 if (slice_value - tol < particle.r[z] < slice_value + tol):
-                    xs.append(particle.r[x])
-                    ys.append(particle.r[y])
+                    xs.append(particle.r[x] * axis_scaling_factor[axis_units])
+                    ys.append(particle.r[y] * axis_scaling_factor[axis_units])
             axes.scatter(xs, ys, **source_kwargs)
 
         return axes
@@ -1687,6 +1685,91 @@ class Model:
                 self.geometry.get_all_materials().values()
             )
 
+    def _create_mgxs_sources(
+        self,
+        groups: openmc.mgxs.EnergyGroups,
+        spatial_dist: openmc.stats.Spatial,
+        source_energy: openmc.stats.Univariate | None = None,
+    ) -> list[openmc.IndependentSource]:
+        """Create a list of independent sources to use with MGXS generation.
+
+        Note that in all cases, a discrete source that is uniform over all
+        energy groups is created (strength = 0.01) to ensure that total cross
+        sections are generated for all energy groups. In the case that the user
+        has provided a source_energy distribution as an argument, an additional
+        source (strength = 0.99) is created using that energy distribution. If
+        the user has not provided a source_energy distribution, but the model
+        has sources defined, and all of those sources are of IndependentSource
+        type, then additional sources are created based on the model's existing
+        sources, keeping their energy distributions but replacing their
+        spatial/angular distributions, with their combined strength being 0.99.
+        If the user has not provided a source_energy distribution and no sources
+        are defined on the model and the run mode is 'eigenvalue', then a
+        default Watt spectrum source (strength = 0.99) is added.
+
+        Parameters
+        ----------
+        groups : openmc.mgxs.EnergyGroups
+            Energy group structure for the MGXS.
+        spatial_dist : openmc.stats.Spatial
+            Spatial distribution to use for all sources.
+        source_energy : openmc.stats.Univariate, optional
+            Energy distribution to use when generating MGXS data, replacing any
+            existing sources in the model.
+
+        Returns
+        -------
+        list[openmc.IndependentSource]
+            A list of independent sources to use for MGXS generation.
+        """
+        # Make a discrete source that is uniform over the bins of the group structure
+        midpoints = []
+        strengths = []
+        for i in range(groups.num_groups):
+            bounds = groups.get_group_bounds(i+1)
+            midpoints.append((bounds[0] + bounds[1]) / 2.0)
+            strengths.append(1.0)
+
+        uniform_energy = openmc.stats.Discrete(x=midpoints, p=strengths)
+        uniform_distribution = openmc.IndependentSource(spatial_dist, energy=uniform_energy, strength=0.01)
+        sources = [uniform_distribution]
+
+        # If the user provided an energy distribution, use that
+        if source_energy is not None:
+            user_energy = openmc.IndependentSource(
+                space=spatial_dist, energy=source_energy, strength=0.99)
+            sources.append(user_energy)
+
+        # If the user did not provide an energy distribution, create sources
+        # based on what is in their model, keeping the energy spectrum but
+        # replacing the spatial/angular distributions. We only do this if ALL
+        # sources are of IndependentSource type, as we can't pull the energy
+        # distribution from e.g. CompiledSource or FileSource types.
+        else:
+            if self.settings.source is not None:
+                for src in self.settings.source:
+                    if not isinstance(src, openmc.IndependentSource):
+                        break
+                else:
+                    n_user_sources = len(self.settings.source)
+                    for src in self.settings.source:
+                        # Create a new IndependentSource with adjusted strength, space, and angle
+                        user_source = openmc.IndependentSource(
+                            space=spatial_dist,
+                            energy=src.energy,
+                            strength=0.99 / n_user_sources
+                        )
+                        sources.append(user_source)
+            else:
+                # No user sources defined. If we are in eigenvalue mode, then use the default Watt spectrum.
+                if self.settings.run_mode == 'eigenvalue':
+                    watt_energy = openmc.stats.Watt()
+                    watt_source = openmc.IndependentSource(
+                        space=spatial_dist, energy=watt_energy, strength=0.99)
+                    sources.append(watt_source)
+
+        return sources
+
     def _generate_infinite_medium_mgxs(
         self,
         groups: openmc.mgxs.EnergyGroups,
@@ -1694,6 +1777,7 @@ class Model:
         mgxs_path: PathLike,
         correction: str | None,
         directory: PathLike,
+        source_energy: openmc.stats.Univariate | None = None,
     ):
         """Generate a MGXS library by running multiple OpenMC simulations, each
         representing an infinite medium simulation of a single isolated
@@ -1702,6 +1786,20 @@ class Model:
         method that ignores all spatial self shielding effects and all resonance
         shielding effects between materials.
 
+        Note that in all cases, a discrete source that is uniform over all
+        energy groups is created (strength = 0.01) to ensure that total cross
+        sections are generated for all energy groups. In the case that the user
+        has provided a source_energy distribution as an argument, an additional
+        source (strength = 0.99) is created using that energy distribution. If
+        the user has not provided a source_energy distribution, but the model
+        has sources defined, and all of those sources are of IndependentSource
+        type, then additional sources are created based on the model's existing
+        sources, keeping their energy distributions but replacing their
+        spatial/angular distributions, with their combined strength being 0.99.
+        If the user has not provided a source_energy distribution and no sources
+        are defined on the model and the run mode is 'eigenvalue', then a
+        default Watt spectrum source (strength = 0.99) is added.
+
         Parameters
         ----------
         groups : openmc.mgxs.EnergyGroups
@@ -1715,9 +1813,10 @@ class Model:
             "P0".
         directory : str
             Directory to run the simulation in, so as to contain XML files.
+        source_energy : openmc.stats.Univariate, optional
+            Energy distribution to use when generating MGXS data, replacing any
+            existing sources in the model.
         """
-        warnings.warn("The infinite medium method of generating MGXS may hang "
-                      "if a material has a k-infinity > 1.0.")
         mgxs_sets = []
         for material in self.materials:
             model = openmc.Model()
@@ -1728,20 +1827,16 @@ class Model:
             # Settings
             model.settings.batches = 100
             model.settings.particles = nparticles
+
+            model.settings.source = self._create_mgxs_sources(
+                groups,
+                spatial_dist=openmc.stats.Point(),
+                source_energy=source_energy
+            )
+
             model.settings.run_mode = 'fixed source'
+            model.settings.create_fission_neutrons = False
 
-            # Make a discrete source that is uniform over the bins of the group structure
-            n_groups = groups.num_groups
-            midpoints = []
-            strengths = []
-            for i in range(n_groups):
-                bounds = groups.get_group_bounds(i+1)
-                midpoints.append((bounds[0] + bounds[1]) / 2.0)
-                strengths.append(1.0)
-
-            energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths)
-            model.settings.source = openmc.IndependentSource(
-                space=openmc.stats.Point(), energy=energy_distribution)
             model.settings.output = {'summary': True, 'tallies': False}
 
             # Geometry
@@ -1891,6 +1986,7 @@ class Model:
         mgxs_path: PathLike,
         correction: str | None,
         directory: PathLike,
+        source_energy: openmc.stats.Univariate | None = None,
     ) -> None:
         """Generate MGXS assuming a stochastic "sandwich" of materials in a layered
         slab geometry. While geometry-specific spatial shielding effects are not
@@ -1915,6 +2011,23 @@ class Model:
             "P0".
         directory : str
             Directory to run the simulation in, so as to contain XML files.
+        source_energy : openmc.stats.Univariate, optional
+            Energy distribution to use when generating MGXS data, replacing any
+            existing sources in the model. In all cases, a discrete source that
+            is uniform over all energy groups is created (strength = 0.01) to
+            ensure that total cross sections are generated for all energy
+            groups. In the case that the user has provided a source_energy
+            distribution as an argument, an additional source (strength = 0.99)
+            is created using that energy distribution. If the user has not
+            provided a source_energy distribution, but the model has sources
+            defined, and all of those sources are of IndependentSource type,
+            then additional sources are created based on the model's existing
+            sources, keeping their energy distributions but replacing their
+            spatial/angular distributions, with their combined strength being
+            0.99. If the user has not provided a source_energy distribution and
+            no sources are defined on the model and the run mode is
+            'eigenvalue', then a default Watt spectrum source (strength = 0.99)
+            is added.
         """
         model = openmc.Model()
         model.materials = self.materials
@@ -1924,24 +2037,20 @@ class Model:
         model.settings.inactive = 100
         model.settings.particles = nparticles
         model.settings.output = {'summary': True, 'tallies': False}
-        model.settings.run_mode = self.settings.run_mode
 
         # Stochastic slab geometry
         model.geometry, spatial_distribution = Model._create_stochastic_slab_geometry(
             model.materials)
 
-        # Make a discrete source that is uniform over the bins of the group structure
-        n_groups = groups.num_groups
-        midpoints = []
-        strengths = []
-        for i in range(n_groups):
-            bounds = groups.get_group_bounds(i+1)
-            midpoints.append((bounds[0] + bounds[1]) / 2.0)
-            strengths.append(1.0)
+        # Define the sources
+        model.settings.source = self._create_mgxs_sources(
+            groups,
+            spatial_dist=spatial_distribution,
+            source_energy=source_energy
+        )
 
-        energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths)
-        model.settings.source = [openmc.IndependentSource(
-            space=spatial_distribution, energy=energy_distribution, strength=1.0)]
+        model.settings.run_mode = 'fixed source'
+        model.settings.create_fission_neutrons = False
 
         model.settings.output = {'summary': True, 'tallies': False}
 
@@ -2099,6 +2208,7 @@ class Model:
         overwrite_mgxs_library: bool = False,
         mgxs_path: PathLike = "mgxs.h5",
         correction: str | None = None,
+        source_energy: openmc.stats.Univariate | None = None,
     ):
         """Convert all materials from continuous energy to multigroup.
 
@@ -2112,11 +2222,33 @@ class Model:
         groups : openmc.mgxs.EnergyGroups or str, optional
             Energy group structure for the MGXS or the name of the group
             structure (based on keys from openmc.mgxs.GROUP_STRUCTURES).
+        nparticles : int, optional
+            Number of particles to simulate per batch when generating MGXS.
+        overwrite_mgxs_library : bool, optional
+            Whether to overwrite an existing MGXS library file.
         mgxs_path : str, optional
-            Filename of the mgxs.h5 library file.
+            Path to the mgxs.h5 library file.
         correction : str, optional
             Transport correction to apply to the MGXS. Options are None and
             "P0".
+        source_energy : openmc.stats.Univariate, optional
+            Energy distribution to use when generating MGXS data, replacing any
+            existing sources in the model. In all cases, a discrete source that
+            is uniform over all energy groups is created (strength = 0.01) to
+            ensure that total cross sections are generated for all energy
+            groups. In the case that the user has provided a source_energy
+            distribution as an argument, an additional source (strength = 0.99)
+            is created using that energy distribution. If the user has not
+            provided a source_energy distribution, but the model has sources
+            defined, and all of those sources are of IndependentSource type,
+            then additional sources are created based on the model's existing
+            sources, keeping their energy distributions but replacing their
+            spatial/angular distributions, with their combined strength being
+            0.99. If the user has not provided a source_energy distribution and
+            no sources are defined on the model and the run mode is
+            'eigenvalue', then a default Watt spectrum source (strength = 0.99)
+            is added. Note that this argument is only used when using the
+            "stochastic_slab" or "infinite_medium" MGXS generation methods.
         """
         if isinstance(groups, str):
             groups = openmc.mgxs.EnergyGroups(groups)
@@ -2146,13 +2278,13 @@ class Model:
             if not Path(mgxs_path).is_file() or overwrite_mgxs_library:
                 if method == "infinite_medium":
                     self._generate_infinite_medium_mgxs(
-                        groups, nparticles, mgxs_path, correction, tmpdir)
+                        groups, nparticles, mgxs_path, correction, tmpdir, source_energy)
                 elif method == "material_wise":
                     self._generate_material_wise_mgxs(
                         groups, nparticles, mgxs_path, correction, tmpdir)
                 elif method == "stochastic_slab":
                     self._generate_stochastic_slab_mgxs(
-                        groups, nparticles, mgxs_path, correction, tmpdir)
+                        groups, nparticles, mgxs_path, correction, tmpdir, source_energy)
                 else:
                     raise ValueError(
                         f'MGXS generation method "{method}" not recognized')

openmc/plots.py

@@ -1,7 +1,8 @@
-from collections.abc import Iterable, Mapping
+from collections.abc import Iterable, Mapping, Sequence
 from numbers import Integral, Real
 from pathlib import Path
 from textwrap import dedent
+import warnings
 
 import h5py
 import lxml.etree as ET
@@ -354,6 +355,86 @@ def voxel_to_vtk(voxel_file: PathLike, output: PathLike = 'plot.vti'):
     return output
 
 
+def id_map_to_rgb(
+    id_map: np.ndarray,
+    color_by: str = 'cell',
+    colors: dict | None = None,
+    overlap_color: Sequence[int] | str = (255, 0, 0)
+) -> np.ndarray:
+    """Convert ID map array to RGB image array.
+
+    Parameters
+    ----------
+    id_map : numpy.ndarray
+        Array with shape (v_pixels, h_pixels, 3) containing cell IDs,
+        cell instances, and material IDs
+    color_by : {'cell', 'material'}
+        Whether to color by cell or material
+    colors : dict, optional
+        Dictionary mapping cells/materials to colors
+    overlap_color : sequence of int or str, optional
+        Color to use for overlaps. Defaults to red (255, 0, 0).
+
+    Returns
+    -------
+    numpy.ndarray
+        RGB image array with shape (v_pixels, h_pixels, 3) with values
+        in range [0, 1] for matplotlib
+    """
+    # Initialize RGB array with white background (values between 0 and 1 for matplotlib)
+    img = np.ones(id_map.shape, dtype=float)
+    
+    # Get the appropriate index based on color_by
+    if color_by == 'cell':
+        id_index = 0  # Cell IDs are in the first channel
+    elif color_by == 'material':
+        id_index = 2  # Material IDs are in the third channel
+    else:
+        raise ValueError("color_by must be either 'cell' or 'material'")
+    
+    # Get all unique IDs in the plot
+    unique_ids = np.unique(id_map[:, :, id_index])
+    
+    # Generate default colors if not provided
+    if colors is None:
+        colors = {}
+    
+    # Convert colors dict to use IDs as keys
+    color_map = {}
+    for key, color in colors.items():
+        if isinstance(key, (openmc.Cell, openmc.Material)):
+            color_map[key.id] = color
+        else:
+            color_map[key] = color
+    
+    # Generate random colors for IDs not in color_map
+    rng = np.random.RandomState(1)
+    for uid in unique_ids:
+        if uid > 0 and uid not in color_map:
+            color_map[uid] = rng.randint(0, 256, (3,))
+    
+    # Apply colors to each pixel
+    for uid in unique_ids:
+        if uid == -1:  # Background/void
+            continue
+        elif uid == -3:  # Overlap (only present if color_overlaps was True)
+            if isinstance(overlap_color, str):
+                rgb = _SVG_COLORS[overlap_color.lower()]
+            else:
+                rgb = overlap_color
+            mask = id_map[:, :, id_index] == uid
+            img[mask] = np.array(rgb) / 255.0
+        elif uid in color_map:
+            color = color_map[uid]
+            if isinstance(color, str):
+                rgb = _SVG_COLORS[color.lower()]
+            else:
+                rgb = color
+            mask = id_map[:, :, id_index] == uid
+            img[mask] = np.array(rgb) / 255.0
+    
+    return img
+
 class PlotBase(IDManagerMixin):
     """
     Parameters
@@ -626,14 +707,15 @@ class PlotBase(IDManagerMixin):
         return element
 
 
-class Plot(PlotBase):
-    """Definition of a finite region of space to be plotted.
+class SlicePlot(PlotBase):
+    """Definition of a 2D slice plot of the geometry.
 
-    OpenMC is capable of generating two-dimensional slice plots, or
-    three-dimensional voxel or projection plots. Colors that are used in plots can be given as
-    RGB tuples, e.g. (255, 255, 255) would be white, or by a string indicating a
+    Colors that are used in plots can be given as RGB tuples, e.g.
+    (255, 255, 255) would be white, or by a string indicating a
     valid `SVG color <https://www.w3.org/TR/SVG11/types.html#ColorKeywords>`_.
 
+    .. versionadded:: 0.15.4
+
     Parameters
     ----------
     plot_id : int
@@ -648,7 +730,7 @@ class Plot(PlotBase):
     name : str
         Name of the plot
     pixels : Iterable of int
-        Number of pixels to use in each direction
+        Number of pixels to use in each direction (2 values)
     filename : str
         Path to write the plot to
     color_by : {'cell', 'material'}
@@ -671,11 +753,9 @@ class Plot(PlotBase):
     level : int
         Universe depth to plot at
     width : Iterable of float
-        Width of the plot in each basis direction
+        Width of the plot in each basis direction (2 values)
     origin : tuple or list of ndarray
-        Origin (center) of the plot
-    type : {'slice', 'voxel'}
-        The type of the plot
+        Origin (center) of the plot (3 values)
     basis : {'xy', 'xz', 'yz'}
         The basis directions for the plot
     meshlines : dict
@@ -688,10 +768,37 @@ class Plot(PlotBase):
         super().__init__(plot_id, name)
         self._width = [4.0, 4.0]
         self._origin = [0., 0., 0.]
-        self._type = 'slice'
         self._basis = 'xy'
         self._meshlines = None
 
+    @property
+    def type(self):
+        warnings.warn(
+            "The 'type' attribute is deprecated and will be removed in a future version. "
+            "This is a SlicePlot instance.",
+            FutureWarning, stacklevel=2
+        )
+        return 'slice'
+
+    @type.setter
+    def type(self, value):
+        raise TypeError(
+            "Setting plot.type is no longer supported. "
+            "Use openmc.SlicePlot() for 2D slice plots or openmc.VoxelPlot() for 3D voxel plots."
+        )
+
+    @property
+    def pixels(self):
+        return self._pixels
+
+    @pixels.setter
+    def pixels(self, pixels):
+        cv.check_type('plot pixels', pixels, Iterable, Integral)
+        cv.check_length('plot pixels', pixels, 2, 2)
+        for dim in pixels:
+            cv.check_greater_than('plot pixels', dim, 0)
+        self._pixels = pixels
+
     @property
     def width(self):
         return self._width
@@ -699,7 +806,7 @@ class Plot(PlotBase):
     @width.setter
     def width(self, width):
         cv.check_type('plot width', width, Iterable, Real)
-        cv.check_length('plot width', width, 2, 3)
+        cv.check_length('plot width', width, 2, 2)
         self._width = width
 
     @property
@@ -712,15 +819,6 @@ class Plot(PlotBase):
         cv.check_length('plot origin', origin, 3)
         self._origin = origin
 
-    @property
-    def type(self):
-        return self._type
-
-    @type.setter
-    def type(self, plottype):
-        cv.check_value('plot type', plottype, ['slice', 'voxel'])
-        self._type = plottype
-
     @property
     def basis(self):
         return self._basis
@@ -763,11 +861,10 @@ class Plot(PlotBase):
         self._meshlines = meshlines
 
     def __repr__(self):
-        string = 'Plot\n'
+        string = 'SlicePlot\n'
         string += '{: <16}=\t{}\n'.format('\tID', self._id)
         string += '{: <16}=\t{}\n'.format('\tName', self._name)
         string += '{: <16}=\t{}\n'.format('\tFilename', self._filename)
-        string += '{: <16}=\t{}\n'.format('\tType', self._type)
         string += '{: <16}=\t{}\n'.format('\tBasis', self._basis)
         string += '{: <16}=\t{}\n'.format('\tWidth', self._width)
         string += '{: <16}=\t{}\n'.format('\tOrigin', self._origin)
@@ -883,7 +980,7 @@ class Plot(PlotBase):
                 self._colors[domain] = (r, g, b)
 
     def to_xml_element(self):
-        """Return XML representation of the slice/voxel plot
+        """Return XML representation of the slice plot
 
         Returns
         -------
@@ -893,10 +990,8 @@ class Plot(PlotBase):
         """
 
         element = super().to_xml_element()
-        element.set("type", self._type)
-
-        if self._type == 'slice':
-            element.set("basis", self._basis)
+        element.set("type", "slice")
+        element.set("basis", self._basis)
 
         subelement = ET.SubElement(element, "origin")
         subelement.text = ' '.join(map(str, self._origin))
@@ -942,8 +1037,8 @@ class Plot(PlotBase):
 
         Returns
         -------
-        openmc.Plot
-            Plot object
+        openmc.SlicePlot
+            SlicePlot object
 
         """
         plot_id = int(get_text(elem, "id"))
@@ -952,9 +1047,7 @@ class Plot(PlotBase):
         if "filename" in elem.keys():
             plot.filename = get_text(elem, "filename")
         plot.color_by = get_text(elem, "color_by")
-        plot.type = get_text(elem, "type")
-        if plot.type == 'slice':
-            plot.basis = get_text(elem, "basis")
+        plot.basis = get_text(elem, "basis")
 
         plot.origin = tuple(get_elem_list(elem, "origin", float))
         plot.width = tuple(get_elem_list(elem, "width", float))
@@ -1036,9 +1129,215 @@ class Plot(PlotBase):
         # Return produced image
         return _get_plot_image(self, cwd)
 
+
+
+class VoxelPlot(PlotBase):
+    """Definition of a 3D voxel plot of the geometry.
+
+    Colors that are used in plots can be given as RGB tuples, e.g.
+    (255, 255, 255) would be white, or by a string indicating a
+    valid `SVG color <https://www.w3.org/TR/SVG11/types.html#ColorKeywords>`_.
+
+    .. versionadded:: 0.15.1
+
+    Parameters
+    ----------
+    plot_id : int
+        Unique identifier for the plot
+    name : str
+        Name of the plot
+
+    Attributes
+    ----------
+    id : int
+        Unique identifier
+    name : str
+        Name of the plot
+    pixels : Iterable of int
+        Number of pixels to use in each direction (3 values)
+    filename : str
+        Path to write the plot to
+    color_by : {'cell', 'material'}
+        Indicate whether the plot should be colored by cell or by material
+    background : Iterable of int or str
+        Color of the background
+    mask_components : Iterable of openmc.Cell or openmc.Material or int
+        The cells or materials (or corresponding IDs) to mask
+    mask_background : Iterable of int or str
+        Color to apply to all cells/materials listed in mask_components
+    show_overlaps : bool
+        Indicate whether or not overlapping regions are shown
+    overlap_color : Iterable of int or str
+        Color to apply to overlapping regions
+    colors : dict
+        Dictionary indicating that certain cells/materials should be
+        displayed with a particular color. The keys can be of type
+        :class:`~openmc.Cell`, :class:`~openmc.Material`, or int (ID for a
+        cell/material).
+    level : int
+        Universe depth to plot at
+    width : Iterable of float
+        Width of the plot in each dimension (3 values)
+    origin : tuple or list of ndarray
+        Origin (center) of the plot (3 values)
+
+    """
+
+    def __init__(self, plot_id=None, name=''):
+        super().__init__(plot_id, name)
+        self._width = [4.0, 4.0, 4.0]
+        self._origin = [0., 0., 0.]
+        self._pixels = [400, 400, 400]
+
+    @property
+    def pixels(self):
+        return self._pixels
+
+    @pixels.setter
+    def pixels(self, pixels):
+        cv.check_type('plot pixels', pixels, Iterable, Integral)
+        cv.check_length('plot pixels', pixels, 3, 3)
+        for dim in pixels:
+            cv.check_greater_than('plot pixels', dim, 0)
+        self._pixels = pixels
+
+    @property
+    def width(self):
+        return self._width
+
+    @width.setter
+    def width(self, width):
+        cv.check_type('plot width', width, Iterable, Real)
+        cv.check_length('plot width', width, 3, 3)
+        self._width = width
+
+    @property
+    def origin(self):
+        return self._origin
+
+    @origin.setter
+    def origin(self, origin):
+        cv.check_type('plot origin', origin, Iterable, Real)
+        cv.check_length('plot origin', origin, 3)
+        self._origin = origin
+
+    def __repr__(self):
+        string = 'VoxelPlot\n'
+        string += '{: <16}=\t{}\n'.format('\tID', self._id)
+        string += '{: <16}=\t{}\n'.format('\tName', self._name)
+        string += '{: <16}=\t{}\n'.format('\tFilename', self._filename)
+        string += '{: <16}=\t{}\n'.format('\tWidth', self._width)
+        string += '{: <16}=\t{}\n'.format('\tOrigin', self._origin)
+        string += '{: <16}=\t{}\n'.format('\tPixels', self._pixels)
+        string += '{: <16}=\t{}\n'.format('\tColor by', self._color_by)
+        string += '{: <16}=\t{}\n'.format('\tBackground', self._background)
+        string += '{: <16}=\t{}\n'.format('\tMask components',
+                                          self._mask_components)
+        string += '{: <16}=\t{}\n'.format('\tMask background',
+                                          self._mask_background)
+        string += '{: <16}=\t{}\n'.format('\tOverlap Color',
+                                          self._overlap_color)
+        string += '{: <16}=\t{}\n'.format('\tColors', self._colors)
+        string += '{: <16}=\t{}\n'.format('\tLevel', self._level)
+        return string
+
+    def to_xml_element(self):
+        """Return XML representation of the voxel plot
+
+        Returns
+        -------
+        element : lxml.etree._Element
+            XML element containing plot data
+
+        """
+
+        element = super().to_xml_element()
+        element.set("type", "voxel")
+
+        subelement = ET.SubElement(element, "origin")
+        subelement.text = ' '.join(map(str, self._origin))
+
+        subelement = ET.SubElement(element, "width")
+        subelement.text = ' '.join(map(str, self._width))
+
+        if self._colors:
+            self._colors_to_xml(element)
+
+        if self._show_overlaps:
+            subelement = ET.SubElement(element, "show_overlaps")
+            subelement.text = "true"
+
+            if self._overlap_color is not None:
+                color = self._overlap_color
+                if isinstance(color, str):
+                    color = _SVG_COLORS[color.lower()]
+                subelement = ET.SubElement(element, "overlap_color")
+                subelement.text = ' '.join(str(x) for x in color)
+
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate plot object from an XML element
+
+        Parameters
+        ----------
+        elem : lxml.etree._Element
+            XML element
+
+        Returns
+        -------
+        openmc.VoxelPlot
+            VoxelPlot object
+
+        """
+        plot_id = int(get_text(elem, "id"))
+        name = get_text(elem, 'name', '')
+        plot = cls(plot_id, name)
+        if "filename" in elem.keys():
+            plot.filename = get_text(elem, "filename")
+        plot.color_by = get_text(elem, "color_by")
+
+        plot.origin = tuple(get_elem_list(elem, "origin", float))
+        plot.width = tuple(get_elem_list(elem, "width", float))
+        plot.pixels = tuple(get_elem_list(elem, "pixels"))
+        background = get_elem_list(elem, "background")
+        if background is not None:
+            plot._background = tuple(background)
+
+        # Set plot colors
+        colors = {}
+        for color_elem in elem.findall("color"):
+            uid = int(get_text(color_elem, "id"))
+            colors[uid] = tuple(get_elem_list(color_elem, "rgb", int))
+        plot.colors = colors
+
+        # Set masking information
+        mask_elem = elem.find("mask")
+        if mask_elem is not None:
+            plot.mask_components = get_elem_list(mask_elem, "components", int)
+            background = get_elem_list(mask_elem, "background", int)
+            if background is not None:
+                plot.mask_background = tuple(background)
+
+        # show overlaps
+        overlap = get_text(elem, "show_overlaps")
+        if overlap is not None:
+            plot.show_overlaps = (overlap in ('true', '1'))
+        overlap_color = get_elem_list(elem, "overlap_color", int)
+        if overlap_color is not None:
+            plot.overlap_color = tuple(overlap_color)
+
+        # Set universe level
+        level = get_text(elem, "level")
+        if level is not None:
+            plot.level = int(level)
+
+        return plot
+
     def to_vtk(self, output: PathLike | None = None,
                openmc_exec: str = 'openmc', cwd: str = '.'):
-        """Render plot as an voxel image
+        """Render plot as a voxel image
 
         This method runs OpenMC in plotting mode to produce a .vti file.
 
@@ -1059,10 +1358,6 @@ class Plot(PlotBase):
             Path of the .vti file produced
 
         """
-        if self.type != 'voxel':
-            raise ValueError(
-                'Generating a VTK file only works for voxel plots')
-
         # Create plots.xml
         Plots([self]).export_to_xml(cwd)
 
@@ -1082,6 +1377,20 @@ class Plot(PlotBase):
         return voxel_to_vtk(h5_voxel_file, output)
 
 
+def Plot(plot_id=None, name=''):
+    """Legacy Plot class for backward compatibility.
+
+    .. deprecated:: 0.15.4
+        Use :class:`SlicePlot` for 2D slice plots or :class:`VoxelPlot` for 3D voxel plots.
+
+    """
+    warnings.warn(
+        "The Plot class is deprecated. Use SlicePlot for 2D slice plots "
+        "or VoxelPlot for 3D voxel plots.", FutureWarning
+    )
+    return SlicePlot(plot_id, name)
+
+
 class RayTracePlot(PlotBase):
     """Definition of a camera's view of OpenMC geometry
 
@@ -1737,16 +2046,16 @@ class SolidRayTracePlot(RayTracePlot):
 
 
 class Plots(cv.CheckedList):
-    """Collection of Plots used for an OpenMC simulation.
+    """Collection of plots used for an OpenMC simulation.
 
     This class corresponds directly to the plots.xml input file. It can be
     thought of as a normal Python list where each member is inherits from
     :class:`PlotBase`. It behaves like a list as the following example
     demonstrates:
 
-    >>> xz_plot = openmc.Plot()
-    >>> big_plot = openmc.Plot()
-    >>> small_plot = openmc.Plot()
+    >>> xz_plot = openmc.SlicePlot()
+    >>> big_plot = openmc.VoxelPlot()
+    >>> small_plot = openmc.SlicePlot()
     >>> p = openmc.Plots((xz_plot, big_plot))
     >>> p.append(small_plot)
     >>> small_plot = p.pop()
@@ -1782,7 +2091,7 @@ class Plots(cv.CheckedList):
         ----------
         index : int
             Index in list
-        plot : openmc.Plot
+        plot : openmc.PlotBase
             Plot to insert
 
         """
@@ -1903,8 +2212,13 @@ class Plots(cv.CheckedList):
                 plots.append(WireframeRayTracePlot.from_xml_element(e))
             elif plot_type == 'solid_raytrace':
                 plots.append(SolidRayTracePlot.from_xml_element(e))
-            elif plot_type in ('slice', 'voxel'):
-                plots.append(Plot.from_xml_element(e))
+            elif plot_type == 'slice':
+                plots.append(SlicePlot.from_xml_element(e))
+            elif plot_type == 'voxel':
+                plots.append(VoxelPlot.from_xml_element(e))
+            elif plot_type is None:
+                # For backward compatibility, assume slice if no type specified
+                plots.append(SlicePlot.from_xml_element(e))
             else:
                 raise ValueError("Unknown plot type: {}".format(plot_type))
         return plots

openmc/settings.py

@@ -182,7 +182,7 @@ class Settings:
         Options for configuring the random ray solver. Acceptable keys are:
 
         :distance_inactive:
-            Indicates the total active distance in [cm] a ray should travel
+            Indicates the total inactive distance in [cm] a ray should travel
         :distance_active:
             Indicates the total active distance in [cm] a ray should travel
         :ray_source:
@@ -1278,9 +1278,12 @@ class Settings:
         return self._weight_windows_file
 
     @weight_windows_file.setter
-    def weight_windows_file(self, value: PathLike):
-        cv.check_type('weight windows file', value, PathLike)
-        self._weight_windows_file = input_path(value)
+    def weight_windows_file(self, value: PathLike | None):
+        if value is None:
+            self._weight_windows_file = None
+        else:
+            cv.check_type('weight windows file', value, PathLike)
+            self._weight_windows_file = input_path(value)
 
     @property
     def weight_window_generators(self) -> list[WeightWindowGenerator]:

openmc/statepoint.py

@@ -723,7 +723,7 @@ class StatePoint:
                     cell = cells[cell_id]
                     if not cell._paths:
                         summary.geometry.determine_paths()
-                    tally_filter.paths = cell.paths
+                    tally_filter._paths = cell.paths
 
         self._summary = summary
 

openmc/surface.py

@@ -123,9 +123,8 @@ class Surface(IDManagerMixin, ABC):
     boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
-        freely pass through the surface. Note that periodic boundary conditions
-        can only be applied to x-, y-, and z-planes, and only axis-aligned
-        periodicity is supported.
+        freely pass through the surface. Note that only axis-aligned
+        periodicity is supported around the x-, y-, and z-axes.
     albedo : float, optional
         Albedo of the surfaces as a ratio of particle weight after interaction
         with the surface to the initial weight. Values must be positive. Only
@@ -822,8 +821,7 @@ class XPlane(PlaneMixin, Surface):
     boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
-        freely pass through the surface. Only axis-aligned periodicity is
-        supported, i.e., x-planes can only be paired with x-planes.
+        freely pass through the surface.
     albedo : float, optional
         Albedo of the surfaces as a ratio of particle weight after interaction
         with the surface to the initial weight. Values must be positive. Only
@@ -887,8 +885,7 @@ class YPlane(PlaneMixin, Surface):
     boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
-        freely pass through the surface. Only axis-aligned periodicity is
-        supported, i.e., y-planes can only be paired with y-planes.
+        freely pass through the surface.
     albedo : float, optional
         Albedo of the surfaces as a ratio of particle weight after interaction
         with the surface to the initial weight. Values must be positive. Only
@@ -952,8 +949,7 @@ class ZPlane(PlaneMixin, Surface):
     boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
-        freely pass through the surface. Only axis-aligned periodicity is
-        supported, i.e., z-planes can only be paired with z-planes.
+        freely pass through the surface.
     albedo : float, optional
         Albedo of the surfaces as a ratio of particle weight after interaction
         with the surface to the initial weight. Values must be positive. Only

openmc/tallies.py

@@ -12,11 +12,11 @@ import lxml.etree as ET
 import h5py
 import numpy as np
 import pandas as pd
-import scipy.sparse as sps
 from scipy.stats import chi2, norm
 
 import openmc
 import openmc.checkvalue as cv
+from ._sparse_compat import lil_array
 from ._xml import clean_indentation, get_elem_list, get_text
 from .mixin import IDManagerMixin
 from .mesh import MeshBase
@@ -290,7 +290,7 @@ class Tally(IDManagerMixin):
 
     @property
     def num_nuclides(self):
-        return len(self._nuclides)
+        return max(len(self._nuclides), 1)
 
     @property
     def scores(self):
@@ -393,6 +393,11 @@ class Tally(IDManagerMixin):
             group = f[f'tallies/tally {self.id}']
             self._num_realizations = int(group['n_realizations'][()])
 
+            for filt in self.filters:
+                if isinstance(filt, openmc.DistribcellFilter):
+                    filter_group = f[f'tallies/filters/filter {filt.id}']
+                    filt._num_bins = int(filter_group['n_bins'][()])
+
             # Update nuclides
             nuclide_names = group['nuclides'][()]
             self._nuclides = [name.decode().strip() for name in nuclide_names]
@@ -430,10 +435,10 @@ class Tally(IDManagerMixin):
 
             # Convert NumPy arrays to SciPy sparse LIL matrices
             if self.sparse:
-                self._sum = sps.lil_matrix(self._sum.flatten(), self._sum.shape)
-                self._sum_sq = sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape)
-                self._sum_third = sps.lil_matrix(self._sum_third.flatten(), self._sum_third.shape)
-                self._sum_fourth = sps.lil_matrix(self.sum_fourth.flatten(), self._sum_fourth.shape)
+                self._sum = lil_array(self._sum.flatten(), self._sum.shape)
+                self._sum_sq = lil_array(self._sum_sq.flatten(), self._sum_sq.shape)
+                self._sum_third = lil_array(self._sum_third.flatten(), self._sum_third.shape)
+                self._sum_fourth = lil_array(self._sum_fourth.flatten(), self._sum_fourth.shape)
 
             # Read simulation time (needed for figure of merit)
             self._simulation_time = f["runtime"]["simulation"][()]
@@ -529,8 +534,7 @@ class Tally(IDManagerMixin):
 
             # Convert NumPy array to SciPy sparse LIL matrix
             if self.sparse:
-                self._mean = sps.lil_matrix(self._mean.flatten(),
-                                            self._mean.shape)
+                self._mean = lil_array(self._mean.flatten(), self._mean.shape)
 
         if self.sparse:
             return np.reshape(self._mean.toarray(), self.shape)
@@ -551,8 +555,7 @@ class Tally(IDManagerMixin):
 
             # Convert NumPy array to SciPy sparse LIL matrix
             if self.sparse:
-                self._std_dev = sps.lil_matrix(self._std_dev.flatten(),
-                                               self._std_dev.shape)
+                self._std_dev = lil_array(self._std_dev.flatten(), self._std_dev.shape)
 
             self.with_batch_statistics = True
 
@@ -583,7 +586,7 @@ class Tally(IDManagerMixin):
             self._vov[mask] = numerator[mask]/denominator[mask] - 1.0/n
 
             if self.sparse:
-                self._vov = sps.lil_matrix(self._vov.flatten(), self._vov.shape)
+                self._vov = lil_array(self._vov.flatten(), self._vov.shape)
 
         if self.sparse:
             return np.reshape(self._vov.toarray(), self.shape)
@@ -958,22 +961,17 @@ class Tally(IDManagerMixin):
         # Convert NumPy arrays to SciPy sparse LIL matrices
         if sparse and not self.sparse:
             if self._sum is not None:
-                self._sum = sps.lil_matrix(self._sum.flatten(), self._sum.shape)
+                self._sum = lil_array(self._sum.flatten(), self._sum.shape)
             if self._sum_sq is not None:
-                self._sum_sq = sps.lil_matrix(self._sum_sq.flatten(),
-                                              self._sum_sq.shape)
+                self._sum_sq = lil_array(self._sum_sq.flatten(), self._sum_sq.shape)
             if self._sum_third is not None:
-                self._sum_third = sps.lil_matrix(self._sum_third.flatten(),
-                                                 self._sum_third.shape)
+                self._sum_third = lil_array(self._sum_third.flatten(), self._sum_third.shape)
             if self._sum_fourth is not None:
-                self._sum_fourth = sps.lil_matrix(self._sum_fourth.flatten(),
-                                                  self._sum_fourth.shape)
+                self._sum_fourth = lil_array(self._sum_fourth.flatten(), self._sum_fourth.shape)
             if self._mean is not None:
-                self._mean = sps.lil_matrix(self._mean.flatten(),
-                                            self._mean.shape)
+                self._mean = lil_array(self._mean.flatten(), self._mean.shape)
             if self._std_dev is not None:
-                self._std_dev = sps.lil_matrix(self._std_dev.flatten(),
-                                               self._std_dev.shape)
+                self._std_dev = lil_array(self._std_dev.flatten(), self._std_dev.shape)
 
             self._sparse = True
 
@@ -3704,43 +3702,17 @@ class Tallies(cv.CheckedList):
             if possible. Defaults to False.
 
         """
-        if not isinstance(tally, Tally):
-            msg = f'Unable to add a non-Tally "{tally}" to the Tallies instance'
-            raise TypeError(msg)
-
         if merge:
-            merged = False
-
             # Look for a tally to merge with this one
             for i, tally2 in enumerate(self):
-
                 # If a mergeable tally is found
                 if tally2.can_merge(tally):
                     # Replace tally2 with the merged tally
                     merged_tally = tally2.merge(tally)
                     self[i] = merged_tally
-                    merged = True
-                    break
+                    return
 
-            # If no mergeable tally was found, simply add this tally
-            if not merged:
-                super().append(tally)
-
-        else:
-            super().append(tally)
-
-    def insert(self, index, item):
-        """Insert tally before index
-
-        Parameters
-        ----------
-        index : int
-            Index in list
-        item : openmc.Tally
-            Tally to insert
-
-        """
-        super().insert(index, item)
+        super().append(tally)
 
     def merge_tallies(self):
         """Merge any mergeable tallies together. Note that n-way merges are

src/boundary_condition.cpp

@@ -158,63 +158,44 @@ void TranslationalPeriodicBC::handle_particle(
 // RotationalPeriodicBC implementation
 //==============================================================================
 
-RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
+RotationalPeriodicBC::RotationalPeriodicBC(
+  int i_surf, int j_surf, PeriodicAxis axis)
   : PeriodicBC(i_surf, j_surf)
 {
   Surface& surf1 {*model::surfaces[i_surf_]};
   Surface& surf2 {*model::surfaces[j_surf_]};
 
-  // Check the type of the first surface
-  bool surf1_is_xyplane;
-  if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf1)) {
-    surf1_is_xyplane = true;
-  } else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf1)) {
-    surf1_is_xyplane = true;
-  } else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf1)) {
-    surf1_is_xyplane = false;
-  } else {
-    throw std::invalid_argument(fmt::format(
-      "Surface {} is an invalid type for "
-      "rotational periodic BCs. Only x-planes, y-planes, or general planes "
-      "(that are perpendicular to z) are supported for these BCs.",
-      surf1.id_));
-  }
-
-  // Check the type of the second surface
-  bool surf2_is_xyplane;
-  if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf2)) {
-    surf2_is_xyplane = true;
-  } else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf2)) {
-    surf2_is_xyplane = true;
-  } else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf2)) {
-    surf2_is_xyplane = false;
-  } else {
-    throw std::invalid_argument(fmt::format(
-      "Surface {} is an invalid type for "
-      "rotational periodic BCs. Only x-planes, y-planes, or general planes "
-      "(that are perpendicular to z) are supported for these BCs.",
-      surf2.id_));
+  // below convention for right handed coordinate system
+  switch (axis) {
+  case x:
+    zero_axis_idx_ = 0; // x component of plane must be zero
+    axis_1_idx_ = 1;    // y component independent
+    axis_2_idx_ = 2;    // z component dependent
+    break;
+  case y:
+    // for a right handed coordinate system, z should be the independent axis
+    // but this would cause the y-rotation case to be different than the other
+    // two. using a left handed coordinate system and a negative rotation the
+    // compute angle and rotation matrix behavior mimics that of the x and z
+    // cases
+    zero_axis_idx_ = 1; // y component of plane must be zero
+    axis_1_idx_ = 0;    // x component independent
+    axis_2_idx_ = 2;    // z component dependent
+    break;
+  case z:
+    zero_axis_idx_ = 2; // z component of plane must be zero
+    axis_1_idx_ = 0;    // x component independent
+    axis_2_idx_ = 1;    // y component dependent
+    break;
+  default:
+    throw std::invalid_argument(
+      fmt::format("You've specified an axis that is not x, y, or z."));
   }
 
   // Compute the surface normal vectors and make sure they are perpendicular
-  // to the z-axis
+  // to the correct axis
   Direction norm1 = surf1.normal({0, 0, 0});
   Direction norm2 = surf2.normal({0, 0, 0});
-  if (std::abs(norm1.z) > FP_PRECISION) {
-    throw std::invalid_argument(fmt::format(
-      "Rotational periodic BCs are only "
-      "supported for rotations about the z-axis, but surface {} is not "
-      "perpendicular to the z-axis.",
-      surf1.id_));
-  }
-  if (std::abs(norm2.z) > FP_PRECISION) {
-    throw std::invalid_argument(fmt::format(
-      "Rotational periodic BCs are only "
-      "supported for rotations about the z-axis, but surface {} is not "
-      "perpendicular to the z-axis.",
-      surf2.id_));
-  }
-
   // Make sure both surfaces intersect the origin
   if (std::abs(surf1.evaluate({0, 0, 0})) > FP_COINCIDENT) {
     throw std::invalid_argument(fmt::format(
@@ -231,15 +212,8 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
       surf2.id_));
   }
 
-  // Compute the BC rotation angle.  Here it is assumed that both surface
-  // normal vectors point inwards---towards the valid geometry region.
-  // Consequently, the rotation angle is not the difference between the two
-  // normals, but is instead the difference between one normal and one
-  // anti-normal.  (An incident ray on one surface must be an outgoing ray on
-  // the other surface after rotation hence the anti-normal.)
-  double theta1 = std::atan2(norm1.y, norm1.x);
-  double theta2 = std::atan2(norm2.y, norm2.x) + PI;
-  angle_ = theta2 - theta1;
+  angle_ = compute_periodic_rotation(norm1[axis_2_idx_], norm1[axis_1_idx_],
+    norm2[axis_2_idx_], norm2[axis_1_idx_]);
 
   // Warn the user if the angle does not evenly divide a circle
   double rem = std::abs(std::remainder((2 * PI / angle_), 1.0));
@@ -251,6 +225,20 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
   }
 }
 
+double RotationalPeriodicBC::compute_periodic_rotation(
+  double rise_1, double run_1, double rise_2, double run_2) const
+{
+  // Compute the BC rotation angle.  Here it is assumed that both surface
+  // normal vectors point inwards---towards the valid geometry region.
+  // Consequently, the rotation angle is not the difference between the two
+  // normals, but is instead the difference between one normal and one
+  // anti-normal.  (An incident ray on one surface must be an outgoing ray on
+  // the other surface after rotation hence the anti-normal.)
+  double theta1 = std::atan2(rise_1, run_1);
+  double theta2 = std::atan2(rise_2, run_2) + PI;
+  return theta2 - theta1;
+}
+
 void RotationalPeriodicBC::handle_particle(
   Particle& p, const Surface& surf) const
 {
@@ -278,10 +266,16 @@ void RotationalPeriodicBC::handle_particle(
   Direction u = p.u();
   double cos_theta = std::cos(theta);
   double sin_theta = std::sin(theta);
-  Position new_r = {
-    cos_theta * r.x - sin_theta * r.y, sin_theta * r.x + cos_theta * r.y, r.z};
-  Direction new_u = {
-    cos_theta * u.x - sin_theta * u.y, sin_theta * u.x + cos_theta * u.y, u.z};
+
+  Position new_r;
+  new_r[zero_axis_idx_] = r[zero_axis_idx_];
+  new_r[axis_1_idx_] = cos_theta * r[axis_1_idx_] - sin_theta * r[axis_2_idx_];
+  new_r[axis_2_idx_] = sin_theta * r[axis_1_idx_] + cos_theta * r[axis_2_idx_];
+
+  Direction new_u;
+  new_u[zero_axis_idx_] = u[zero_axis_idx_];
+  new_u[axis_1_idx_] = cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_];
+  new_u[axis_2_idx_] = sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_];
 
   // Handle the effects of the surface albedo on the particle's weight.
   BoundaryCondition::handle_albedo(p, surf);

src/cell.cpp

@@ -57,7 +57,7 @@ void Cell::set_rotation(const vector<double>& rot)
     fatal_error(fmt::format("Non-3D rotation vector applied to cell {}", id_));
   }
 
-  // Compute and store the rotation matrix.
+  // Compute and store the inverse rotation matrix for the angles given.
   rotation_.clear();
   rotation_.reserve(rot.size() == 9 ? 9 : 12);
   if (rot.size() == 3) {
@@ -1334,14 +1334,14 @@ extern "C" int openmc_cell_bounding_box(
   bbox = c->bounding_box();
 
   // set lower left corner values
-  llc[0] = bbox.xmin;
-  llc[1] = bbox.ymin;
-  llc[2] = bbox.zmin;
+  llc[0] = bbox.min.x;
+  llc[1] = bbox.min.y;
+  llc[2] = bbox.min.z;
 
   // set upper right corner values
-  urc[0] = bbox.xmax;
-  urc[1] = bbox.ymax;
-  urc[2] = bbox.zmax;
+  urc[0] = bbox.max.x;
+  urc[1] = bbox.max.y;
+  urc[2] = bbox.max.z;
 
   return 0;
 }

src/collision_track.cpp

@@ -78,10 +78,10 @@ void write_collision_track_bank(hid_t group_id,
   hid_t banktype = h5_collision_track_banktype();
 #ifdef OPENMC_MPI
   write_bank_dataset("collision_track_bank", group_id, collision_track_bank,
-    bank_index, banktype, mpi::collision_track_site);
+    bank_index, banktype, banktype, mpi::collision_track_site);
 #else
   write_bank_dataset("collision_track_bank", group_id, collision_track_bank,
-    bank_index, banktype);
+    bank_index, banktype, banktype);
 #endif
 
   H5Tclose(banktype);

src/dagmc.cpp

@@ -753,7 +753,7 @@ BoundingBox DAGCell::bounding_box() const
   double min[3], max[3];
   rval = dagmc_ptr_->getobb(vol, min, max);
   MB_CHK_ERR_CONT(rval);
-  return {min[0], max[0], min[1], max[1], min[2], max[2]};
+  return {{min[0], min[1], min[2]}, {max[0], max[1], max[2]}};
 }
 
 //==============================================================================

src/distribution_angle.cpp

@@ -7,6 +7,7 @@
 
 #include "openmc/endf.h"
 #include "openmc/hdf5_interface.h"
+#include "openmc/math_functions.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 #include "openmc/vector.h" // for vector
@@ -64,23 +65,10 @@ AngleDistribution::AngleDistribution(hid_t group)
 
 double AngleDistribution::sample(double E, uint64_t* seed) const
 {
-  // Determine number of incoming energies
-  auto n = energy_.size();
-
-  // Find energy bin and calculate interpolation factor -- if the energy is
-  // outside the range of the tabulated energies, choose the first or last bins
+  // Find energy bin and calculate interpolation factor
   int i;
   double r;
-  if (E < energy_[0]) {
-    i = 0;
-    r = 0.0;
-  } else if (E > energy_[n - 1]) {
-    i = n - 2;
-    r = 1.0;
-  } else {
-    i = lower_bound_index(energy_.begin(), energy_.end(), E);
-    r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]);
-  }
+  get_energy_index(energy_, E, i, r);
 
   // Sample between the ith and (i+1)th bin
   if (r > prn(seed))

src/eigenvalue.cpp

@@ -403,6 +403,16 @@ void calculate_average_keff()
         t_value *
         std::sqrt(
           (simulation::k_sum[1] / n - std::pow(simulation::keff, 2)) / (n - 1));
+
+      // In some cases (such as an infinite medium problem), random ray
+      // may estimate k exactly and in an unvarying manner between iterations.
+      // In this case, the floating point roundoff between the division and the
+      // power operations may cause an extremely small negative value to occur
+      // inside the sqrt operation, leading to NaN. If this occurs, we check for
+      // it and set the std dev to zero.
+      if (!std::isfinite(simulation::keff_std)) {
+        simulation::keff_std = 0.0;
+      }
     }
   }
 }

src/finalize.cpp

@@ -142,7 +142,7 @@ int openmc_finalize()
   settings::uniform_source_sampling = false;
   settings::ufs_on = false;
   settings::urr_ptables_on = true;
-  settings::verbosity = 7;
+  settings::verbosity = -1;
   settings::weight_cutoff = 0.25;
   settings::weight_survive = 1.0;
   settings::weight_windows_file.clear();

src/geometry.cpp

@@ -480,14 +480,14 @@ extern "C" int openmc_global_bounding_box(double* llc, double* urc)
   auto bbox = model::universes.at(model::root_universe)->bounding_box();
 
   // set lower left corner values
-  llc[0] = bbox.xmin;
-  llc[1] = bbox.ymin;
-  llc[2] = bbox.zmin;
+  llc[0] = bbox.min.x;
+  llc[1] = bbox.min.y;
+  llc[2] = bbox.min.z;
 
   // set upper right corner values
-  urc[0] = bbox.xmax;
-  urc[1] = bbox.ymax;
-  urc[2] = bbox.zmax;
+  urc[0] = bbox.max.x;
+  urc[1] = bbox.max.y;
+  urc[2] = bbox.max.z;
 
   return 0;
 }

src/initialize.cpp

@@ -226,6 +226,15 @@ int parse_command_line(int argc, char* argv[])
         i += 1;
         settings::n_particles = std::stoll(argv[i]);
 
+      } else if (arg == "-q" || arg == "--verbosity") {
+        i += 1;
+        settings::verbosity = std::stoi(argv[i]);
+        if (settings::verbosity > 10 || settings::verbosity < 1) {
+          auto msg = fmt::format("Invalid verbosity: {}.", settings::verbosity);
+          strcpy(openmc_err_msg, msg.c_str());
+          return OPENMC_E_INVALID_ARGUMENT;
+        }
+
       } else if (arg == "-e" || arg == "--event") {
         settings::event_based = true;
       } else if (arg == "-r" || arg == "--restart") {
@@ -376,8 +385,10 @@ bool read_model_xml()
   auto settings_root = root.child("settings");
 
   // Verbosity
-  if (check_for_node(settings_root, "verbosity")) {
+  if (check_for_node(settings_root, "verbosity") && settings::verbosity == -1) {
     settings::verbosity = std::stoi(get_node_value(settings_root, "verbosity"));
+  } else if (settings::verbosity == -1) {
+    settings::verbosity = 7;
   }
 
   // To this point, we haven't displayed any output since we didn't know what

src/math_functions.cpp

@@ -919,4 +919,19 @@ std::complex<double> w_derivative(std::complex<double> z, int order)
   }
 }
 
+// Helper function to get index and interpolation function on an incident energy
+// grid
+void get_energy_index(
+  const vector<double>& energies, double E, int& i, double& f)
+{
+  // Get index and interpolation factor for linear-linear energy grid
+  i = 0;
+  f = 0.0;
+  if (E >= energies.front()) {
+    i = lower_bound_index(energies.begin(), energies.end(), E);
+    if (i + 1 < energies.size())
+      f = (E - energies[i]) / (energies[i + 1] - energies[i]);
+  }
+}
+
 } // namespace openmc

src/mcpl_interface.cpp

@@ -63,7 +63,7 @@ using mcpl_read_fpt = const mcpl_particle_repr_t* (*)(mcpl_file_t* file_handle);
 using mcpl_close_file_fpt = void (*)(mcpl_file_t* file_handle);
 
 using mcpl_hdr_add_data_fpt = void (*)(mcpl_outfile_t* file_handle,
-  const char* key, int32_t ldata, const char* data);
+  const char* key, uint32_t datalength, const char* data);
 using mcpl_create_outfile_fpt = mcpl_outfile_t* (*)(const char* filename);
 using mcpl_hdr_set_srcname_fpt = void (*)(
   mcpl_outfile_t* outfile_handle, const char* srcname);
@@ -150,13 +150,20 @@ struct McplApi {
       load_symbol_platform("mcpl_create_outfile"));
     hdr_set_srcname = reinterpret_cast<mcpl_hdr_set_srcname_fpt>(
       load_symbol_platform("mcpl_hdr_set_srcname"));
-    hdr_add_data = reinterpret_cast<mcpl_hdr_add_data_fpt>(
-      load_symbol_platform("mcpl_hdr_add_data"));
     add_particle = reinterpret_cast<mcpl_add_particle_fpt>(
       load_symbol_platform("mcpl_add_particle"));
     close_outfile = reinterpret_cast<mcpl_close_outfile_fpt>(
       load_symbol_platform("mcpl_close_outfile"));
 
+    // Try to load mcpl_hdr_add_data (available in MCPL >= 2.1.0)
+    // Set to nullptr if not available for graceful fallback
+    try {
+      hdr_add_data = reinterpret_cast<mcpl_hdr_add_data_fpt>(
+        load_symbol_platform("mcpl_hdr_add_data"));
+    } catch (const std::runtime_error&) {
+      hdr_add_data = nullptr;
+    }
+
     // Try to load mcpl_hdr_add_stat_sum (available in MCPL >= 2.1.0)
     // Set to nullptr if not available for graceful fallback
     try {

src/mesh.cpp

@@ -51,6 +51,7 @@
 #include "libmesh/mesh_modification.h"
 #include "libmesh/mesh_tools.h"
 #include "libmesh/numeric_vector.h"
+#include "libmesh/replicated_mesh.h"
 #endif
 
 #ifdef OPENMC_DAGMC_ENABLED
@@ -358,9 +359,10 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
   std::array<int, 3> n_rays = {nx, ny, nz};
 
   // Determine effective width of rays
-  Position width((nx > 0) ? (bbox.xmax - bbox.xmin) / nx : 0.0,
-    (ny > 0) ? (bbox.ymax - bbox.ymin) / ny : 0.0,
-    (nz > 0) ? (bbox.zmax - bbox.zmin) / nz : 0.0);
+  Position width = bbox.max - bbox.min;
+  width.x = (nx > 0) ? width.x / nx : 0.0;
+  width.y = (ny > 0) ? width.y / ny : 0.0;
+  width.z = (nz > 0) ? width.z / nz : 0.0;
 
   // Set flag for mesh being contained within model
   bool out_of_model = false;
@@ -379,15 +381,15 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
     for (int axis = 0; axis < 3; ++axis) {
       // Set starting position and direction
       site.r = {0.0, 0.0, 0.0};
-      site.r[axis] = bbox.min()[axis];
+      site.r[axis] = bbox.min[axis];
       site.u = {0.0, 0.0, 0.0};
       site.u[axis] = 1.0;
 
       // Determine width of rays and number of rays in other directions
       int ax1 = (axis + 1) % 3;
       int ax2 = (axis + 2) % 3;
-      double min1 = bbox.min()[ax1];
-      double min2 = bbox.min()[ax2];
+      double min1 = bbox.min[ax1];
+      double min2 = bbox.min[ax2];
       double d1 = width[ax1];
       double d2 = width[ax2];
       int n1 = n_rays[ax1];
@@ -432,7 +434,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
           while (true) {
             // Ray trace from r_start to r_end
             Position r0 = p.r();
-            double max_distance = bbox.max()[axis] - r0[axis];
+            double max_distance = bbox.max[axis] - r0[axis];
 
             // Find the distance to the nearest boundary
             BoundaryInfo boundary = distance_to_boundary(p);
@@ -2414,14 +2416,14 @@ extern "C" int openmc_mesh_bounding_box(int32_t index, double* ll, double* ur)
   BoundingBox bbox = model::meshes[index]->bounding_box();
 
   // set lower left corner values
-  ll[0] = bbox.xmin;
-  ll[1] = bbox.ymin;
-  ll[2] = bbox.zmin;
+  ll[0] = bbox.min.x;
+  ll[1] = bbox.min.y;
+  ll[2] = bbox.min.z;
 
   // set upper right corner values
-  ur[0] = bbox.xmax;
-  ur[1] = bbox.ymax;
-  ur[2] = bbox.zmax;
+  ur[0] = bbox.max.x;
+  ur[1] = bbox.max.y;
+  ur[2] = bbox.max.z;
   return 0;
 }
 
@@ -3435,7 +3437,7 @@ LibMesh::LibMesh(hid_t group) : UnstructuredMesh(group)
 // create the mesh from a pointer to a libMesh Mesh
 LibMesh::LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier)
 {
-  if (!dynamic_cast<libMesh::ReplicatedMesh*>(&input_mesh)) {
+  if (!input_mesh.is_replicated()) {
     fatal_error("At present LibMesh tallies require a replicated mesh. Please "
                 "ensure 'input_mesh' is a libMesh::ReplicatedMesh.");
   }

src/output.cpp

@@ -281,6 +281,7 @@ void print_usage()
       "  -t, --track            Write tracks for all particles (up to "
       "max_tracks)\n"
       "  -e, --event            Run using event-based parallelism\n"
+      "  -q, --verbosity        Output verbosity\n"
       "  -v, --version          Show version information\n"
       "  -h, --help             Show this message\n");
   }

src/random_ray/flat_source_domain.cpp

@@ -124,7 +124,9 @@ void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
         double chi = chi_[material * negroups_ + g_out];
 
         scatter_source += sigma_s * scalar_flux;
-        fission_source += nu_sigma_f * scalar_flux * chi;
+        if (settings::create_fission_neutrons) {
+          fission_source += nu_sigma_f * scalar_flux * chi;
+        }
       }
       srh.source(g_out) =
         (scatter_source + fission_source * inverse_k_eff) / sigma_t;
@@ -369,6 +371,7 @@ void FlatSourceDomain::compute_k_eff()
   // Adds entropy value to shared entropy vector in openmc namespace.
   simulation::entropy.push_back(H);
 
+  fission_rate_ = fission_rate_new;
   k_eff_ = k_eff_new;
 }
 
@@ -519,12 +522,33 @@ void FlatSourceDomain::reset_tally_volumes()
 // simulation
 double FlatSourceDomain::compute_fixed_source_normalization_factor() const
 {
-  // If we are not in fixed source mode, then there are no external sources
-  // so no normalization is needed.
-  if (settings::run_mode != RunMode::FIXED_SOURCE || adjoint_) {
+  // Eigenvalue mode normalization
+  if (settings::run_mode == RunMode::EIGENVALUE) {
+    // Normalize fluxes by total number of fission neutrons produced. This
+    // ensures consistent scaling of the eigenvector such that its magnitude is
+    // comparable to the eigenvector produced by the Monte Carlo solver.
+    // Multiplying by the eigenvalue is unintuitive, but it is necessary.
+    // If the eigenvalue is 1.2, per starting source neutron, you will
+    // generate 1.2 neutrons. Thus if we normalize to generating only ONE
+    // neutron in total for the whole domain, then we don't actually have enough
+    // flux to generate the required 1.2 neutrons. We only know the flux
+    // required to generate 1 neutron (which would have required less than one
+    // starting neutron). Thus, you have to scale the flux up by the eigenvalue
+    // such that 1.2 neutrons are generated, so as to be consistent with the
+    // bookkeeping in MC which is all done per starting source neutron (not per
+    // neutron produced).
+    return k_eff_ / (fission_rate_ * simulation_volume_);
+  }
+
+  // If we are in adjoint mode of a fixed source problem, the external
+  // source is already normalized, such that all resulting fluxes are
+  // also normalized.
+  if (adjoint_) {
     return 1.0;
   }
 
+  // Fixed source mode normalization
+
   // Step 1 is to sum over all source regions and energy groups to get the
   // total external source strength in the simulation.
   double simulation_external_source_strength = 0.0;

src/random_ray/linear_source_domain.cpp

@@ -68,9 +68,11 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
 
         // Compute source terms for flat and linear components of the flux
         scatter_flat += sigma_s * flux_flat;
-        fission_flat += nu_sigma_f * flux_flat * chi;
         scatter_linear += sigma_s * flux_linear;
-        fission_linear += nu_sigma_f * flux_linear * chi;
+        if (settings::create_fission_neutrons) {
+          fission_flat += nu_sigma_f * flux_flat * chi;
+          fission_linear += nu_sigma_f * flux_linear * chi;
+        }
       }
 
       // Compute the flat source term

src/secondary_correlated.cpp

@@ -10,6 +10,7 @@
 
 #include "openmc/endf.h"
 #include "openmc/hdf5_interface.h"
+#include "openmc/math_functions.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 
@@ -156,21 +157,10 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group)
 void CorrelatedAngleEnergy::sample(
   double E_in, double& E_out, double& mu, uint64_t* seed) const
 {
-  // Find energy bin and calculate interpolation factor -- if the energy is
-  // outside the range of the tabulated energies, choose the first or last bins
-  auto n_energy_in = energy_.size();
+  // Find energy bin and calculate interpolation factor
   int i;
   double r;
-  if (E_in < energy_[0]) {
-    i = 0;
-    r = 0.0;
-  } else if (E_in > energy_[n_energy_in - 1]) {
-    i = n_energy_in - 2;
-    r = 1.0;
-  } else {
-    i = lower_bound_index(energy_.begin(), energy_.end(), E_in);
-    r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]);
-  }
+  get_energy_index(energy_, E_in, i, r);
 
   // Sample between the ith and [i+1]th bin
   int l = r > prn(seed) ? i + 1 : i;

src/secondary_kalbach.cpp

@@ -9,6 +9,7 @@
 #include "xtensor/xview.hpp"
 
 #include "openmc/hdf5_interface.h"
+#include "openmc/math_functions.h"
 #include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
@@ -117,21 +118,10 @@ KalbachMann::KalbachMann(hid_t group)
 void KalbachMann::sample(
   double E_in, double& E_out, double& mu, uint64_t* seed) const
 {
-  // Find energy bin and calculate interpolation factor -- if the energy is
-  // outside the range of the tabulated energies, choose the first or last bins
-  auto n_energy_in = energy_.size();
+  // Find energy bin and calculate interpolation factor
   int i;
   double r;
-  if (E_in < energy_[0]) {
-    i = 0;
-    r = 0.0;
-  } else if (E_in > energy_[n_energy_in - 1]) {
-    i = n_energy_in - 2;
-    r = 1.0;
-  } else {
-    i = lower_bound_index(energy_.begin(), energy_.end(), E_in);
-    r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]);
-  }
+  get_energy_index(energy_, E_in, i, r);
 
   // Sample between the ith and [i+1]th bin
   int l = r > prn(seed) ? i + 1 : i;

src/secondary_thermal.cpp

@@ -1,6 +1,7 @@
 #include "openmc/secondary_thermal.h"
 
 #include "openmc/hdf5_interface.h"
+#include "openmc/math_functions.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 
@@ -11,20 +12,6 @@
 
 namespace openmc {
 
-// Helper function to get index on incident energy grid
-void get_energy_index(
-  const vector<double>& energies, double E, int& i, double& f)
-{
-  // Get index and interpolation factor for elastic grid
-  i = 0;
-  f = 0.0;
-  if (E >= energies.front()) {
-    i = lower_bound_index(energies.begin(), energies.end(), E);
-    if (i + 1 < energies.size())
-      f = (E - energies[i]) / (energies[i + 1] - energies[i]);
-  }
-}
-
 //==============================================================================
 // CoherentElasticAE implementation
 //==============================================================================

src/settings.cpp

@@ -145,7 +145,7 @@ int trace_gen;
 int64_t trace_particle;
 vector<array<int, 3>> track_identifiers;
 int trigger_batch_interval {1};
-int verbosity {7};
+int verbosity {-1};
 double weight_cutoff {0.25};
 double weight_survive {1.0};
 
@@ -396,8 +396,10 @@ void read_settings_xml()
   xml_node root = doc.document_element();
 
   // Verbosity
-  if (check_for_node(root, "verbosity")) {
+  if (check_for_node(root, "verbosity") && verbosity == -1) {
     verbosity = std::stoi(get_node_value(root, "verbosity"));
+  } else if (verbosity == -1) {
+    verbosity = 7;
   }
 
   // To this point, we haven't displayed any output since we didn't know what
@@ -545,6 +547,20 @@ void read_settings_xml(pugi::xml_node root)
     } else if (rel_max_lost_particles <= 0.0 || rel_max_lost_particles >= 1.0) {
       fatal_error("Relative max lost particles must be between zero and one.");
     }
+
+    // Check for user value for the number of generation of the Iterated Fission
+    // Probability (IFP) method
+    if (check_for_node(root, "ifp_n_generation")) {
+      ifp_n_generation = std::stoi(get_node_value(root, "ifp_n_generation"));
+      if (ifp_n_generation <= 0) {
+        fatal_error("'ifp_n_generation' must be greater than 0.");
+      }
+      // Avoid tallying 0 if IFP logs are not complete when active cycles start
+      if (ifp_n_generation > n_inactive) {
+        fatal_error("'ifp_n_generation' must be lower than or equal to the "
+                    "number of inactive cycles.");
+      }
+    }
   }
 
   // Copy plotting random number seed if specified
@@ -1130,20 +1146,6 @@ void read_settings_xml(pugi::xml_node root)
     temperature_range[1] = range.at(1);
   }
 
-  // Check for user value for the number of generation of the Iterated Fission
-  // Probability (IFP) method
-  if (check_for_node(root, "ifp_n_generation")) {
-    ifp_n_generation = std::stoi(get_node_value(root, "ifp_n_generation"));
-    if (ifp_n_generation <= 0) {
-      fatal_error("'ifp_n_generation' must be greater than 0.");
-    }
-    // Avoid tallying 0 if IFP logs are not complete when active cycles start
-    if (ifp_n_generation > n_inactive) {
-      fatal_error("'ifp_n_generation' must be lower than or equal to the "
-                  "number of inactive cycles.");
-    }
-  }
-
   // Check for tabular_legendre options
   if (check_for_node(root, "tabular_legendre")) {
     // Get pointer to tabular_legendre node

src/state_point.cpp

@@ -554,7 +554,7 @@ extern "C" int openmc_statepoint_load(const char* filename)
   return 0;
 }
 
-hid_t h5banktype()
+hid_t h5banktype(bool memory)
 {
   // Create compound type for position
   hid_t postype = H5Tcreate(H5T_COMPOUND, sizeof(struct Position));
@@ -569,7 +569,10 @@ hid_t h5banktype()
   // - openmc/statepoint.py
   // - docs/source/io_formats/statepoint.rst
   // - docs/source/io_formats/source.rst
-  hid_t banktype = H5Tcreate(H5T_COMPOUND, sizeof(struct SourceSite));
+  auto n = sizeof(SourceSite);
+  if (!memory)
+    n = 2 * sizeof(struct Position) + 3 * sizeof(double) + 3 * sizeof(int);
+  hid_t banktype = H5Tcreate(H5T_COMPOUND, n);
   H5Tinsert(banktype, "r", HOFFSET(SourceSite, r), postype);
   H5Tinsert(banktype, "u", HOFFSET(SourceSite, u), postype);
   H5Tinsert(banktype, "E", HOFFSET(SourceSite, E), H5T_NATIVE_DOUBLE);
@@ -641,17 +644,19 @@ void write_h5_source_point(const char* filename, span<SourceSite> source_bank,
 void write_source_bank(hid_t group_id, span<SourceSite> source_bank,
   const vector<int64_t>& bank_index)
 {
-  hid_t banktype = h5banktype();
+  hid_t membanktype = h5banktype(true);
+  hid_t filebanktype = h5banktype(false);
 
 #ifdef OPENMC_MPI
-  write_bank_dataset("source_bank", group_id, source_bank, bank_index, banktype,
-    mpi::source_site);
+  write_bank_dataset("source_bank", group_id, source_bank, bank_index,
+    membanktype, filebanktype, mpi::source_site);
 #else
-  write_bank_dataset(
-    "source_bank", group_id, source_bank, bank_index, banktype);
+  write_bank_dataset("source_bank", group_id, source_bank, bank_index,
+    membanktype, filebanktype);
 #endif
 
-  H5Tclose(banktype);
+  H5Tclose(membanktype);
+  H5Tclose(filebanktype);
 }
 
 // Determine member names of a compound HDF5 datatype
@@ -672,7 +677,7 @@ std::string dtype_member_names(hid_t dtype_id)
 void read_source_bank(
   hid_t group_id, vector<SourceSite>& sites, bool distribute)
 {
-  hid_t banktype = h5banktype();
+  hid_t banktype = h5banktype(true);
 
   // Open the dataset
   hid_t dset = H5Dopen(group_id, "source_bank", H5P_DEFAULT);

src/surface.cpp

@@ -251,9 +251,9 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
 BoundingBox SurfaceXPlane::bounding_box(bool pos_side) const
 {
   if (pos_side) {
-    return {x0_, INFTY, -INFTY, INFTY, -INFTY, INFTY};
+    return {{x0_, -INFTY, -INFTY}, {INFTY, INFTY, INFTY}};
   } else {
-    return {-INFTY, x0_, -INFTY, INFTY, -INFTY, INFTY};
+    return {{-INFTY, -INFTY, -INFTY}, {x0_, INFTY, INFTY}};
   }
 }
 
@@ -291,9 +291,9 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
 BoundingBox SurfaceYPlane::bounding_box(bool pos_side) const
 {
   if (pos_side) {
-    return {-INFTY, INFTY, y0_, INFTY, -INFTY, INFTY};
+    return {{-INFTY, y0_, -INFTY}, {INFTY, INFTY, INFTY}};
   } else {
-    return {-INFTY, INFTY, -INFTY, y0_, -INFTY, INFTY};
+    return {{-INFTY, -INFTY, -INFTY}, {INFTY, y0_, INFTY}};
   }
 }
 
@@ -331,9 +331,9 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
 BoundingBox SurfaceZPlane::bounding_box(bool pos_side) const
 {
   if (pos_side) {
-    return {-INFTY, INFTY, -INFTY, INFTY, z0_, INFTY};
+    return {{-INFTY, -INFTY, z0_}, {INFTY, INFTY, INFTY}};
   } else {
-    return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, z0_};
+    return {{-INFTY, -INFTY, -INFTY}, {INFTY, INFTY, z0_}};
   }
 }
 
@@ -492,8 +492,8 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
 BoundingBox SurfaceXCylinder::bounding_box(bool pos_side) const
 {
   if (!pos_side) {
-    return {-INFTY, INFTY, y0_ - radius_, y0_ + radius_, z0_ - radius_,
-      z0_ + radius_};
+    return {{-INFTY, y0_ - radius_, z0_ - radius_},
+      {INFTY, y0_ + radius_, z0_ + radius_}};
   } else {
     return {};
   }
@@ -535,8 +535,8 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
 BoundingBox SurfaceYCylinder::bounding_box(bool pos_side) const
 {
   if (!pos_side) {
-    return {x0_ - radius_, x0_ + radius_, -INFTY, INFTY, z0_ - radius_,
-      z0_ + radius_};
+    return {{x0_ - radius_, -INFTY, z0_ - radius_},
+      {x0_ + radius_, INFTY, z0_ + radius_}};
   } else {
     return {};
   }
@@ -579,8 +579,8 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
 BoundingBox SurfaceZCylinder::bounding_box(bool pos_side) const
 {
   if (!pos_side) {
-    return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, -INFTY,
-      INFTY};
+    return {{x0_ - radius_, y0_ - radius_, -INFTY},
+      {x0_ + radius_, y0_ + radius_, INFTY}};
   } else {
     return {};
   }
@@ -657,8 +657,8 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const
 BoundingBox SurfaceSphere::bounding_box(bool pos_side) const
 {
   if (!pos_side) {
-    return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_,
-      z0_ - radius_, z0_ + radius_};
+    return {{x0_ - radius_, y0_ - radius_, z0_ - radius_},
+      {x0_ + radius_, y0_ + radius_, z0_ + radius_}};
   } else {
     return {};
   }
@@ -1334,8 +1334,44 @@ void read_surfaces(pugi::xml_node node)
       surf1.bc_ = make_unique<TranslationalPeriodicBC>(i_surf, j_surf);
       surf2.bc_ = make_unique<TranslationalPeriodicBC>(i_surf, j_surf);
     } else {
-      surf1.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf);
-      surf2.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf);
+      // check that both normals have at least one 0 component
+      if (std::abs(norm1.x) > FP_PRECISION &&
+          std::abs(norm1.y) > FP_PRECISION &&
+          std::abs(norm1.z) > FP_PRECISION) {
+        fatal_error(fmt::format(
+          "The normal ({}) of the periodic surface ({}) does not contain any "
+          "component with a zero value. A RotationalPeriodicBC requires one "
+          "component which is zero for both plane normals.",
+          norm1, i_surf));
+      }
+      if (std::abs(norm2.x) > FP_PRECISION &&
+          std::abs(norm2.y) > FP_PRECISION &&
+          std::abs(norm2.z) > FP_PRECISION) {
+        fatal_error(fmt::format(
+          "The normal ({}) of the periodic surface ({}) does not contain any "
+          "component with a zero value. A RotationalPeriodicBC requires one "
+          "component which is zero for both plane normals.",
+          norm2, j_surf));
+      }
+      // find common zero component, which indicates the periodic axis
+      RotationalPeriodicBC::PeriodicAxis axis;
+      if (std::abs(norm1.x) <= FP_PRECISION &&
+          std::abs(norm2.x) <= FP_PRECISION) {
+        axis = RotationalPeriodicBC::PeriodicAxis::x;
+      } else if (std::abs(norm1.y) <= FP_PRECISION &&
+                 std::abs(norm2.y) <= FP_PRECISION) {
+        axis = RotationalPeriodicBC::PeriodicAxis::y;
+      } else if (std::abs(norm1.z) <= FP_PRECISION &&
+                 std::abs(norm2.z) <= FP_PRECISION) {
+        axis = RotationalPeriodicBC::PeriodicAxis::z;
+      } else {
+        fatal_error(fmt::format(
+          "There is no component which is 0.0 in both normal vectors. This "
+          "indicates that the two planes are not periodic about the X, Y, or Z "
+          "axis, which is not supported."));
+      }
+      surf1.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf, axis);
+      surf2.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf, axis);
     }
 
     // If albedo data is present in albedo map, set the boundary albedo.

src/tallies/filter_mesh.cpp

@@ -6,6 +6,7 @@
 #include "openmc/constants.h"
 #include "openmc/error.h"
 #include "openmc/mesh.h"
+#include "openmc/position.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -30,6 +31,10 @@ void MeshFilter::from_xml(pugi::xml_node node)
   if (check_for_node(node, "translation")) {
     set_translation(get_node_array<double>(node, "translation"));
   }
+  // Read the rotation transform.
+  if (check_for_node(node, "rotation")) {
+    set_rotation(get_node_array<double>(node, "rotation"));
+  }
 }
 
 void MeshFilter::get_all_bins(
@@ -45,6 +50,12 @@ void MeshFilter::get_all_bins(
     last_r -= translation();
     r -= translation();
   }
+  // apply rotation if present
+  if (!rotation_.empty()) {
+    last_r = last_r.rotate(rotation_);
+    r = r.rotate(rotation_);
+    u = u.rotate(rotation_);
+  }
 
   if (estimator != TallyEstimator::TRACKLENGTH) {
     auto bin = model::meshes[mesh_]->get_bin(r);
@@ -65,6 +76,9 @@ void MeshFilter::to_statepoint(hid_t filter_group) const
   if (translated_) {
     write_dataset(filter_group, "translation", translation_);
   }
+  if (rotated_) {
+    write_dataset(filter_group, "rotation", rotation_);
+  }
 }
 
 std::string MeshFilter::text_label(int bin) const
@@ -93,6 +107,40 @@ void MeshFilter::set_translation(const double translation[3])
   this->set_translation({translation[0], translation[1], translation[2]});
 }
 
+void MeshFilter::set_rotation(const vector<double>& rot)
+{
+  rotated_ = true;
+
+  // Compute and store the inverse rotation matrix for the angles given.
+  rotation_.clear();
+  rotation_.reserve(rot.size() == 9 ? 9 : 12);
+  if (rot.size() == 3) {
+    double phi = -rot[0] * PI / 180.0;
+    double theta = -rot[1] * PI / 180.0;
+    double psi = -rot[2] * PI / 180.0;
+    rotation_.push_back(std::cos(theta) * std::cos(psi));
+    rotation_.push_back(-std::cos(phi) * std::sin(psi) +
+                        std::sin(phi) * std::sin(theta) * std::cos(psi));
+    rotation_.push_back(std::sin(phi) * std::sin(psi) +
+                        std::cos(phi) * std::sin(theta) * std::cos(psi));
+    rotation_.push_back(std::cos(theta) * std::sin(psi));
+    rotation_.push_back(std::cos(phi) * std::cos(psi) +
+                        std::sin(phi) * std::sin(theta) * std::sin(psi));
+    rotation_.push_back(-std::sin(phi) * std::cos(psi) +
+                        std::cos(phi) * std::sin(theta) * std::sin(psi));
+    rotation_.push_back(-std::sin(theta));
+    rotation_.push_back(std::sin(phi) * std::cos(theta));
+    rotation_.push_back(std::cos(phi) * std::cos(theta));
+
+    // When user specifies angles, write them at end of vector
+    rotation_.push_back(rot[0]);
+    rotation_.push_back(rot[1]);
+    rotation_.push_back(rot[2]);
+  } else {
+    std::copy(rot.begin(), rot.end(), std::back_inserter(rotation_));
+  }
+}
+
 //==============================================================================
 // C-API functions
 //==============================================================================
@@ -201,4 +249,48 @@ extern "C" int openmc_mesh_filter_set_translation(
   return 0;
 }
 
+//! Return the rotation matrix of a mesh filter
+extern "C" int openmc_mesh_filter_get_rotation(
+  int32_t index, double rot[], size_t* n)
+{
+  // Make sure this is a valid index to an allocated filter
+  if (int err = verify_filter(index))
+    return err;
+
+  // Check the filter type
+  const auto& filter = model::tally_filters[index];
+  if (filter->type() != FilterType::MESH) {
+    set_errmsg("Tried to get a rotation from a non-mesh filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+  // Get rotation from the mesh filter and set value
+  auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());
+  *n = mesh_filter->rotation().size();
+  std::memcpy(rot, mesh_filter->rotation().data(),
+    *n * sizeof(mesh_filter->rotation()[0]));
+  return 0;
+}
+
+//! Set the flattened rotation matrix of a mesh filter
+extern "C" int openmc_mesh_filter_set_rotation(
+  int32_t index, const double rot[], size_t rot_len)
+{
+  // Make sure this is a valid index to an allocated filter
+  if (int err = verify_filter(index))
+    return err;
+
+  const auto& filter = model::tally_filters[index];
+  // Check the filter type
+  if (filter->type() != FilterType::MESH) {
+    set_errmsg("Tried to set a rotation from a non-mesh filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+
+  // Get a pointer to the filter and downcast
+  auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());
+  std::vector<double> vec_rot(rot, rot + rot_len);
+  mesh_filter->set_rotation(vec_rot);
+  return 0;
+}
+
 } // namespace openmc

src/tallies/tally.cpp

@@ -215,7 +215,7 @@ Tally::Tally(pugi::xml_node node)
                       "number of inactive cycles.");
         }
         settings::ifp_on = true;
-      } else {
+      } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
         fatal_error(
           "Iterated Fission Probability can only be used in an eigenvalue "
           "calculation.");

src/universe.cpp

@@ -61,7 +61,7 @@ bool Universe::find_cell(GeometryState& p) const
 
 BoundingBox Universe::bounding_box() const
 {
-  BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
+  BoundingBox bbox = BoundingBox::inverted();
   if (cells_.size() == 0) {
     return {};
   } else {

src/weight_windows.cpp

@@ -933,7 +933,8 @@ void WeightWindowsGenerator::create_tally()
   for (const auto& f : model::tally_filters) {
     if (f->type() == FilterType::MESH) {
       const auto* mesh_filter = dynamic_cast<MeshFilter*>(f.get());
-      if (mesh_filter->mesh() == mesh_idx && !mesh_filter->translated()) {
+      if (mesh_filter->mesh() == mesh_idx && !mesh_filter->translated() &&
+          !mesh_filter->rotated()) {
         ww_tally->add_filter(f.get());
         found_mesh_filter = true;
         break;

tests/regression_tests/distribmat/test.py

@@ -73,7 +73,7 @@ class DistribmatTestHarness(PyAPITestHarness):
         # Plots
         ####################
 
-        plot1 = openmc.Plot(plot_id=1)
+        plot1 = openmc.SlicePlot(plot_id=1)
         plot1.basis = 'xy'
         plot1.color_by = 'cell'
         plot1.filename = 'cellplot'
@@ -81,7 +81,7 @@ class DistribmatTestHarness(PyAPITestHarness):
         plot1.width = (7, 7)
         plot1.pixels = (400, 400)
 
-        plot2 = openmc.Plot(plot_id=2)
+        plot2 = openmc.SlicePlot(plot_id=2)
         plot2.basis = 'xy'
         plot2.color_by = 'material'
         plot2.filename = 'matplot'

tests/regression_tests/filter_rotations/inputs_true.dat

@@ -0,0 +1,59 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="10.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+    </material>
+    <material id="2">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="Zr90" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="1 -2 3 -4 10 -9" universe="1"/>
+    <cell id="2" material="2" region="(-1 | 2 | -3 | 4) (5 -6 7 -8) 10 -9" universe="1"/>
+    <surface id="1" name="minimum x" type="x-plane" coeffs="-5.0"/>
+    <surface id="2" name="maximum x" type="x-plane" coeffs="5.0"/>
+    <surface id="3" name="minimum y" type="y-plane" coeffs="-5.0"/>
+    <surface id="4" name="maximum y" type="y-plane" coeffs="5.0"/>
+    <surface id="5" name="minimum x" type="x-plane" boundary="reflective" coeffs="-10.0"/>
+    <surface id="6" name="maximum x" type="x-plane" boundary="reflective" coeffs="10.0"/>
+    <surface id="7" name="minimum y" type="y-plane" boundary="reflective" coeffs="-10.0"/>
+    <surface id="8" name="maximum y" type="y-plane" boundary="reflective" coeffs="10.0"/>
+    <surface id="9" type="z-plane" boundary="vacuum" coeffs="10.0"/>
+    <surface id="10" type="z-plane" boundary="vacuum" coeffs="-10.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>5</batches>
+    <inactive>0</inactive>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>3 4 5</dimension>
+      <lower_left>-9 -9 -9</lower_left>
+      <upper_right>9 9 9</upper_right>
+    </mesh>
+    <mesh id="2">
+      <dimension>3 4 5</dimension>
+      <lower_left>-9 -9 -9</lower_left>
+      <upper_right>9 9 9</upper_right>
+    </mesh>
+    <filter id="1" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="mesh" rotation="0 0 10">
+      <bins>2</bins>
+    </filter>
+    <tally id="1">
+      <filters>1</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="2">
+      <filters>2</filters>
+      <scores>total</scores>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/filter_rotations/results_true.dat

@@ -0,0 +1,244 @@
+k-combined:
+7.729082E-01 3.775399E-02
+tally 1:
+5.296804E-02
+5.661701E-04
+8.356446E-02
+1.412139E-03
+5.041335E-02
+5.143568E-04
+1.299348E-01
+3.467618E-03
+3.929702E-01
+3.147038E-02
+1.379707E-01
+3.888484E-03
+1.405034E-01
+4.473799E-03
+3.785796E-01
+2.940585E-02
+1.422010E-01
+4.113723E-03
+5.647073E-02
+6.735251E-04
+7.911154E-02
+1.329137E-03
+5.160755E-02
+5.361448E-04
+6.669424E-02
+9.090832E-04
+1.008621E-01
+2.134534E-03
+6.808932E-02
+9.355993E-04
+1.873006E-01
+7.135961E-03
+6.221575E-01
+7.819842E-02
+1.856653E-01
+6.954762E-03
+2.014929E-01
+8.327845E-03
+5.853251E-01
+6.945708E-02
+1.709645E-01
+5.917124E-03
+7.214913E-02
+1.058962E-03
+1.027720E-01
+2.138475E-03
+6.099853E-02
+7.493941E-04
+6.892071E-02
+9.630680E-04
+1.035459E-01
+2.173883E-03
+6.973870E-02
+9.904237E-04
+2.125703E-01
+9.112659E-03
+9.012205E-01
+2.163546E-01
+2.066426E-01
+8.617414E-03
+2.258950E-01
+1.039607E-02
+9.476792E-01
+2.350708E-01
+2.225585E-01
+1.017898E-02
+7.111503E-02
+1.036847E-03
+1.117012E-01
+2.530040E-03
+6.870474E-02
+9.551035E-04
+5.738897E-02
+6.699030E-04
+9.522335E-02
+1.835769E-03
+6.570917E-02
+8.656870E-04
+1.945592E-01
+7.593336E-03
+5.514753E-01
+6.122981E-02
+2.144202E-01
+9.421739E-03
+1.971631E-01
+7.944046E-03
+6.088996E-01
+7.442954E-02
+1.965447E-01
+7.765628E-03
+7.005494E-02
+1.012891E-03
+1.010633E-01
+2.084095E-03
+6.145926E-02
+7.694351E-04
+4.999479E-02
+5.129164E-04
+7.238243E-02
+1.062921E-03
+4.902309E-02
+4.852193E-04
+1.324655E-01
+3.642431E-03
+3.305312E-01
+2.265726E-02
+1.332993E-01
+3.728385E-03
+1.547469E-01
+4.894837E-03
+3.625944E-01
+2.747313E-02
+1.435761E-01
+4.334405E-03
+5.789603E-02
+7.065383E-04
+7.589559E-02
+1.205386E-03
+5.210018E-02
+5.790843E-04
+tally 2:
+4.597932E-02
+4.246849E-04
+8.494738E-02
+1.467238E-03
+5.155072E-02
+5.373129E-04
+1.479395E-01
+4.468164E-03
+3.931843E-01
+3.141943E-02
+1.264081E-01
+3.243689E-03
+1.229742E-01
+3.276501E-03
+3.768427E-01
+2.927971E-02
+1.574908E-01
+5.079304E-03
+5.621044E-02
+6.877866E-04
+8.490616E-02
+1.510673E-03
+4.709600E-02
+4.576632E-04
+5.382674E-02
+5.890600E-04
+1.116560E-01
+2.599817E-03
+6.648634E-02
+8.862677E-04
+1.996861E-01
+8.137198E-03
+6.218616E-01
+7.805532E-02
+1.672379E-01
+5.667797E-03
+1.841275E-01
+6.936896E-03
+5.887874E-01
+7.031665E-02
+1.872574E-01
+7.086584E-03
+7.574698E-02
+1.160992E-03
+1.100173E-01
+2.453972E-03
+5.427866E-02
+5.955556E-04
+5.644318E-02
+6.503634E-04
+1.058628E-01
+2.254338E-03
+8.012794E-02
+1.297032E-03
+2.292757E-01
+1.064490E-02
+9.001874E-01
+2.153097E-01
+1.921874E-01
+7.537343E-03
+2.058642E-01
+8.537136E-03
+9.442653E-01
+2.359241E-01
+2.419900E-01
+1.202201E-02
+7.563876E-02
+1.169874E-03
+1.165428E-01
+2.750405E-03
+5.646354E-02
+6.408683E-04
+4.963879E-02
+4.975915E-04
+1.005478E-01
+2.046400E-03
+7.041191E-02
+1.002590E-03
+1.959123E-01
+7.727902E-03
+5.625596E-01
+6.366248E-02
+1.987313E-01
+8.108832E-03
+1.922194E-01
+7.602752E-03
+6.062527E-01
+7.384124E-02
+2.039506E-01
+8.379010E-03
+7.019905E-02
+1.034239E-03
+1.081383E-01
+2.344040E-03
+5.125142E-02
+5.430397E-04
+4.230495E-02
+3.666730E-04
+7.722275E-02
+1.208856E-03
+4.998937E-02
+5.049355E-04
+1.461632E-01
+4.456332E-03
+3.315459E-01
+2.288314E-02
+1.238273E-01
+3.205844E-03
+1.355317E-01
+3.754923E-03
+3.646520E-01
+2.766067E-02
+1.543013E-01
+5.044487E-03
+6.082173E-02
+7.636985E-04
+8.490569E-02
+1.506309E-03
+4.046126E-02
+3.495832E-04

tests/regression_tests/filter_rotations/test.py

@@ -0,0 +1,72 @@
+import numpy as np
+
+import openmc
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def model():
+
+    model = openmc.model.Model()
+
+    fuel = openmc.Material()
+    fuel.set_density('g/cm3', 10.0)
+    fuel.add_nuclide('U235', 1.0)
+    zr = openmc.Material()
+    zr.set_density('g/cm3', 1.0)
+    zr.add_nuclide('Zr90', 1.0)
+    model.materials.extend([fuel, zr])
+
+    box1 = openmc.model.RectangularPrism(10.0, 10.0)
+    box2 = openmc.model.RectangularPrism(20.0, 20.0, boundary_type='reflective')
+    top = openmc.ZPlane(z0=10.0, boundary_type='vacuum')
+    bottom = openmc.ZPlane(z0=-10.0, boundary_type='vacuum')
+    cell1 = openmc.Cell(fill=fuel, region=-box1 & +bottom & -top)
+    cell2 = openmc.Cell(fill=zr, region=+box1 & -box2 & +bottom & -top)
+    model.geometry = openmc.Geometry([cell1, cell2])
+
+    model.settings.batches = 5
+    model.settings.inactive = 0
+    model.settings.particles = 1000
+
+    rotation = np.array((0, 0, 10))
+
+    llc = np.array([-9, -9, -9])
+    urc = np.array([9, 9, 9])
+
+    mesh_dims = (3, 4, 5)
+
+    filters = []
+
+    # un-rotated meshes
+    reg_mesh = openmc.RegularMesh()
+    reg_mesh.dimension = mesh_dims
+    reg_mesh.lower_left = llc
+    reg_mesh.upper_right = urc
+
+    filters.append(openmc.MeshFilter(reg_mesh))
+
+    # rotated meshes
+    rotated_reg_mesh = openmc.RegularMesh()
+    rotated_reg_mesh.dimension = mesh_dims
+    rotated_reg_mesh.lower_left = llc
+    rotated_reg_mesh.upper_right = urc
+
+    filters.append(openmc.MeshFilter(rotated_reg_mesh))
+    filters[-1].rotation = rotation
+
+    # Create tallies
+    for f in filters:
+        tally = openmc.Tally()
+        tally.filters = [f]
+        tally.scores = ['total']
+        model.tallies.append(tally)
+
+    return model
+
+
+def test_filter_mesh_rotations(model):
+    harness = PyAPITestHarness('statepoint.5.h5', model)
+    harness.main()

tests/regression_tests/periodic_cyls/test.py

@@ -0,0 +1,91 @@
+import openmc
+import numpy as np
+import pytest
+from openmc.utility_funcs import change_directory
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def xcyl_model():
+    model = openmc.Model()
+    # Define materials
+    fuel = openmc.Material()
+    fuel.add_nuclide('U235', 0.2)
+    fuel.add_nuclide('U238', 0.8)
+    fuel.set_density('g/cc', 19.1)
+    model.materials = openmc.Materials([fuel])
+
+    # Define geometry
+    # finite cylinder
+    x_min = openmc.XPlane(x0=0.0, boundary_type='reflective')
+    x_max = openmc.XPlane(x0=20.0, boundary_type='reflective')
+    x_cyl = openmc.XCylinder(r=20.0,boundary_type='vacuum')
+    # slice cylinder for periodic BC
+    periodic_bounding_yplane = openmc.YPlane(y0=0, boundary_type='periodic')
+    periodic_bounding_plane = openmc.Plane(
+        a=0.0, b=-np.sqrt(3) / 3, c=1, boundary_type='periodic',
+    )
+    sixth_cyl_cell = openmc.Cell(1, fill=fuel, region = 
+        +x_min &- x_max & -x_cyl & +periodic_bounding_yplane & +periodic_bounding_plane)
+    periodic_bounding_yplane.periodic_surface = periodic_bounding_plane
+    periodic_bounding_plane.periodic_surface = periodic_bounding_yplane
+
+    model.geometry = openmc.Geometry([sixth_cyl_cell])
+
+
+    # Define settings
+    model.settings.particles = 1000
+    model.settings.batches = 4
+    model.settings.inactive = 0
+    model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
+        (0, 0, 0), (20, 20, 20))
+    )
+    return model
+
+@pytest.fixture
+def ycyl_model():
+    model = openmc.Model()
+    # Define materials
+    fuel = openmc.Material()
+    fuel.add_nuclide('U235', 0.2)
+    fuel.add_nuclide('U238', 0.8)
+    fuel.set_density('g/cc', 19.1)
+    model.materials = openmc.Materials([fuel])
+
+    # Define geometry
+    # finite cylinder
+    y_min = openmc.YPlane(y0=0.0, boundary_type='reflective')
+    y_max = openmc.YPlane(y0=20.0, boundary_type='reflective')
+    y_cyl = openmc.YCylinder(r=20.0,boundary_type='vacuum')
+    # slice cylinder for periodic BC
+    periodic_bounding_xplane = openmc.XPlane(x0=0, boundary_type='periodic')
+    periodic_bounding_plane = openmc.Plane(
+        a=-np.sqrt(3) / 3, b=0.0, c=1, boundary_type='periodic',
+    )
+    sixth_cyl_cell = openmc.Cell(1, fill=fuel, region = 
+        +y_min &- y_max & -y_cyl & +periodic_bounding_xplane & +periodic_bounding_plane)
+    periodic_bounding_xplane.periodic_surface = periodic_bounding_plane
+    periodic_bounding_plane.periodic_surface = periodic_bounding_xplane
+    model.geometry = openmc.Geometry([sixth_cyl_cell])
+
+
+    # Define settings
+    model.settings.particles = 1000
+    model.settings.batches = 4
+    model.settings.inactive = 0
+    model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
+        (0, 0, 0), (20, 20, 20))
+    )
+    return model
+
+def test_xcyl(xcyl_model):
+    with change_directory("xcyl_model"):
+        openmc.reset_auto_ids()
+        harness = PyAPITestHarness('statepoint.4.h5', xcyl_model)
+        harness.main()
+
+def test_ycyl(ycyl_model):
+    with change_directory("ycyl_model"):
+        openmc.reset_auto_ids()
+        harness = PyAPITestHarness('statepoint.4.h5', ycyl_model)
+        harness.main()

tests/regression_tests/periodic_cyls/xcyl_model/inputs_true.dat

@@ -0,0 +1,29 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="19.1" units="g/cc"/>
+      <nuclide name="U235" ao="0.2"/>
+      <nuclide name="U238" ao="0.8"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="1 -2 -3 4 5" universe="1"/>
+    <surface id="1" type="x-plane" boundary="reflective" coeffs="0.0"/>
+    <surface id="2" type="x-plane" boundary="reflective" coeffs="20.0"/>
+    <surface id="3" type="x-cylinder" boundary="vacuum" coeffs="0.0 0.0 20.0"/>
+    <surface id="4" type="y-plane" boundary="periodic" coeffs="0" periodic_surface_id="5"/>
+    <surface id="5" type="plane" boundary="periodic" coeffs="0.0 -0.5773502691896257 1 0.0" periodic_surface_id="4"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>4</batches>
+    <inactive>0</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>0 0 0 20 20 20</parameters>
+      </space>
+    </source>
+  </settings>
+  </model>

tests/regression_tests/periodic_cyls/xcyl_model/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+1.082283E+00 6.676373E-02

tests/regression_tests/periodic_cyls/ycyl_model/inputs_true.dat

@@ -0,0 +1,29 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="2" depletable="true">
+      <density value="19.1" units="g/cc"/>
+      <nuclide name="U235" ao="0.2"/>
+      <nuclide name="U238" ao="0.8"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="2" region="6 -7 -8 9 10" universe="2"/>
+    <surface id="6" type="y-plane" boundary="reflective" coeffs="0.0"/>
+    <surface id="7" type="y-plane" boundary="reflective" coeffs="20.0"/>
+    <surface id="8" type="y-cylinder" boundary="vacuum" coeffs="0.0 0.0 20.0"/>
+    <surface id="9" type="x-plane" boundary="periodic" coeffs="0" periodic_surface_id="10"/>
+    <surface id="10" type="plane" boundary="periodic" coeffs="-0.5773502691896257 0.0 1 0.0" periodic_surface_id="9"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>4</batches>
+    <inactive>0</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>0 0 0 20 20 20</parameters>
+      </space>
+    </source>
+  </settings>
+  </model>

tests/regression_tests/periodic_cyls/ycyl_model/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+1.082652E+00 3.316031E-02

tests/regression_tests/random_ray_adjoint_k_eff/results_true.dat

@@ -1,171 +1,171 @@
 k-combined:
 1.006640E+00 1.812969E-03
 tally 1:
-6.684129E+00
-8.939821E+00
-2.685967E+00
-1.443592E+00
+1.208044E+00
+2.920182E-01
+4.854426E-01
+4.715453E-02
 0.000000E+00
 0.000000E+00
-6.358774E+00
-8.091444E+00
-9.687217E-01
-1.878029E-01
+1.149242E+00
+2.643067E-01
+1.750801E-01
+6.134563E-03
 0.000000E+00
 0.000000E+00
-5.963160E+00
-7.117108E+00
-1.932332E-01
-7.473914E-03
+1.077743E+00
+2.324814E-01
+3.492371E-02
+2.441363E-04
 0.000000E+00
 0.000000E+00
-5.137593E+00
-5.283310E+00
-1.714616E-01
-5.884834E-03
-1.086218E-06
-2.361752E-13
-4.857253E+00
-4.719856E+00
-5.689580E-02
-6.476286E-04
-2.989356E-03
-1.787808E-06
-4.830516E+00
-4.666801E+00
-7.203015E-03
-1.037676E-05
-3.620020E+00
-2.620927E+00
-5.161382E+00
-5.328124E+00
-6.786255E-02
-9.210763E-04
-5.531943E+00
-6.120553E+00
-5.414034E+00
-5.864661E+00
+9.285362E-01
+1.725808E-01
+3.098889E-02
+1.922297E-04
+1.963161E-07
+7.714727E-15
+8.778641E-01
+1.541719E-01
+1.028293E-02
+2.115448E-05
+5.402741E-04
+5.839789E-08
+8.730274E-01
+1.524358E-01
+1.301813E-03
+3.389450E-07
+6.542525E-01
+8.560964E-02
+9.328247E-01
+1.740366E-01
+1.226491E-02
+3.008584E-05
+9.997969E-01
+1.999204E-01
+9.784958E-01
+1.915688E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-5.632338E+00
-6.347626E+00
+1.017952E+00
+2.073461E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-5.682608E+00
-6.462382E+00
+1.027039E+00
+2.110955E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-5.310716E+00
-5.645180E+00
+9.598240E-01
+1.844004E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-4.945409E+00
-4.893171E+00
+8.937969E-01
+1.598332E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-4.842688E+00
-4.690352E+00
+8.752275E-01
+1.532052E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-5.117198E+00
-5.237280E+00
+9.248400E-01
+1.710699E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-6.938711E+00
-9.633345E+00
-2.835258E+00
-1.608212E+00
+1.254054E+00
+3.146708E-01
+5.124223E-01
+5.253093E-02
 0.000000E+00
 0.000000E+00
-6.549505E+00
-8.584036E+00
-1.015138E+00
-2.061993E-01
+1.183712E+00
+2.803961E-01
+1.834683E-01
+6.735381E-03
 0.000000E+00
 0.000000E+00
-6.050651E+00
-7.327711E+00
-1.992816E-01
-7.948424E-03
+1.093555E+00
+2.393604E-01
+3.601678E-02
+2.596341E-04
 0.000000E+00
 0.000000E+00
-5.113981E+00
-5.234801E+00
-1.732323E-01
-6.006619E-03
-1.097435E-06
-2.410627E-13
-4.837033E+00
-4.680541E+00
-5.760042E-02
-6.637112E-04
-3.026377E-03
-1.832205E-06
-4.827049E+00
-4.660105E+00
-7.319913E-03
-1.071647E-05
-3.678770E+00
-2.706730E+00
-5.175337E+00
-5.356957E+00
-6.923046E-02
-9.586177E-04
-5.643451E+00
-6.370016E+00
-6.693323E+00
-8.964322E+00
-2.753307E+00
-1.516683E+00
+9.242694E-01
+1.709967E-01
+3.130894E-02
+1.962084E-04
+1.983437E-07
+7.874400E-15
+8.742100E-01
+1.528879E-01
+1.041026E-02
+2.167970E-05
+5.469644E-04
+5.984780E-08
+8.724009E-01
+1.522171E-01
+1.322938E-03
+3.500386E-07
+6.648691E-01
+8.841161E-02
+9.353464E-01
+1.749781E-01
+1.251209E-02
+3.131171E-05
+1.019947E+00
+2.080664E-01
+1.209708E+00
+2.928214E-01
+4.976146E-01
+4.954258E-02
 0.000000E+00
 0.000000E+00
-6.358384E+00
-8.090233E+00
-9.912008E-01
-1.965868E-01
+1.149174E+00
+2.642694E-01
+1.791431E-01
+6.421530E-03
 0.000000E+00
 0.000000E+00
-5.957484E+00
-7.103246E+00
-1.974033E-01
-7.798286E-03
+1.076718E+00
+2.320295E-01
+3.567737E-02
+2.547314E-04
 0.000000E+00
 0.000000E+00
-5.130744E+00
-5.268844E+00
-1.749233E-01
-6.123348E-03
-1.108148E-06
-2.457474E-13
-4.857340E+00
-4.720019E+00
-5.816659E-02
-6.768049E-04
-3.056125E-03
-1.868351E-06
-4.830629E+00
-4.667018E+00
-7.366289E-03
-1.085264E-05
-3.702077E+00
-2.741125E+00
-5.164864E+00
-5.335279E+00
-6.947917E-02
-9.655086E-04
-5.663725E+00
-6.415806E+00
+9.272977E-01
+1.721077E-01
+3.161443E-02
+2.000179E-04
+2.002789E-07
+8.027290E-15
+8.778794E-01
+1.541769E-01
+1.051257E-02
+2.210729E-05
+5.523400E-04
+6.102818E-08
+8.730477E-01
+1.524429E-01
+1.331320E-03
+3.544869E-07
+6.690816E-01
+8.953516E-02
+9.334544E-01
+1.742707E-01
+1.255707E-02
+3.153710E-05
+1.023613E+00
+2.095641E-01

tests/regression_tests/random_ray_auto_convert/infinite_medium/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-7.797820E-01 1.054725E-02
+7.479770E-01 1.624548E-02

tests/regression_tests/random_ray_auto_convert/material_wise/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-7.356667E-01 6.637270E-03
+7.372542E-01 6.967831E-03

tests/regression_tests/random_ray_auto_convert/stochastic_slab/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-7.551716E-01 8.117378E-03
+6.413334E-01 2.083132E-02

tests/regression_tests/random_ray_auto_convert/test.py

@@ -27,7 +27,7 @@ def test_random_ray_auto_convert(method):
 
         # Convert to a multi-group model
         model.convert_to_multigroup(
-            method=method, groups='CASMO-2', nparticles=30,
+            method=method, groups='CASMO-2', nparticles=100,
             overwrite_mgxs_library=False, mgxs_path="mgxs.h5"
         )
 

tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/model/inputs_true.dat

@@ -0,0 +1,61 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="UO2__2_4__" depletable="true">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="UO2__2_4__"/>
+    </material>
+    <material id="2" name="Zircaloy">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="Zircaloy"/>
+    </material>
+    <material id="3" name="Hot_borated_water">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="Hot_borated_water"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" name="Fuel" material="1" region="-1" universe="0"/>
+    <cell id="2" name="Cladding" material="2" region="1 -2" universe="0"/>
+    <cell id="3" name="Water" material="3" region="2 3 -4 5 -6" universe="0"/>
+    <surface id="1" name="Fuel OR" type="z-cylinder" coeffs="0 0 0.39218"/>
+    <surface id="2" name="Clad OR" type="z-cylinder" coeffs="0 0 0.4572"/>
+    <surface id="3" name="left" type="x-plane" boundary="reflective" coeffs="-0.63"/>
+    <surface id="4" name="right" type="x-plane" boundary="reflective" coeffs="0.63"/>
+    <surface id="5" name="bottom" type="y-plane" boundary="reflective" coeffs="-0.63"/>
+    <surface id="6" name="top" type="y-plane" boundary="reflective" coeffs="0.63"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <energy type="discrete">
+        <parameters>7000000.0 1.0</parameters>
+      </energy>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <source type="independent" strength="1.0" particle="neutron">
+        <space type="box">
+          <parameters>-0.63 -0.63 -1.0 0.63 0.63 1.0</parameters>
+        </space>
+      </source>
+      <distance_inactive>30.0</distance_inactive>
+      <distance_active>150.0</distance_active>
+      <source_region_meshes>
+        <mesh id="1">
+          <domain id="0" type="universe"/>
+        </mesh>
+      </source_region_meshes>
+      <source_shape>linear</source_shape>
+    </random_ray>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <upper_right>0.63 0.63</upper_right>
+    </mesh>
+  </settings>
+  </model>

tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/model/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+7.657815E-01 2.317564E-02

tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/user/inputs_true.dat

@@ -0,0 +1,64 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="UO2__2_4__" depletable="true">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="UO2__2_4__"/>
+    </material>
+    <material id="2" name="Zircaloy">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="Zircaloy"/>
+    </material>
+    <material id="3" name="Hot_borated_water">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="Hot_borated_water"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" name="Fuel" material="1" region="-1" universe="0"/>
+    <cell id="2" name="Cladding" material="2" region="1 -2" universe="0"/>
+    <cell id="3" name="Water" material="3" region="2 3 -4 5 -6" universe="0"/>
+    <surface id="1" name="Fuel OR" type="z-cylinder" coeffs="0 0 0.39218"/>
+    <surface id="2" name="Clad OR" type="z-cylinder" coeffs="0 0 0.4572"/>
+    <surface id="3" name="left" type="x-plane" boundary="reflective" coeffs="-0.63"/>
+    <surface id="4" name="right" type="x-plane" boundary="reflective" coeffs="0.63"/>
+    <surface id="5" name="bottom" type="y-plane" boundary="reflective" coeffs="-0.63"/>
+    <surface id="6" name="top" type="y-plane" boundary="reflective" coeffs="0.63"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <source type="independent" strength="1.0" particle="neutron">
+        <space type="box">
+          <parameters>-0.63 -0.63 -1.0 0.63 0.63 1.0</parameters>
+        </space>
+      </source>
+      <distance_inactive>30.0</distance_inactive>
+      <distance_active>150.0</distance_active>
+      <source_region_meshes>
+        <mesh id="1">
+          <domain id="0" type="universe"/>
+        </mesh>
+      </source_region_meshes>
+      <source_shape>linear</source_shape>
+    </random_ray>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <upper_right>0.63 0.63</upper_right>
+    </mesh>
+  </settings>
+  </model>

tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/user/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+7.827784E-01 2.062954E-02

tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/model/inputs_true.dat

@@ -0,0 +1,61 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="UO2__2_4__" depletable="true">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="UO2__2_4__"/>
+    </material>
+    <material id="2" name="Zircaloy">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="Zircaloy"/>
+    </material>
+    <material id="3" name="Hot_borated_water">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="Hot_borated_water"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" name="Fuel" material="1" region="-1" universe="0"/>
+    <cell id="2" name="Cladding" material="2" region="1 -2" universe="0"/>
+    <cell id="3" name="Water" material="3" region="2 3 -4 5 -6" universe="0"/>
+    <surface id="1" name="Fuel OR" type="z-cylinder" coeffs="0 0 0.39218"/>
+    <surface id="2" name="Clad OR" type="z-cylinder" coeffs="0 0 0.4572"/>
+    <surface id="3" name="left" type="x-plane" boundary="reflective" coeffs="-0.63"/>
+    <surface id="4" name="right" type="x-plane" boundary="reflective" coeffs="0.63"/>
+    <surface id="5" name="bottom" type="y-plane" boundary="reflective" coeffs="-0.63"/>
+    <surface id="6" name="top" type="y-plane" boundary="reflective" coeffs="0.63"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <energy type="discrete">
+        <parameters>7000000.0 1.0</parameters>
+      </energy>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <source type="independent" strength="1.0" particle="neutron">
+        <space type="box">
+          <parameters>-0.63 -0.63 -1.0 0.63 0.63 1.0</parameters>
+        </space>
+      </source>
+      <distance_inactive>30.0</distance_inactive>
+      <distance_active>150.0</distance_active>
+      <source_region_meshes>
+        <mesh id="1">
+          <domain id="0" type="universe"/>
+        </mesh>
+      </source_region_meshes>
+      <source_shape>linear</source_shape>
+    </random_ray>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <upper_right>0.63 0.63</upper_right>
+    </mesh>
+  </settings>
+  </model>

tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/model/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+7.479571E-01 2.398563E-02