fix typo

.github/workflows/ci.yml

@@ -1,4 +1,4 @@
-name: CI
+name: Tests and Coverage
 
 on:
   # allows us to run workflows manually
@@ -21,7 +21,25 @@ env:
   GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
 
 jobs:
+  filter-changes:
+    runs-on: ubuntu-latest
+    outputs:
+      source_changed: ${{ steps.filter.outputs.source_changed }}
+    steps:
+      - name: Check out the repository
+        uses: actions/checkout@v4
+      - name: Examine changed files
+        id: filter
+        uses: dorny/paths-filter@668c092af3649c4b664c54e4b704aa46782f6f7c # latest master commit, not released yet
+        with:
+          filters: |
+            source_changed:
+              - '!docs/**'
+              - '!**/*.md'
+          predicate-quantifier: 'every'
   main:
+    needs: filter-changes
+    if: ${{ needs.filter-changes.outputs.source_changed == 'true' }}
     runs-on: ubuntu-22.04
     strategy:
       matrix:
@@ -205,12 +223,33 @@ jobs:
           flag-name: C++ and Python
           path-to-lcov: coverage.lcov
 
-  finish:
-    needs: main
+  coverage:
+    needs: [filter-changes, main]
+    if: ${{ always() }}
     runs-on: ubuntu-latest
     steps:
       - name: Coveralls Finished
+        if: ${{ needs.filter-changes.outputs.source_changed == 'true' }}
         uses: coverallsapp/github-action@v2
         with:
           github-token: ${{ secrets.GITHUB_TOKEN }}
           parallel-finished: true
+
+  ci-pass:
+    needs: [filter-changes, main, coverage]
+    name: Check CI status
+    if: ${{ always() }}
+    runs-on: ubuntu-latest
+    steps:
+      - name: Check CI status
+        run: |
+          if [[ "${{ needs.filter-changes.outputs.source_changed }}" == "false" ]]; then
+            echo "Documentation-only change - CI skipped successfully"
+            exit 0
+          fi
+          if [[ "${{ needs.main.result }}" == "success" && "${{ needs.coverage.result }}" == "success" ]]; then
+            echo "CI passed"
+            exit 0
+          fi
+          echo "CI failed"
+          exit 1

.gitignore

@@ -30,7 +30,7 @@ docs/source/_images/*.aux
 docs/source/pythonapi/generated/
 
 # Source build
-build
+build*/
 
 # build from src/utils/setup.py
 src/utils/build

AGENTS.md

@@ -31,6 +31,24 @@ C++17/Python codebase where:
 - **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)
 
+## Git Branching Workflow
+
+OpenMC uses a git flow branching model with two primary branches:
+
+- **`develop` branch**: The main development branch where all ongoing development takes place. This is the **primary branch against which pull requests are submitted and merged**. This branch is not guaranteed to be stable and may contain work-in-progress features.
+- **`master` branch**: The stable release branch containing the latest stable release of OpenMC. This branch only receives merges from `develop` when the development team decides a release should occur.
+
+### Instructions for Code Review
+
+When analyzing code changes on a feature or bugfix branch (e.g., when a user asks "what do you think of these changes?"), **compare the branch changes against `develop`, not `master`**. Pull requests are submitted to merge into `develop`, so differences relative to `develop` represent the actual proposed changes. Comparing against `master` will include unrelated changes from other features that have already been merged to `develop`.
+
+### Workflow for contributors
+
+1. Create a feature/bugfix branch off `develop`
+2. Make changes and commit to the feature branch
+3. Open a pull request to merge the feature branch into `develop`
+4. A committer reviews and merges the PR into `develop`
+
 ## Critical Build & Test Workflows
 
 ### Build Dependencies

CMakeLists.txt

@@ -372,6 +372,7 @@ list(APPEND libopenmc_SOURCES
   src/particle.cpp
   src/particle_data.cpp
   src/particle_restart.cpp
+  src/particle_type.cpp
   src/photon.cpp
   src/physics.cpp
   src/physics_common.cpp

docs/source/io_formats/collision_track.rst

@@ -10,7 +10,7 @@ may also be written after each batch when multiple files are requested
 (``collision_track.N.h5``) or when the run is performed in parallel. The file
 contains the information needed to reconstruct each recorded collision.
 
-The current revision of the collision track file format is 1.0.
+The current revision of the collision track file format is 1.1.
 
 **/**
 
@@ -37,9 +37,9 @@ The current revision of the collision track file format is 1.0.
               - ``material_id`` (*int*) -- ID of the material containing the collision site.
               - ``universe_id`` (*int*) -- ID of the universe containing the collision site.
               - ``n_collision`` (*int*) -- Collision counter for the particle history.
-              - ``particle`` (*int*) -- Particle type (0=neutron, 1=photon, 2=electron, 3=positron).
-              - ``parent_id`` (*int64*) -- Unique ID of the parent particle.
-              - ``progeny_id`` (*int64*) -- Progeny ID of the particle.
+              - ``particle`` (*int32_t*) -- Particle type (PDG number).
+              - ``parent_id`` (*int64_t*) -- Unique ID of the parent particle.
+              - ``progeny_id`` (*int64_t*) -- Progeny ID of the particle.
 
 In an MPI run, OpenMC writes the combined dataset by gathering collision-track
 entries from all ranks before flushing them to disk, so the final file appears

docs/source/io_formats/depletion_results.rst

@@ -34,6 +34,7 @@ The current version of the depletion results file format is 1.2.
 
 :Attributes: - **index** (*int*) -- Index used in results for this material
              - **volume** (*double*) -- Volume of this material in [cm^3]
+             - **name** (*char[]*) -- Name of this material
 
 **/nuclides/<name>/**
 

docs/source/io_formats/particle_restart.rst

@@ -4,7 +4,7 @@
 Particle Restart File Format
 ============================
 
-The current version of the particle restart file format is 2.0.
+The current version of the particle restart file format is 2.1.
 
 **/**
 
@@ -26,8 +26,7 @@ The current version of the particle restart file format is 2.0.
            - **run_mode** (*char[]*) -- Run mode used, either 'fixed source',
              'eigenvalue', or 'particle restart'.
            - **id** (*int8_t*) -- Unique identifier of the particle.
-           - **type** (*int*) -- Particle type (0=neutron, 1=photon, 2=electron,
-             3=positron)
+           - **type** (*int32_t*) -- Particle type (PDG number)
            - **weight** (*double*) -- Weight of the particle.
            - **energy** (*double*) -- Energy of the particle in eV for
              continuous-energy mode, or the energy group of the particle for

docs/source/io_formats/settings.rst

@@ -721,7 +721,10 @@ attributes/sub-elements:
     is present.
 
   :particle:
-    The source particle type, either ``neutron`` or ``photon``.
+    The source particle type, specified as a PDG number or a string alias (e.g.,
+    ``neutron``/``n``, ``photon``/``gamma``, ``electron``, ``positron``,
+    ``proton``/``p``, ``deuteron``/``d``, ``triton``/``t``, ``alpha``, or GNDS
+    nuclide names like ``Fe57``).
 
     *Default*: neutron
 
@@ -1537,7 +1540,8 @@ sub-elements/attributes:
     *Default*: None
 
   :particle_type:
-    The particle that the weight windows will apply to (e.g., 'neutron')
+    The particle that the weight windows will apply to, specified as a PDG
+    code or string (e.g., ``neutron``).
 
     *Default*: 'neutron'
 
@@ -1597,7 +1601,8 @@ mesh-based weight windows.
     *Default*: None
 
   :particle_type:
-    The particle that the weight windows will apply to (e.g., 'neutron')
+    The particle that the weight windows will apply to, specified as a PDG
+    code or string (e.g., ``neutron``).
 
     *Default*: neutron
 

docs/source/io_formats/source.rst

@@ -15,6 +15,8 @@ following the same format.
 **/**
 
 :Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
+             - **version** (*int[2]*) -- Major and minor version of the source
+               file format.
 
 :Datasets:
 
@@ -22,5 +24,5 @@ following the same format.
              particle. The compound type has fields ``r``, ``u``, ``E``,
              ``time``, ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``,
              which represent the position, direction, energy, time, weight,
-             delayed group, surface ID, and particle type (0=neutron, 1=photon,
-             2=electron, 3=positron), respectively.
+             delayed group, surface ID, and particle type (PDG number),
+             respectively.

docs/source/io_formats/statepoint.rst

@@ -4,7 +4,7 @@
 State Point File Format
 =======================
 
-The current version of the statepoint file format is 18.1.
+The current version of the statepoint file format is 18.2.
 
 **/**
 
@@ -56,8 +56,8 @@ The current version of the statepoint file format is 18.1.
              ``time``, ``wgt``, ``delayed_group``, ``surf_id``, and
              ``particle``, which represent the position, direction, energy,
              time, weight, delayed group, surface ID, and particle type
-             (0=neutron, 1=photon, 2=electron, 3=positron), respectively. Only
-             present when `run_mode` is 'eigenvalue'.
+             (PDG number), respectively. Only present when `run_mode` is
+             'eigenvalue'.
 
 **/tallies/**
 

docs/source/io_formats/tallies.rst

@@ -318,8 +318,8 @@ should be set to:
   they use ``energy`` and ``y``.
 
 :particle:
-  A list of integers indicating the type of particles to tally ('neutron' = 1,
-  'photon' = 2, 'electron' = 3, 'positron' = 4).
+  A list of particle identifiers to tally, specified as strings (e.g.,
+  ``neutron``, ``photon``, ``He4``) or as integer PDG numbers.
 
 ------------------
 ``<mesh>`` Element

docs/source/io_formats/track.rst

@@ -4,7 +4,7 @@
 Track File Format
 =================
 
-The current revision of the particle track file format is 3.0.
+The current revision of the particle track file format is 3.1.
 
 **/**
 
@@ -32,6 +32,5 @@ The current revision of the particle track file format is 3.0.
                             the array for each primary/secondary particle. The
                             last offset should match the total size of the
                             array.
-                          - **particles** (*int[]*) -- Particle type for each
-                            primary/secondary particle (0=neutron, 1=photon,
-                            2=electron, 3=positron).
+                          - **particles** (*int32_t[]*) -- Particle type for
+                            each primary/secondary particle (PDG number).

docs/source/usersguide/random_ray.rst

@@ -513,6 +513,7 @@ Supported scores:
     - total
     - fission
     - nu-fission
+    - kappa-fission
     - events
 
 Supported Estimators:

docs/source/usersguide/settings.rst

@@ -400,7 +400,7 @@ below.
     {
       openmc::SourceSite particle;
       // weight
-      particle.particle = openmc::ParticleType::neutron;
+      particle.particle = openmc::ParticleType::neutron();
       particle.wgt = 1.0;
       // position
       double angle = 2.0 * M_PI * openmc::prn(seed);
@@ -477,7 +477,7 @@ parameters to the source class when it is created:
     {
       openmc::SourceSite particle;
       // weight
-      particle.particle = openmc::ParticleType::neutron;
+      particle.particle = openmc::ParticleType::neutron();
       particle.wgt = 1.0;
       // position
       particle.r.x = 0.0;

examples/custom_source/source_ring.cpp

@@ -1,17 +1,16 @@
-#include <cmath> // for M_PI
+#include <cmath>  // for M_PI
 #include <memory> // for unique_ptr
 
+#include "openmc/particle.h"
 #include "openmc/random_lcg.h"
 #include "openmc/source.h"
-#include "openmc/particle.h"
 
-class RingSource : public openmc::Source
-{
+class RingSource : public openmc::Source {
   openmc::SourceSite sample(uint64_t* seed) const
   {
     openmc::SourceSite particle;
     // particle type
-    particle.particle = openmc::ParticleType::neutron;
+    particle.particle = openmc::ParticleType::neutron();
     // position
     double angle = 2.0 * M_PI * openmc::prn(seed);
     double radius = 3.0;
@@ -25,10 +24,11 @@ class RingSource : public openmc::Source
   }
 };
 
-// A function to create a unique pointer to an instance of this class when generated
-// via a plugin call using dlopen/dlsym.
-// You must have external C linkage here otherwise dlopen will not find the file
-extern "C" std::unique_ptr<RingSource> openmc_create_source(std::string parameters)
+// A function to create a unique pointer to an instance of this class when
+// generated via a plugin call using dlopen/dlsym. You must have external C
+// linkage here otherwise dlopen will not find the file
+extern "C" std::unique_ptr<RingSource> openmc_create_source(
+  std::string parameters)
 {
   return std::make_unique<RingSource>();
 }

examples/parameterized_custom_source/parameterized_source_ring.cpp

@@ -1,63 +1,65 @@
-#include <cmath> // for M_PI
+#include <cmath>  // for M_PI
 #include <memory> // for unique_ptr
 #include <unordered_map>
 
+#include "openmc/particle.h"
 #include "openmc/random_lcg.h"
 #include "openmc/source.h"
-#include "openmc/particle.h"
 
 class RingSource : public openmc::Source {
-  public:
-    RingSource(double radius, double energy) : radius_(radius), energy_(energy) { }
+public:
+  RingSource(double radius, double energy) : radius_(radius), energy_(energy) {}
 
-    // Defines a function that can create a unique pointer to a new instance of this class
-    // by extracting the parameters from the provided string.
-    static std::unique_ptr<RingSource> from_string(std::string parameters)
-    {
-      std::unordered_map<std::string, std::string> parameter_mapping;
+  // Defines a function that can create a unique pointer to a new instance of
+  // this class by extracting the parameters from the provided string.
+  static std::unique_ptr<RingSource> from_string(std::string parameters)
+  {
+    std::unordered_map<std::string, std::string> parameter_mapping;
 
-      std::stringstream ss(parameters);
-      std::string parameter;
-      while (std::getline(ss, parameter, ',')) {
-        parameter.erase(0, parameter.find_first_not_of(' '));
-        std::string key = parameter.substr(0, parameter.find_first_of('='));
-        std::string value = parameter.substr(parameter.find_first_of('=') + 1, parameter.length());
-        parameter_mapping[key] = value;
-      }
-
-      double radius = std::stod(parameter_mapping["radius"]);
-      double energy = std::stod(parameter_mapping["energy"]);
-      return std::make_unique<RingSource>(radius, energy);
+    std::stringstream ss(parameters);
+    std::string parameter;
+    while (std::getline(ss, parameter, ',')) {
+      parameter.erase(0, parameter.find_first_not_of(' '));
+      std::string key = parameter.substr(0, parameter.find_first_of('='));
+      std::string value =
+        parameter.substr(parameter.find_first_of('=') + 1, parameter.length());
+      parameter_mapping[key] = value;
     }
 
-    // Samples from an instance of this class.
-    openmc::SourceSite sample(uint64_t* seed) const
-    {
-      openmc::SourceSite particle;
-      // particle type
-      particle.particle = openmc::ParticleType::neutron;
-      // position
-      double angle = 2.0 * M_PI * openmc::prn(seed);
-      double radius = this->radius_;
-      particle.r.x = radius * std::cos(angle);
-      particle.r.y = radius * std::sin(angle);
-      particle.r.z = 0.0;
-      // angle
-      particle.u = {1.0, 0.0, 0.0};
-      particle.E = this->energy_;
+    double radius = std::stod(parameter_mapping["radius"]);
+    double energy = std::stod(parameter_mapping["energy"]);
+    return std::make_unique<RingSource>(radius, energy);
+  }
 
-      return particle;
-    }
+  // Samples from an instance of this class.
+  openmc::SourceSite sample(uint64_t* seed) const
+  {
+    openmc::SourceSite particle;
+    // particle type
+    particle.particle = openmc::ParticleType::neutron();
+    // position
+    double angle = 2.0 * M_PI * openmc::prn(seed);
+    double radius = this->radius_;
+    particle.r.x = radius * std::cos(angle);
+    particle.r.y = radius * std::sin(angle);
+    particle.r.z = 0.0;
+    // angle
+    particle.u = {1.0, 0.0, 0.0};
+    particle.E = this->energy_;
 
-  private:
-    double radius_;
-    double energy_;
+    return particle;
+  }
+
+private:
+  double radius_;
+  double energy_;
 };
 
-// A function to create a unique pointer to an instance of this class when generated
-// via a plugin call using dlopen/dlsym.
-// You must have external C linkage here otherwise dlopen will not find the file
-extern "C" std::unique_ptr<RingSource> openmc_create_source(std::string parameters)
+// A function to create a unique pointer to an instance of this class when
+// generated via a plugin call using dlopen/dlsym. You must have external C
+// linkage here otherwise dlopen will not find the file
+extern "C" std::unique_ptr<RingSource> openmc_create_source(
+  std::string parameters)
 {
   return RingSource::from_string(parameters);
 }

include/openmc/capi.h

@@ -116,7 +116,7 @@ int openmc_mesh_set_id(int32_t index, int32_t id);
 int openmc_mesh_get_n_elements(int32_t index, size_t* n);
 int openmc_mesh_get_volumes(int32_t index, double* volumes);
 int openmc_mesh_material_volumes(int32_t index, int nx, int ny, int nz,
-  int max_mats, int32_t* materials, double* volumes);
+  int max_mats, int32_t* materials, double* volumes, double* bboxes);
 int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
 int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
 int openmc_new_filter(const char* type, int32_t* index);
@@ -201,8 +201,8 @@ int openmc_weight_windows_set_energy_bounds(
   int32_t index, double* e_bounds, size_t e_bounds_size);
 int openmc_weight_windows_get_energy_bounds(
   int32_t index, const double** e_bounds, size_t* e_bounds_size);
-int openmc_weight_windows_set_particle(int32_t index, int particle);
-int openmc_weight_windows_get_particle(int32_t index, int* particle);
+int openmc_weight_windows_set_particle(int32_t index, int32_t particle);
+int openmc_weight_windows_get_particle(int32_t index, int32_t* particle);
 int openmc_weight_windows_get_bounds(int32_t index, const double** lower_bounds,
   const double** upper_bounds, size_t* size);
 int openmc_weight_windows_set_bounds(int32_t index, const double* lower_bounds,
@@ -227,7 +227,7 @@ int openmc_zernike_filter_set_order(int32_t index, int order);
 int openmc_zernike_filter_set_params(
   int32_t index, const double* x, const double* y, const double* r);
 
-int openmc_particle_filter_get_bins(int32_t idx, int bins[]);
+int openmc_particle_filter_get_bins(int32_t idx, int32_t bins[]);
 
 //! Sets the mesh and energy grid for CMFD reweight
 //! \param[in] meshtyally_id id of CMFD Mesh Tally

include/openmc/cell.h

@@ -123,11 +123,11 @@ private:
   //! BoundingBox if the particle is in a complex cell.
   BoundingBox bounding_box_complex(vector<int32_t> postfix) const;
 
-  //! Enfource precedence: Parenthases, Complement, Intersection, Union
-  void add_precedence();
+  //! Enforce precedence between intersections and unions
+  void enforce_precedence();
 
   //! Add parenthesis to enforce precedence
-  int64_t add_parentheses(int64_t start);
+  void add_parentheses(int64_t start);
 
   //! Remove complement operators from the expression
   void remove_complement_ops();

include/openmc/constants.h

@@ -25,16 +25,16 @@ using double_4dvec = vector<vector<vector<vector<double>>>>;
 constexpr int HDF5_VERSION[] {3, 0};
 
 // Version numbers for binary files
-constexpr array<int, 2> VERSION_STATEPOINT {18, 1};
-constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
-constexpr array<int, 2> VERSION_TRACK {3, 0};
+constexpr array<int, 2> VERSION_STATEPOINT {18, 2};
+constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 1};
+constexpr array<int, 2> VERSION_TRACK {3, 1};
 constexpr array<int, 2> VERSION_SUMMARY {6, 1};
 constexpr array<int, 2> VERSION_VOLUME {1, 0};
 constexpr array<int, 2> VERSION_VOXEL {2, 0};
 constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
 constexpr array<int, 2> VERSION_PROPERTIES {1, 1};
 constexpr array<int, 2> VERSION_WEIGHT_WINDOWS {1, 0};
-constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 0};
+constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 1};
 
 // ============================================================================
 // ADJUSTABLE PARAMETERS

include/openmc/distribution_spatial.h

@@ -167,6 +167,7 @@ private:
 class SpatialBox : public SpatialDistribution {
 public:
   explicit SpatialBox(pugi::xml_node node, bool fission = false);
+  SpatialBox(Position lower_left, Position upper_right, bool fission = false);
 
   //! Sample a position from the distribution
   //! \param seed Pseudorandom number seed pointer

include/openmc/material.h

@@ -14,6 +14,7 @@
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/ncrystal_interface.h"
 #include "openmc/particle.h"
+#include "openmc/settings.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -110,9 +111,12 @@ public:
   //! \return Density in [atom/b-cm]
   double density() const { return density_; }
 
-  //! Get density in [g/cm^3]
+  //! Get density in [g/cm^3].
   //! \return Density in [g/cm^3]
-  double density_gpcc() const { return density_gpcc_; }
+  double density_gpcc() const
+  {
+    return settings::run_CE ? density_gpcc_ : density();
+  }
 
   //! Get charge density in [e/b-cm]
   //! \return Charge density in [e/b-cm]

include/openmc/mesh.h

@@ -4,6 +4,7 @@
 #ifndef OPENMC_MESH_H
 #define OPENMC_MESH_H
 
+#include <set>
 #include <unordered_map>
 
 #include "hdf5.h"
@@ -87,8 +88,12 @@ namespace detail {
 
 class MaterialVolumes {
 public:
+  MaterialVolumes(int32_t* mats, double* vols, double* bboxes, int table_size)
+    : materials_(mats), volumes_(vols), bboxes_(bboxes), table_size_(table_size)
+  {}
+
   MaterialVolumes(int32_t* mats, double* vols, int table_size)
-    : materials_(mats), volumes_(vols), table_size_(table_size)
+    : MaterialVolumes(mats, vols, nullptr, table_size)
   {}
 
   //! Add volume for a given material in a mesh element
@@ -96,8 +101,11 @@ public:
   //! \param[in] index_elem Index of the mesh element
   //! \param[in] index_material Index of the material within the model
   //! \param[in] volume Volume to add
-  void add_volume(int index_elem, int index_material, double volume);
-  void add_volume_unsafe(int index_elem, int index_material, double volume);
+  //! \param[in] bbox Bounding box to union into the result (optional)
+  void add_volume(int index_elem, int index_material, double volume,
+    const BoundingBox* bbox = nullptr);
+  void add_volume_unsafe(int index_elem, int index_material, double volume,
+    const BoundingBox* bbox = nullptr);
 
   // Accessors
   int32_t& materials(int i, int j) { return materials_[i * table_size_ + j]; }
@@ -112,11 +120,23 @@ public:
     return volumes_[i * table_size_ + j];
   }
 
+  double& bboxes(int i, int j, int k)
+  {
+    return bboxes_[(i * table_size_ + j) * 6 + k];
+  }
+  const double& bboxes(int i, int j, int k) const
+  {
+    return bboxes_[(i * table_size_ + j) * 6 + k];
+  }
+
+  bool has_bboxes() const { return bboxes_ != nullptr; }
+
   bool table_full() const { return table_full_; }
 
 private:
   int32_t* materials_;      //!< material index (bins, table_size)
   double* volumes_;         //!< volume in [cm^3] (bins, table_size)
+  double* bboxes_;          //!< bounding boxes (bins, table_size, 6)
   int table_size_;          //!< Size of hash table for each mesh element
   bool table_full_ {false}; //!< Whether the hash table is full
 };
@@ -239,6 +259,19 @@ public:
   void material_volumes(int nx, int ny, int nz, int max_materials,
     int32_t* materials, double* volumes) const;
 
+  //! Determine volume and bounding boxes of materials within each mesh element
+  //
+  //! \param[in] nx Number of samples in x direction
+  //! \param[in] ny Number of samples in y direction
+  //! \param[in] nz Number of samples in z direction
+  //! \param[in] max_materials Maximum number of materials in a single mesh
+  //!                          element
+  //! \param[inout] materials Array storing material indices
+  //! \param[inout] volumes Array storing volumes
+  //! \param[inout] bboxes Array storing bounding boxes (n_elems, table_size, 6)
+  void material_volumes(int nx, int ny, int nz, int max_materials,
+    int32_t* materials, double* volumes, double* bboxes) const;
+
   //! Determine bounding box of mesh
   //
   //! \return Bounding box of mesh
@@ -969,7 +1002,7 @@ public:
 
   Position sample_element(int32_t bin, uint64_t* seed) const override;
 
-  int get_bin(Position r) const override;
+  virtual int get_bin(Position r) const override;
 
   int n_bins() const override;
 
@@ -1007,16 +1040,21 @@ public:
 
 protected:
   // Methods
-
   //! Translate a bin value to an element reference
   virtual const libMesh::Elem& get_element_from_bin(int bin) const;
 
   //! Translate an element pointer to a bin index
   virtual int get_bin_from_element(const libMesh::Elem* elem) const;
 
+  // Data members
   libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set
                          //!< during intialization
+  vector<unique_ptr<libMesh::PointLocatorBase>>
+    pl_;                      //!< per-thread point locators
+  libMesh::BoundingBox bbox_; //!< bounding box of the mesh
+
 private:
+  // Methods
   void initialize() override;
   void set_mesh_pointer_from_filename(const std::string& filename);
   void build_eqn_sys();
@@ -1025,8 +1063,6 @@ private:
   unique_ptr<libMesh::MeshBase> unique_m_ =
     nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is
              //!< created inside OpenMC
-  vector<unique_ptr<libMesh::PointLocatorBase>>
-    pl_; //!< per-thread point locators
   unique_ptr<libMesh::EquationSystems>
     equation_systems_; //!< pointer to the libMesh EquationSystems
                        //!< instance
@@ -1035,7 +1071,6 @@ private:
   std::unordered_map<std::string, unsigned int>
     variable_map_; //!< mapping of variable names (tally scores) to libMesh
                    //!< variable numbers
-  libMesh::BoundingBox bbox_; //!< bounding box of the mesh
   libMesh::dof_id_type
     first_element_id_; //!< id of the first element in the mesh
 };
@@ -1043,8 +1078,9 @@ private:
 class AdaptiveLibMesh : public LibMesh {
 public:
   // Constructor
-  AdaptiveLibMesh(
-    libMesh::MeshBase& input_mesh, double length_multiplier = 1.0);
+  AdaptiveLibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0,
+    const std::set<libMesh::subdomain_id_type>& block_ids =
+      std::set<libMesh::subdomain_id_type>());
 
   // Overridden methods
   int n_bins() const override;
@@ -1056,6 +1092,8 @@ public:
 
   void write(const std::string& filename) const override;
 
+  int get_bin(Position r) const override;
+
 protected:
   // Overridden methods
   int get_bin_from_element(const libMesh::Elem* elem) const override;
@@ -1064,6 +1102,9 @@ protected:
 
 private:
   // Data members
+  const std::set<libMesh::subdomain_id_type>
+    block_ids_;               //!< subdomains of the mesh to tally on
+  const bool block_restrict_; //!< whether a subset of the mesh is being used
   const libMesh::dof_id_type num_active_; //!< cached number of active elements
 
   std::vector<libMesh::dof_id_type>

include/openmc/nuclide.h

@@ -163,7 +163,7 @@ bool multipole_in_range(const Nuclide& nuc, double E);
 namespace data {
 
 // Minimum/maximum transport energy for each particle type. Order corresponds to
-// that of the ParticleType enum
+// transport_index() for supported transport particles.
 extern array<double, 4> energy_min;
 extern array<double, 4> energy_max;
 

include/openmc/particle.h

@@ -126,10 +126,6 @@ public:
 //! Functions
 //============================================================================
 
-std::string particle_type_to_str(ParticleType type);
-
-ParticleType str_to_particle_type(std::string str);
-
 void add_surf_source_to_bank(Particle& p, const Surface& surf);
 
 } // namespace openmc

include/openmc/particle_data.h

@@ -3,6 +3,7 @@
 
 #include "openmc/array.h"
 #include "openmc/constants.h"
+#include "openmc/particle_type.h"
 #include "openmc/position.h"
 #include "openmc/random_lcg.h"
 #include "openmc/tallies/filter_match.h"
@@ -30,9 +31,6 @@ constexpr double CACHE_INVALID {-1.0};
 //==========================================================================
 // Aliases and type definitions
 
-//! Particle types
-enum class ParticleType { neutron, photon, electron, positron };
-
 //! Saved ("banked") state of a particle
 //! NOTE: This structure's MPI type is built in initialize_mpi() of
 //! initialize.cpp. Any changes made to the struct here must also be
@@ -496,7 +494,7 @@ private:
   MacroXS macro_xs_;
   CacheDataMG mg_xs_cache_;
 
-  ParticleType type_ {ParticleType::neutron};
+  ParticleType type_;
 
   double E_;
   double E_last_;

include/openmc/particle_type.h

@@ -0,0 +1,172 @@
+//==============================================================================
+// ParticleType class definition
+//==============================================================================
+
+#ifndef OPENMC_PARTICLE_TYPE_H
+#define OPENMC_PARTICLE_TYPE_H
+
+#include <cstddef>
+#include <cstdint>
+#include <string>
+#include <string_view>
+#include <type_traits>
+
+#include "openmc/constants.h"
+
+namespace openmc {
+
+//------------------------------------------------------------------------------
+// PDG constants (canonical particle identity as simple integers)
+//------------------------------------------------------------------------------
+
+inline constexpr int32_t PDG_NEUTRON = 2112;
+inline constexpr int32_t PDG_PHOTON = 22;
+inline constexpr int32_t PDG_ELECTRON = 11;
+inline constexpr int32_t PDG_POSITRON = -11;
+inline constexpr int32_t PDG_PROTON = 2212;
+inline constexpr int32_t PDG_DEUTERON = 1000010020;
+inline constexpr int32_t PDG_TRITON = 1000010030;
+inline constexpr int32_t PDG_ALPHA = 1000020040;
+
+//------------------------------------------------------------------------------
+// ParticleType class (standard-layout, trivially copyable)
+//------------------------------------------------------------------------------
+
+class ParticleType {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors
+
+  // Default constructor: defaults to neutron
+  constexpr ParticleType() : pdg_number_(PDG_NEUTRON) {}
+
+  // Constructor from PDG number
+  constexpr explicit ParticleType(int32_t pdg_number) : pdg_number_(pdg_number)
+  {}
+
+  // Constructor from particle name string (e.g., "neutron", "photon", "Fe56")
+  explicit ParticleType(std::string_view str);
+
+  // Constructor from Z, A, and metastable state for nuclear particles
+  constexpr ParticleType(int Z, int A, int m = 0)
+    : pdg_number_(1000000000 + Z * 10000 + A * 10 + m)
+  {}
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  // Accessor for the underlying PDG number
+  constexpr int32_t pdg_number() const { return pdg_number_; }
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  // Convert to string representation
+  std::string str() const;
+
+  // Check if this represents a nucleus (vs elementary particle)
+  constexpr bool is_nucleus() const
+  {
+    // PDG nuclear codes are >= 1000000000 (100ZZZAAAI format)
+    return pdg_number_ >= 1000000000;
+  }
+
+  // Get transport index (0-3 for transportable particles, C_NONE otherwise)
+  constexpr int transport_index() const;
+
+  // Check if this is a neutron
+  constexpr bool is_neutron() const { return pdg_number_ == PDG_NEUTRON; }
+
+  // Check if this is a photon
+  constexpr bool is_photon() const { return pdg_number_ == PDG_PHOTON; }
+
+  constexpr bool is_transportable() const
+  {
+    return this->transport_index() != C_NONE;
+  }
+
+  //----------------------------------------------------------------------------
+  // Static factory methods
+
+  static constexpr ParticleType neutron() { return ParticleType {PDG_NEUTRON}; }
+  static constexpr ParticleType photon() { return ParticleType {PDG_PHOTON}; }
+  static constexpr ParticleType electron()
+  {
+    return ParticleType {PDG_ELECTRON};
+  }
+  static constexpr ParticleType positron()
+  {
+    return ParticleType {PDG_POSITRON};
+  }
+  static constexpr ParticleType proton() { return ParticleType {PDG_PROTON}; }
+  static constexpr ParticleType deuteron()
+  {
+    return ParticleType {PDG_DEUTERON};
+  }
+  static constexpr ParticleType triton() { return ParticleType {PDG_TRITON}; }
+  static constexpr ParticleType alpha() { return ParticleType {PDG_ALPHA}; }
+
+private:
+  int32_t pdg_number_;
+};
+
+//------------------------------------------------------------------------------
+// Static assertions to ensure standard-layout and trivially copyable
+//------------------------------------------------------------------------------
+
+static_assert(std::is_standard_layout_v<ParticleType>,
+  "ParticleType must be standard-layout");
+static_assert(std::is_trivially_copyable_v<ParticleType>,
+  "ParticleType must be trivially copyable");
+static_assert(sizeof(ParticleType) == sizeof(int32_t),
+  "ParticleType must be same size as int32_t");
+
+//------------------------------------------------------------------------------
+// Comparison operators (free functions for symmetry)
+//------------------------------------------------------------------------------
+
+constexpr bool operator==(ParticleType lhs, ParticleType rhs)
+{
+  return lhs.pdg_number() == rhs.pdg_number();
+}
+
+constexpr bool operator!=(ParticleType lhs, ParticleType rhs)
+{
+  return lhs.pdg_number() != rhs.pdg_number();
+}
+
+constexpr bool operator<(ParticleType lhs, ParticleType rhs)
+{
+  return lhs.pdg_number() < rhs.pdg_number();
+}
+
+//------------------------------------------------------------------------------
+// ParticleType member function implementations (inline)
+//------------------------------------------------------------------------------
+
+constexpr int ParticleType::transport_index() const
+{
+  switch (pdg_number_) {
+  case PDG_NEUTRON:
+    return 0;
+  case PDG_PHOTON:
+    return 1;
+  case PDG_ELECTRON:
+    return 2;
+  case PDG_POSITRON:
+    return 3;
+  default:
+    return C_NONE;
+  }
+}
+
+//------------------------------------------------------------------------------
+// Legacy conversion helpers
+//------------------------------------------------------------------------------
+
+// Legacy enum code (0..3) to ParticleType conversion
+ParticleType legacy_particle_index_to_type(int code);
+
+} // namespace openmc
+
+#endif // OPENMC_PARTICLE_TYPE_H

include/openmc/random_ray/flat_source_domain.h

@@ -107,6 +107,7 @@ public:
   vector<double> nu_sigma_f_;
   vector<double> sigma_f_;
   vector<double> chi_;
+  vector<double> kappa_fission_;
 
   // 3D arrays stored in 1D representing values for all materials x energy
   // groups x energy groups

include/openmc/random_ray/random_ray_simulation.h

@@ -19,9 +19,9 @@ public:
 
   //----------------------------------------------------------------------------
   // Methods
-  void compute_segment_correction_factors();
   void apply_fixed_sources_and_mesh_domains();
   void prepare_fixed_sources_adjoint();
+  void prepare_adjoint_simulation();
   void simulate();
   void output_simulation_results() const;
   void instability_check(
@@ -34,9 +34,15 @@ public:
   // Accessors
   FlatSourceDomain* domain() const { return domain_.get(); }
 
+  //----------------------------------------------------------------------------
+  // Public data members
+
+  // Flag for adjoint simulation;
+  bool adjoint_needed_;
+
 private:
   //----------------------------------------------------------------------------
-  // Data members
+  // Private data members
 
   // Contains all flat source region data
   unique_ptr<FlatSourceDomain> domain_;
@@ -51,6 +57,9 @@ private:
   // Number of energy groups
   int negroups_;
 
+  // Toggle for first simulation
+  bool is_first_simulation_;
+
 }; // class RandomRaySimulation
 
 //============================================================================
@@ -60,6 +69,7 @@ private:
 void openmc_run_random_ray();
 void validate_random_ray_inputs();
 void openmc_reset_random_ray();
+void print_adjoint_header();
 
 } // namespace openmc
 

include/openmc/random_ray/source_region.h

@@ -146,6 +146,7 @@ public:
 
   // Scalar fields
   int* material_;
+  double* density_mult_;
   int* is_small_;
   int* n_hits_;
   int* birthday_;
@@ -195,6 +196,9 @@ public:
   int& material() { return *material_; }
   const int material() const { return *material_; }
 
+  double& density_mult() { return *density_mult_; }
+  const double density_mult() const { return *density_mult_; }
+
   int& is_small() { return *is_small_; }
   const int is_small() const { return *is_small_; }
 
@@ -316,7 +320,9 @@ public:
   //---------------------------------------
   // Scalar fields
 
-  int material_ {0}; //!< Index in openmc::model::materials array
+  int material_ {0};          //!< Index in openmc::model::materials array
+  double density_mult_ {1.0}; //!< A density multiplier queried from the cell
+                              //!< corresponding to the source region.
   OpenMPMutex lock_;
   double volume_ {
     0.0}; //!< Volume (computed from the sum of ray crossing lengths)
@@ -394,6 +400,9 @@ public:
   int& material(int64_t sr) { return material_[sr]; }
   const int material(int64_t sr) const { return material_[sr]; }
 
+  double& density_mult(int64_t sr) { return density_mult_[sr]; }
+  const double density_mult(int64_t sr) const { return density_mult_[sr]; }
+
   int& is_small(int64_t sr) { return is_small_[sr]; }
   const int is_small(int64_t sr) const { return is_small_[sr]; }
 
@@ -625,6 +634,7 @@ private:
 
   // SoA storage for scalar fields (one item per source region)
   vector<int> material_;
+  vector<double> density_mult_;
   vector<int> is_small_;
   vector<int> n_hits_;
   vector<int> mesh_;

include/openmc/reaction_product.h

@@ -10,7 +10,7 @@
 #include "openmc/chain.h"
 #include "openmc/endf.h"
 #include "openmc/memory.h" // for unique_ptr
-#include "openmc/particle.h"
+#include "openmc/particle_type.h"
 #include "openmc/vector.h" // for vector
 
 namespace openmc {

include/openmc/source.h

@@ -12,7 +12,7 @@
 #include "openmc/distribution_multi.h"
 #include "openmc/distribution_spatial.h"
 #include "openmc/memory.h"
-#include "openmc/particle.h"
+#include "openmc/particle_type.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -148,11 +148,11 @@ protected:
 
 private:
   // Data members
-  ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted
-  UPtrSpace space_;                               //!< Spatial distribution
-  UPtrAngle angle_;                               //!< Angular distribution
-  UPtrDist energy_;                               //!< Energy distribution
-  UPtrDist time_;                                 //!< Time distribution
+  ParticleType particle_; //!< Type of particle emitted
+  UPtrSpace space_;       //!< Spatial distribution
+  UPtrAngle angle_;       //!< Angular distribution
+  UPtrDist energy_;       //!< Energy distribution
+  UPtrDist time_;         //!< Time distribution
 };
 
 //==============================================================================

include/openmc/tallies/filter_particle.h

@@ -1,7 +1,7 @@
 #ifndef OPENMC_TALLIES_FILTER_PARTICLE_H
 #define OPENMC_TALLIES_FILTER_PARTICLE_H
 
-#include "openmc/particle.h"
+#include "openmc/particle_type.h"
 #include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"

include/openmc/weight_windows.h

@@ -10,7 +10,7 @@
 #include "openmc/constants.h"
 #include "openmc/memory.h"
 #include "openmc/mesh.h"
-#include "openmc/particle.h"
+#include "openmc/particle_type.h"
 #include "openmc/span.h"
 #include "openmc/tallies/tally.h"
 #include "openmc/vector.h"
@@ -193,10 +193,9 @@ public:
 private:
   //----------------------------------------------------------------------------
   // Data members
-  int32_t id_;    //!< Unique ID
-  int64_t index_; //!< Index into weight windows vector
-  ParticleType particle_type_ {
-    ParticleType::neutron};      //!< Particle type to apply weight windows to
+  int32_t id_;                   //!< Unique ID
+  int64_t index_;                //!< Index into weight windows vector
+  ParticleType particle_type_;   //!< Particle type to apply weight windows to
   vector<double> energy_bounds_; //!< Energy boundaries [eV]
   xt::xtensor<double, 2> lower_ww_; //!< Lower weight window bounds (shape:
                                     //!< energy_bins, mesh_bins (k, j, i))

openmc/data/photon.py

@@ -203,7 +203,9 @@ class AtomicRelaxation(EqualityMixin):
         for subshell, df in transitions.items():
             cv.check_value('subshell', subshell, _SUBSHELLS)
             cv.check_type('transitions', df, pd.DataFrame)
-        self._transitions = transitions
+        self._transitions = {
+            subshell: df.convert_dtypes() for subshell, df in transitions.items()
+        }
 
     @classmethod
     def from_ace(cls, ace):

openmc/deplete/abc.py

@@ -209,6 +209,8 @@ class TransportOperator(ABC):
             simulation.
         full_burn_list : list of int
             All burnable materials in the geometry.
+        name_list : list of str
+            Material names corresponding to materials in burn_list
         """
 
     def finalize(self):

openmc/deplete/openmc_operator.py

@@ -211,7 +211,7 @@ class OpenMCOperator(TransportOperator):
                            "section data.")
                     warn(msg)
             if mat.depletable:
-                burnable_mats.add(str(mat.id))
+                burnable_mats.add((str(mat.id), mat.name))
                 if mat.volume is None:
                     if mat.name is None:
                         msg = ("Volume not specified for depletable material "
@@ -229,9 +229,12 @@ class OpenMCOperator(TransportOperator):
                 "No depletable materials were found in the model.")
 
         # Sort the sets
-        burnable_mats = sorted(burnable_mats, key=int)
+        burnable_mats = sorted(burnable_mats, key=lambda x: int(x[0]))
         model_nuclides = sorted(model_nuclides)
 
+        # Store material names for later use
+        burnable_mats, self.name_list = zip(*burnable_mats)
+
         # Construct a global nuclide dictionary, burned first
         nuclides = list(self.chain.nuclide_dict)
         for nuc in model_nuclides:
@@ -541,6 +544,8 @@ class OpenMCOperator(TransportOperator):
             A list of all material IDs to be burned.  Used for sorting the simulation.
         full_burn_list : list
             List of all burnable material IDs
+        name_list : list of str
+            Material names corresponding to materials in burn_list
 
         """
         nuc_list = self.number.burnable_nuclides
@@ -554,4 +559,4 @@ class OpenMCOperator(TransportOperator):
         volume_list = comm.allgather(volume)
         volume = {k: v for d in volume_list for k, v in d.items()}
 
-        return volume, nuc_list, burn_list, self.burnable_mats
+        return volume, nuc_list, burn_list, self.burnable_mats, self.name_list

openmc/deplete/r2s.py

@@ -269,6 +269,7 @@ class R2SManager:
             # Compute material volume fractions on the mesh
             if mat_vol_kwargs is None:
                 mat_vol_kwargs = {}
+            mat_vol_kwargs.setdefault('bounding_boxes', True)
             self.results['mesh_material_volumes'] = mmv = comm.bcast(
                 self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs))
 
@@ -539,14 +540,15 @@ class R2SManager:
     def get_decay_photon_source_mesh(
         self,
         time_index: int = -1
-    ) -> list[openmc.MeshSource]:
+    ) -> list[openmc.IndependentSource]:
         """Create decay photon source for a mesh-based calculation.
 
-        This function creates N :class:`MeshSource` objects where N is the
-        maximum number of unique materials that appears in a single mesh
-        element. For each mesh element-material combination, and
-        IndependentSource instance is created with a spatial constraint limited
-        the sampled decay photons to the correct region.
+        For each mesh element-material combination, an
+        :class:`~openmc.IndependentSource` is created with a
+        :class:`~openmc.stats.Box` spatial distribution based on the bounding
+        box of the material within the mesh element. A material constraint is
+        also applied so that sampled source sites are limited to the correct
+        region.
 
         When the photon transport model is different from the neutron model, the
         photon MeshMaterialVolumes is used to determine whether an (element,
@@ -559,19 +561,15 @@ class R2SManager:
 
         Returns
         -------
-        list of openmc.MeshSource
-            A list of MeshSource objects, each containing IndependentSource
-            instances for the decay photons in the corresponding mesh element.
+        list of openmc.IndependentSource
+            A list of IndependentSource objects for the decay photons, one for
+            each mesh element-material combination with non-zero source strength.
 
         """
         mat_dict = self.neutron_model._get_all_materials()
 
-        # Some MeshSource objects will have empty positions; create a "null source"
-        # that is used for this case
-        null_source = openmc.IndependentSource(particle='photon', strength=0.0)
-
-        # List to hold sources for each MeshSource (length = N)
-        source_lists = []
+        # List to hold all sources
+        sources = []
 
         # Index in the overall list of activated materials
         index_mat = 0
@@ -594,7 +592,7 @@ class R2SManager:
                     if mat_id is not None
                 }
 
-            for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)):
+            for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
                 # Skip void volume
                 if mat_id is None:
                     continue
@@ -604,30 +602,27 @@ class R2SManager:
                     index_mat += 1
                     continue
 
-                # Check whether a new MeshSource object is needed
-                if j >= len(source_lists):
-                    source_lists.append([null_source]*n_elements)
-
                 # Get activated material composition
                 original_mat = materials[index_mat]
                 activated_mat = results[time_index].get_material(str(original_mat.id))
 
-                # Create decay photon source source
+                # Create decay photon source
                 energy = activated_mat.get_decay_photon_energy()
                 if energy is not None:
                     strength = energy.integral()
-                    source_lists[j][index_elem] = openmc.IndependentSource(
+                    space = openmc.stats.Box(*bbox)
+                    sources.append(openmc.IndependentSource(
+                        space=space,
                         energy=energy,
                         particle='photon',
                         strength=strength,
                         constraints={'domains': [mat_dict[mat_id]]}
-                    )
+                    ))
 
                 # Increment index of activated material
                 index_mat += 1
 
-        # Return list of mesh sources
-        return [openmc.MeshSource(self.domains, sources) for sources in source_lists]
+        return sources
 
     def load_results(self, path: PathLike):
         """Load results from a previous R2S calculation.

openmc/deplete/stepresult.py

@@ -70,6 +70,7 @@ class StepResult:
         self.index_mat = None
         self.index_nuc = None
         self.mat_to_hdf5_ind = None
+        self.name_list = None
 
         self.data = None
 
@@ -138,7 +139,7 @@ class StepResult:
     def n_hdf5_mats(self):
         return len(self.mat_to_hdf5_ind)
 
-    def allocate(self, volume, nuc_list, burn_list, full_burn_list):
+    def allocate(self, volume, nuc_list, burn_list, full_burn_list, name_list=None):
         """Allocate memory for depletion step data
 
         Parameters
@@ -151,12 +152,15 @@ class StepResult:
             A list of all mat IDs to be burned.  Used for sorting the simulation.
         full_burn_list : list of str
             List of all burnable material IDs
+        name_list : list of str, optional
+            Material names corresponding to materials in burn_list
 
         """
         self.volume = copy.deepcopy(volume)
         self.index_nuc = {nuc: i for i, nuc in enumerate(nuc_list)}
         self.index_mat = {mat: i for i, mat in enumerate(burn_list)}
         self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)}
+        self.mat_to_name = dict(zip(burn_list, name_list)) if name_list is not None else {}
 
         # Create storage array
         self.data = np.zeros((self.n_mat, self.n_nuc))
@@ -184,7 +188,7 @@ class StepResult:
 
         # Direct transfer
         direct_attrs = ("time", "k", "source_rate", "index_nuc",
-                        "mat_to_hdf5_ind", "proc_time")
+                        "mat_to_hdf5_ind", "mat_to_name", "proc_time")
         for attr in direct_attrs:
             setattr(new, attr, getattr(self, attr))
         # Get applicable slice of data
@@ -223,6 +227,8 @@ class StepResult:
                 f'mat_id {mat_id} not found in StepResult. Available mat_id '
                 f'values are {list(self.volume.keys())}'
             ) from e
+        if mat_id in self.mat_to_name:
+            material.name = self.mat_to_name[mat_id]
         for nuc, _ in sorted(self.index_nuc.items(), key=lambda x: x[1]):
             atoms = self[mat_id, nuc]
             if atoms <= 0.0:
@@ -313,6 +319,8 @@ class StepResult:
             mat_single_group = mat_group.create_group(mat)
             mat_single_group.attrs["index"] = self.mat_to_hdf5_ind[mat]
             mat_single_group.attrs["volume"] = self.volume[mat]
+            if mat in self.mat_to_name:
+                mat_single_group.attrs["name"] = self.mat_to_name[mat]
 
         nuc_group = handle.create_group("nuclides")
 
@@ -495,6 +503,7 @@ class StepResult:
         results.volume = {}
         results.index_mat = {}
         results.index_nuc = {}
+        results.mat_to_name = {}
         rxn_nuc_to_ind = {}
         rxn_to_ind = {}
 
@@ -504,6 +513,8 @@ class StepResult:
 
             results.volume[mat] = vol
             results.index_mat[mat] = ind
+            if "name" in mat_handle.attrs:
+                results.mat_to_name[mat] = mat_handle.attrs["name"]
 
         for nuc, nuc_handle in handle["/nuclides"].items():
             ind_atom = nuc_handle.attrs["atom number index"]
@@ -569,11 +580,11 @@ class StepResult:
             .. versionadded:: 0.14.0
         """
         # Get indexing terms
-        vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
+        vol_dict, nuc_list, burn_list, full_burn_list, name_list = op.get_results_info()
 
         # Create results
         results = StepResult()
-        results.allocate(vol_dict, nuc_list, burn_list, full_burn_list)
+        results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, name_list)
 
         n_mat = len(burn_list)
 

openmc/element.py

@@ -184,8 +184,6 @@ class Element(str):
                 for nuclide in absent_nuclides:
                     if nuclide in ['O17', 'O18'] and 'O16' in mutual_nuclides:
                         abundances['O16'] += NATURAL_ABUNDANCE[nuclide]
-                    elif nuclide == 'Ta180_m1' and 'Ta180' in library_nuclides:
-                            abundances['Ta180'] = NATURAL_ABUNDANCE[nuclide]
                     elif nuclide == 'Ta180_m1' and 'Ta181' in mutual_nuclides:
                             abundances['Ta181'] += NATURAL_ABUNDANCE[nuclide]
                     elif nuclide == 'W180' and 'W182' in mutual_nuclides:

openmc/examples.py

@@ -4,7 +4,7 @@ import numpy as np
 
 import openmc
 
-
+PINCELL_PITCH = 1.26 # cm
 
 def pwr_pin_cell() -> openmc.Model:
     """Create a PWR pin-cell model.
@@ -51,7 +51,7 @@ def pwr_pin_cell() -> openmc.Model:
     model.materials = (fuel, clad, hot_water)
 
     # Instantiate ZCylinder surfaces
-    pitch = 1.26
+    pitch = PINCELL_PITCH
     fuel_or = openmc.ZCylinder(x0=0, y0=0, r=0.39218, name='Fuel OR')
     clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
     left = openmc.XPlane(x0=-pitch/2, name='left', boundary_type='reflective')
@@ -319,14 +319,14 @@ def pwr_core() -> openmc.Model:
     l100 = openmc.RectLattice(
         name='Fuel assembly (lower half)', lattice_id=100)
     l100.lower_left = (-10.71, -10.71)
-    l100.pitch = (1.26, 1.26)
+    l100.pitch = (PINCELL_PITCH, PINCELL_PITCH)
     l100.universes = np.tile(fuel_cold, (17, 17))
     l100.universes[tube_x, tube_y] = tube_cold
 
     l101 = openmc.RectLattice(
         name='Fuel assembly (upper half)', lattice_id=101)
     l101.lower_left = (-10.71, -10.71)
-    l101.pitch = (1.26, 1.26)
+    l101.pitch = (PINCELL_PITCH, PINCELL_PITCH)
     l101.universes = np.tile(fuel_hot, (17, 17))
     l101.universes[tube_x, tube_y] = tube_hot
 
@@ -350,7 +350,7 @@ def pwr_core() -> openmc.Model:
     # Define core lattices
     l200 = openmc.RectLattice(name='Core lattice (lower half)', lattice_id=200)
     l200.lower_left = (-224.91, -224.91)
-    l200.pitch = (21.42, 21.42)
+    l200.pitch = (17 * PINCELL_PITCH, 17 * PINCELL_PITCH)
     l200.universes = [
         [fa_cw]*21,
         [fa_cw]*21,
@@ -376,7 +376,7 @@ def pwr_core() -> openmc.Model:
 
     l201 = openmc.RectLattice(name='Core lattice (lower half)', lattice_id=201)
     l201.lower_left = (-224.91, -224.91)
-    l201.pitch = (21.42, 21.42)
+    l201.pitch = (17 * PINCELL_PITCH, 17 * PINCELL_PITCH)
     l201.universes = [
         [fa_hw]*21,
         [fa_hw]*21,
@@ -488,7 +488,7 @@ def pwr_assembly() -> openmc.Model:
     clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
 
     # Create boundary planes to surround the geometry
-    pitch = 21.42
+    pitch = 17 * PINCELL_PITCH
     min_x = openmc.XPlane(x0=-pitch/2, boundary_type='reflective')
     max_x = openmc.XPlane(x0=+pitch/2, boundary_type='reflective')
     min_y = openmc.YPlane(y0=-pitch/2, boundary_type='reflective')
@@ -514,7 +514,7 @@ def pwr_assembly() -> openmc.Model:
 
     # Create fuel assembly Lattice
     assembly = openmc.RectLattice(name='Fuel Assembly')
-    assembly.pitch = (pitch/17, pitch/17)
+    assembly.pitch = (PINCELL_PITCH, PINCELL_PITCH)
     assembly.lower_left = (-pitch/2, -pitch/2)
 
     # Create array indices for guide tube locations in lattice
@@ -656,24 +656,21 @@ def slab_mg(num_regions=1, mat_names=None, mgxslib_name='2g.h5') -> openmc.Model
     return model
 
 
-def random_ray_lattice() -> openmc.Model:
-    """Create a 2x2 PWR pincell asymmetrical lattic eexample.
-
-    This model is a 2x2 reflective lattice of fuel pins with one of the lattice
-    locations having just moderator instead of a fuel pin. It uses 7 group
-    cross section data.
+def _generate_c5g7_materials() -> openmc.Materials:
+    """Generate materials utilizing multi-group cross sections based on the
+    the C5G7 Benchmark.
 
     Returns
     -------
-    model : openmc.Model
-        A PWR 2x2 lattice model
+    materials : openmc.Materials
+        Materials object containing UO2 and water materials.
 
+    Data Sources
+    ------------
+    All cross section data are from:
+    Lewis et al., "Benchmark specification for determinisitc 2D/3D MOX fuel
+    assembly transport calculations without spatial homogenization"
     """
-    model = openmc.Model()
-
-    ###########################################################################
-    # Create MGXS data for the problem
-
     # Instantiate the energy group data
     group_edges = [1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6, 1.0e6, 20.0e6]
     groups = openmc.mgxs.EnergyGroups(group_edges)
@@ -704,9 +701,10 @@ def random_ray_lattice() -> openmc.Model:
     uo2_xsdata.set_fission([7.21206e-03, 8.19301e-04, 6.45320e-03,
                             1.85648e-02, 1.78084e-02, 8.30348e-02,
                             2.16004e-01])
-    uo2_xsdata.set_nu_fission([2.005998e-02, 2.027303e-03, 1.570599e-02,
-                               4.518301e-02, 4.334208e-02, 2.020901e-01,
-                               5.257105e-01])
+    nu_fission = np.array([2.005998e-02, 2.027303e-03, 1.570599e-02,
+                           4.518301e-02, 4.334208e-02, 2.020901e-01,
+                           5.257105e-01])
+    uo2_xsdata.set_nu_fission(nu_fission)
     uo2_xsdata.set_chi([5.8791e-01, 4.1176e-01, 3.3906e-04, 1.1761e-07, 0.0000e+00,
                         0.0000e+00, 0.0000e+00])
 
@@ -752,14 +750,28 @@ def random_ray_lattice() -> openmc.Model:
     # Instantiate a Materials collection and export to XML
     materials = openmc.Materials([uo2, water])
     materials.cross_sections = "mgxs.h5"
+    return materials
 
-    ###########################################################################
-    # Define problem geometry
 
+def _generate_subdivided_pin_cell(uo2, water) -> openmc.Universe:
+    """Create a radially and azimuthally subdivided pin cell universe. Helper
+    function for random_ray_pin_cell() and random_ray_lattice()
+
+    Parameters
+    ----------
+    uo2 : openmc.Material
+        UO2 material
+    water : openmc.Material
+        Water material
+
+    Returns
+    -------
+    pincell : openmc.Universe
+        Universe containing an unbounded pin cell
+
+    """
     ########################################
-    # Define an unbounded pincell universe
-
-    pitch = 1.26
+    # Define an unbounded pin cell universe
 
     # Create a surface for the fuel outer radius
     fuel_or = openmc.ZCylinder(r=0.54, name='Fuel OR')
@@ -781,7 +793,7 @@ def random_ray_lattice() -> openmc.Model:
     moderator_c = openmc.Cell(
         fill=water, region=+outer_ring_b, name='moderator outer c')
 
-    # Create pincell universe
+    # Create pin cell universe
     pincell_base = openmc.Universe()
 
     # Register Cells with Universe
@@ -801,19 +813,125 @@ def random_ray_lattice() -> openmc.Model:
     for i in range(8):
         azimuthal_cell = openmc.Cell(name=f'azimuthal_cell_{i}')
         azimuthal_cell.fill = pincell_base
-        azimuthal_cell.region = +azimuthal_planes[i] & -azimuthal_planes[(i+1) % 8]
+        azimuthal_cell.region = + \
+            azimuthal_planes[i] & -azimuthal_planes[(i+1) % 8]
         azimuthal_cells.append(azimuthal_cell)
 
     # Create a geometry with the azimuthal universes
     pincell = openmc.Universe(cells=azimuthal_cells, name='pincell')
 
+    return pincell
+
+
+def random_ray_pin_cell() -> openmc.Model:
+    """Create a PWR pin cell example using C5G7 cross section data.
+    cross section data.
+
+    Returns
+    -------
+    model : openmc.Model
+        A PWR pin cell model
+
+    """
+    model = openmc.Model()
+
+    ###########################################################################
+    # Create Materials for the problem
+    materials = _generate_c5g7_materials()
+    uo2 = materials[0]
+    water = materials[1]
+
+    ###########################################################################
+    # Define problem geometry
+    pincell = _generate_subdivided_pin_cell(uo2, water)
+
+    ########################################
+    # Define cell containing lattice and other stuff
+    pitch = PINCELL_PITCH
+    box = openmc.model.RectangularPrism(pitch, pitch, boundary_type='reflective')
+
+    pincell = openmc.Cell(fill=pincell, region=-box, name='pincell')
+
+    # Create a geometry with the top-level cell
+    geometry = openmc.Geometry([pincell])
+
+    ###########################################################################
+    # Define problem settings
+
+    # Instantiate a Settings object, set all runtime parameters, and export to XML
+    settings = openmc.Settings()
+    settings.energy_mode = "multi-group"
+    settings.batches = 400
+    settings.inactive = 200
+    settings.particles = 100
+
+    # Create an initial uniform spatial source distribution over fissionable zones
+    lower_left = (-pitch / 2, -pitch / 2, -1)
+    upper_right = (pitch / 2, pitch / 2, 1)
+    uniform_dist = openmc.stats.Box(lower_left, upper_right)
+    rr_source = openmc.IndependentSource(space=uniform_dist)
+
+    settings.random_ray['distance_active'] = 100.0
+    settings.random_ray['distance_inactive'] = 20.0
+    settings.random_ray['ray_source'] = rr_source
+    settings.random_ray['volume_normalized_flux_tallies'] = True
+
+    ###########################################################################
+    # Define tallies
+    # Now use the mesh filter in a tally and indicate what scores are desired
+    tally = openmc.Tally(name="Pin tally")
+    tally.scores = ['flux', 'fission', 'nu-fission']
+    tally.estimator = 'analog'
+
+    # Instantiate a Tallies collection and export to XML
+    tallies = openmc.Tallies([tally])
+
+    ###########################################################################
+    #                   Exporting to OpenMC model
+    ###########################################################################
+
+    model.geometry = geometry
+    model.materials = materials
+    model.settings = settings
+    model.tallies = tallies
+    return model
+
+
+def random_ray_lattice() -> openmc.Model:
+    """Create a 2x2 PWR pin cell asymmetrical lattice example.
+
+    This model is a 2x2 reflective lattice of fuel pins with one of the lattice
+    locations having just moderator instead of a fuel pin. It uses C5G7
+    cross section data.
+
+    Returns
+    -------
+    model : openmc.Model
+        A PWR 2x2 lattice model
+
+    """
+    model = openmc.Model()
+
+    ###########################################################################
+    # Create Materials for the problem
+    materials = _generate_c5g7_materials()
+    uo2 = materials[0]
+    water = materials[1]
+
+    ###########################################################################
+    # Define problem geometry
+    pincell = _generate_subdivided_pin_cell(uo2, water)
+
     ########################################
     # Define a moderator lattice universe
 
-    moderator_infinite = openmc.Cell(fill=water, name='moderator infinite')
+    moderator_infinite = openmc.Cell(name='moderator infinite')
+    moderator_infinite.fill = water
+
     mu = openmc.Universe()
     mu.add_cells([moderator_infinite])
 
+    pitch = PINCELL_PITCH
     lattice = openmc.RectLattice()
     lattice.lower_left = [-pitch/2.0, -pitch/2.0]
     lattice.pitch = [pitch/10.0, pitch/10.0]
@@ -837,8 +955,7 @@ def random_ray_lattice() -> openmc.Model:
 
     ########################################
     # Define cell containing lattice and other stuff
-    box = openmc.model.RectangularPrism(
-        pitch*2, pitch*2, boundary_type='reflective')
+    box = openmc.model.RectangularPrism(pitch*2, pitch*2, boundary_type='reflective')
 
     assembly = openmc.Cell(fill=lattice2x2, region=-box, name='assembly')
 

openmc/filter.py

@@ -35,7 +35,6 @@ _CURRENT_NAMES = (
     'z-min out', 'z-min in', 'z-max out', 'z-max in'
 )
 
-_PARTICLES = {'neutron', 'photon', 'electron', 'positron'}
 
 
 class FilterMeta(ABCMeta):
@@ -735,9 +734,8 @@ class ParticleFilter(Filter):
 
     Parameters
     ----------
-    bins : str, or sequence of str
-        The particles to tally represented as strings ('neutron', 'photon',
-        'electron', 'positron').
+    bins : str, int, openmc.ParticleType, or sequence
+        The particle types to tally represented as names, PDG numbers, or types.
     filter_id : int
         Unique identifier for the filter
 
@@ -763,11 +761,16 @@ class ParticleFilter(Filter):
 
     @Filter.bins.setter
     def bins(self, bins):
-        cv.check_type('bins', bins, Sequence, str)
+        if isinstance(bins, (str, Integral, openmc.ParticleType)):
+            bins = [bins]
+        else:
+            cv.check_type('bins', bins, Sequence,
+                          (str, Integral, openmc.ParticleType))
         bins = np.atleast_1d(bins)
-        for edge in bins:
-            cv.check_value('filter bin', edge, _PARTICLES)
-        self._bins = bins
+        normalized = []
+        for entry in bins:
+            normalized.append(str(openmc.ParticleType(entry)))
+        self._bins = np.array(normalized, dtype=str)
 
     @classmethod
     def from_hdf5(cls, group, **kwargs):

openmc/lib/core.py

@@ -28,7 +28,7 @@ class _SourceSite(Structure):
                 ('wgt', c_double),
                 ('delayed_group', c_int),
                 ('surf_id', c_int),
-                ('particle', c_int),
+                ('particle', c_int32),
                 ('parent_nuclide', c_int),
                 ('parent_id', c_int64),
                 ('progeny_id', c_int64)]

openmc/lib/filter.py

@@ -132,6 +132,9 @@ _dll.openmc_meshsurface_filter_set_translation.errcheck = _error_handler
 _dll.openmc_new_filter.argtypes = [c_char_p, POINTER(c_int32)]
 _dll.openmc_new_filter.restype = c_int
 _dll.openmc_new_filter.errcheck = _error_handler
+_dll.openmc_particle_filter_get_bins.argtypes = [c_int32, POINTER(c_int32)]
+_dll.openmc_particle_filter_get_bins.restype = c_int
+_dll.openmc_particle_filter_get_bins.errcheck = _error_handler
 _dll.openmc_spatial_legendre_filter_get_order.argtypes = [c_int32, POINTER(c_int)]
 _dll.openmc_spatial_legendre_filter_get_order.restype = c_int
 _dll.openmc_spatial_legendre_filter_get_order.errcheck = _error_handler
@@ -402,8 +405,8 @@ class MeshFilter(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 
+        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.
 
     """
@@ -454,7 +457,7 @@ class MeshFilter(Filter):
         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()
@@ -598,9 +601,9 @@ class ParticleFilter(Filter):
 
     @property
     def bins(self):
-        particle_i = np.zeros((self.n_bins,), dtype=c_int)
+        particle_i = np.zeros((self.n_bins,), dtype=c_int32)
         _dll.openmc_particle_filter_get_bins(
-            self._index, particle_i.ctypes.data_as(POINTER(c_int)))
+            self._index, particle_i.ctypes.data_as(POINTER(c_int32)))
         return [ParticleType(i) for i in particle_i]
 
 

openmc/lib/mesh.py

@@ -1,5 +1,5 @@
 from collections.abc import Mapping, Sequence
-from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER,
+from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, c_void_p,
                     create_string_buffer, c_size_t)
 from math import sqrt
 import sys
@@ -47,7 +47,8 @@ _dll.openmc_mesh_bounding_box.argtypes = [
 _dll.openmc_mesh_bounding_box.restype = c_int
 _dll.openmc_mesh_bounding_box.errcheck = _error_handler
 _dll.openmc_mesh_material_volumes.argtypes = [
-    c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double]
+    c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double,
+    c_void_p]
 _dll.openmc_mesh_material_volumes.restype = c_int
 _dll.openmc_mesh_material_volumes.errcheck = _error_handler
 _dll.openmc_mesh_get_plot_bins.argtypes = [
@@ -188,6 +189,7 @@ class Mesh(_FortranObjectWithID):
             n_samples: int | tuple[int, int, int] = 10_000,
             max_materials: int = 4,
             output: bool = True,
+            bounding_boxes: bool = False,
     ) -> MeshMaterialVolumes:
         """Determine volume of materials in each mesh element.
 
@@ -213,6 +215,11 @@ class Mesh(_FortranObjectWithID):
             Estimated maximum number of materials in any given mesh element.
         output : bool, optional
             Whether or not to show output.
+        bounding_boxes : bool, optional
+            Whether or not to compute an axis-aligned bounding box for each
+            (mesh element, material) combination. When enabled, the bounding
+            box encloses the ray-estimator prisms used for the volume
+            estimation.
 
         Returns
         -------
@@ -243,23 +250,36 @@ class Mesh(_FortranObjectWithID):
         table_size = slot_factor*max_materials
         materials = np.full((n, table_size), EMPTY_SLOT, dtype=np.int32)
         volumes = np.zeros((n, table_size), dtype=np.float64)
+        bboxes = None
+        if bounding_boxes:
+            bboxes = np.empty((n, table_size, 6), dtype=np.float64)
+            bboxes[..., 0:3] = np.inf
+            bboxes[..., 3:6] = -np.inf
 
         # Run material volume calculation
         while True:
             try:
+                bboxes_ptr = None
+                if bboxes is not None:
+                    bboxes_ptr = bboxes.ctypes.data_as(POINTER(c_double))
                 with quiet_dll(output):
                     _dll.openmc_mesh_material_volumes(
-                        self._index, nx, ny, nz, table_size, materials, volumes)
+                        self._index, nx, ny, nz, table_size, materials,
+                        volumes, bboxes_ptr)
             except AllocationError:
                 # Increase size of result array and try again
                 table_size *= 2
                 materials = np.full((n, table_size), EMPTY_SLOT, dtype=np.int32)
                 volumes = np.zeros((n, table_size), dtype=np.float64)
+                if bounding_boxes:
+                    bboxes = np.empty((n, table_size, 6), dtype=np.float64)
+                    bboxes[..., 0:3] = np.inf
+                    bboxes[..., 3:6] = -np.inf
             else:
                 # If no error, break out of loop
                 break
 
-        return MeshMaterialVolumes(materials, volumes)
+        return MeshMaterialVolumes(materials, volumes, bboxes)
 
     def get_plot_bins(
             self,

openmc/lib/weight_windows.py

@@ -53,11 +53,11 @@ _dll.openmc_weight_windows_get_energy_bounds.argtypes = [c_int32, POINTER(POINTE
 _dll.openmc_weight_windows_get_energy_bounds.restype = c_int
 _dll.openmc_weight_windows_get_energy_bounds.errcheck = _error_handler
 
-_dll.openmc_weight_windows_set_particle.argtypes = [c_int32, c_int]
+_dll.openmc_weight_windows_set_particle.argtypes = [c_int32, c_int32]
 _dll.openmc_weight_windows_set_particle.restype = c_int
 _dll.openmc_weight_windows_set_particle.errcheck = _error_handler
 
-_dll.openmc_weight_windows_get_particle.argtypes = [c_int32, POINTER(c_int)]
+_dll.openmc_weight_windows_get_particle.argtypes = [c_int32, POINTER(c_int32)]
 _dll.openmc_weight_windows_get_particle.restype = c_int
 _dll.openmc_weight_windows_get_particle.errcheck = _error_handler
 
@@ -201,16 +201,13 @@ class WeightWindows(_FortranObjectWithID):
 
     @property
     def particle(self):
-        val = c_int()
+        val = c_int32()
         _dll.openmc_weight_windows_get_particle(self._index, val)
         return ParticleType(val.value)
 
     @particle.setter
     def particle(self, p):
-        if isinstance(p, str):
-            p = ParticleType.from_string(p)
-        else:
-            p = ParticleType(p)
+        p = ParticleType(p)
         _dll.openmc_weight_windows_set_particle(self._index, int(p))
 
     @property
@@ -304,10 +301,10 @@ class WeightWindows(_FortranObjectWithID):
         ----------
         tally : openmc.lib.Tally
             The tally used to create the WeightWindows instance.
-        particle : openmc.ParticleType or str, optional
+        particle : openmc.ParticleType or str or int, optional
             The particle type to use for the WeightWindows instance. Should be
-            specified as an instance of ParticleType or as a string with a value of
-            'neutron' or 'photon'.
+            specified as an instance of ParticleType, a PDG number, or as a
+            name.
 
         Returns
         -------
@@ -317,7 +314,8 @@ class WeightWindows(_FortranObjectWithID):
         Raises
         ------
         ValueError
-            If the particle parameter is not an instance of ParticleType or a string.
+            If the particle parameter is not an instance of ParticleType, a string,
+            or an integer PDG number.
         ValueError
             If the particle parameter is not a valid particle type (i.e., not 'neutron'
             or 'photon').
@@ -328,12 +326,13 @@ class WeightWindows(_FortranObjectWithID):
             If the tally does not have a MeshFilter.
         """
         # do some checks on particle value
-        if not isinstance(particle, (ParticleType, str)):
-            raise ValueError(f"Parameter 'particle' must be {ParticleType} or one of ('neutron', 'photon').")
+        if not isinstance(particle, (ParticleType, str, int)):
+            raise ValueError(
+                f"Parameter 'particle' must be {ParticleType} or one of ('neutron', 'photon')."
+            )
 
         # convert particle type if needed
-        if isinstance(particle, str):
-            particle = ParticleType.from_string(particle)
+        particle = ParticleType(particle)
 
         if particle not in (ParticleType.NEUTRON, ParticleType.PHOTON):
             raise ValueError('Weight windows can only be applied for neutrons or photons')

openmc/mesh.py

@@ -17,6 +17,7 @@ import openmc
 import openmc.checkvalue as cv
 from openmc.checkvalue import PathLike
 from openmc.utility_funcs import change_directory
+from .bounding_box import BoundingBox
 from ._xml import get_elem_list, get_text
 from .mixin import IDManagerMixin
 from .surface import _BOUNDARY_TYPES
@@ -40,6 +41,11 @@ class MeshMaterialVolumes(Mapping):
         Array of shape (elements, max_materials) storing material IDs
     volumes : numpy.ndarray
         Array of shape (elements, max_materials) storing material volumes
+    bboxes : numpy.ndarray, optional
+        Array of shape (elements, max_materials, 6) storing axis-aligned
+        bounding boxes for each (element, material) combination with ordering
+        (xmin, ymin, zmin, xmax, ymax, zmax). Bounding boxes enclose the
+        ray-estimator prisms used to compute volumes.
 
     See Also
     --------
@@ -64,9 +70,30 @@ class MeshMaterialVolumes(Mapping):
     [(2, 31.87963824195591), (1, 6.129949130817542)]
 
     """
-    def __init__(self, materials: np.ndarray, volumes: np.ndarray):
+    def __init__(
+        self,
+        materials: np.ndarray,
+        volumes: np.ndarray,
+        bboxes: np.ndarray | None = None
+    ):
         self._materials = materials
         self._volumes = volumes
+        self._bboxes = bboxes
+
+        if self._bboxes is not None:
+            if self._bboxes.shape[:2] != self._materials.shape:
+                raise ValueError(
+                    'bboxes must have shape (elements, max_materials, 6) '
+                    'matching materials/volumes.'
+                )
+            if self._bboxes.shape[2] != 6:
+                raise ValueError(
+                    'bboxes must have shape (elements, max_materials, 6).'
+                )
+
+    @property
+    def has_bounding_boxes(self) -> bool:
+        return self._bboxes is not None
 
     @property
     def num_elements(self) -> int:
@@ -92,7 +119,11 @@ class MeshMaterialVolumes(Mapping):
             volumes[indices] = self._volumes[indices, i]
         return volumes
 
-    def by_element(self, index_elem: int) -> list[tuple[int | None, float]]:
+    def by_element(
+        self,
+        index_elem: int,
+        include_bboxes: bool = False
+    ) -> list[tuple[int | None, float] | tuple[int | None, float, BoundingBox | None]]:
         """Get a list of volumes for each material within a specific element.
 
         Parameters
@@ -102,15 +133,32 @@ class MeshMaterialVolumes(Mapping):
 
         Returns
         -------
-        list of tuple of (material ID, volume)
+        list of tuple
+            If ``include_bboxes`` is False (default), returns tuples of
+            (material ID, volume). If ``include_bboxes`` is True, returns
+            tuples of (material ID, volume, bounding box).
 
         """
         table_size = self._volumes.shape[1]
-        return [
-            (m if m > -1 else None, self._volumes[index_elem, i])
-            for i in range(table_size)
-            if (m := self._materials[index_elem, i]) != -2
-        ]
+        if include_bboxes and self._bboxes is None:
+            raise ValueError('Bounding boxes were not computed for this object.')
+
+        results = []
+        for i in range(table_size):
+            m = self._materials[index_elem, i]
+            if m == -2:
+                continue
+            mat_id = m if m > -1 else None
+            vol = self._volumes[index_elem, i]
+
+            if include_bboxes:
+                vals = self._bboxes[index_elem, i]
+                bbox = BoundingBox(vals[0:3], vals[3:6])
+                results.append((mat_id, vol, bbox))
+            else:
+                results.append((mat_id, vol))
+
+        return results
 
     def save(self, filename: PathLike):
         """Save material volumes to a .npz file.
@@ -120,8 +168,10 @@ class MeshMaterialVolumes(Mapping):
         filename : path-like
             Filename where data will be saved
         """
-        np.savez_compressed(
-            filename, materials=self._materials, volumes=self._volumes)
+        kwargs = {'materials': self._materials, 'volumes': self._volumes}
+        if self._bboxes is not None:
+            kwargs['bboxes'] = self._bboxes
+        np.savez_compressed(filename, **kwargs)
 
     @classmethod
     def from_npz(cls, filename: PathLike) -> MeshMaterialVolumes:
@@ -134,7 +184,8 @@ class MeshMaterialVolumes(Mapping):
 
         """
         filedata = np.load(filename)
-        return cls(filedata['materials'], filedata['volumes'])
+        bboxes = filedata['bboxes'] if 'bboxes' in filedata.files else None
+        return cls(filedata['materials'], filedata['volumes'], bboxes)
 
 
 class MeshBase(IDManagerMixin, ABC):
@@ -153,11 +204,17 @@ class MeshBase(IDManagerMixin, ABC):
         Unique identifier for the mesh
     name : str
         Name of the mesh
+    lower_left : Iterable of float
+        The lower-left coordinates
+    upper_right : Iterable of float
+        The upper-right coordinates
     bounding_box : openmc.BoundingBox
         Axis-aligned bounding box of the mesh as defined by the upper-right and
         lower-left coordinates.
     indices : Iterable of tuple
         An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, 1), ...]
+    n_elements : int
+        Number of elements in the mesh
     """
 
     next_id = 1
@@ -179,6 +236,16 @@ class MeshBase(IDManagerMixin, ABC):
             self._name = name
         else:
             self._name = ''
+            
+    @property
+    @abstractmethod
+    def lower_left(self):
+        pass
+        
+    @property
+    @abstractmethod
+    def upper_right(self):
+        pass
 
     @property
     def bounding_box(self) -> openmc.BoundingBox:
@@ -188,6 +255,11 @@ class MeshBase(IDManagerMixin, ABC):
     @abstractmethod
     def indices(self):
         pass
+        
+    @property
+    @abstractmethod
+    def n_elements(self):
+        pass
 
     def __repr__(self):
         string = type(self).__name__ + '\n'
@@ -366,6 +438,7 @@ class MeshBase(IDManagerMixin, ABC):
             model: openmc.Model,
             n_samples: int | tuple[int, int, int] = 10_000,
             max_materials: int = 4,
+            bounding_boxes: bool = False,
             **kwargs
     ) -> MeshMaterialVolumes:
         """Determine volume of materials in each mesh element.
@@ -388,6 +461,11 @@ class MeshBase(IDManagerMixin, ABC):
             the x, y, and z dimensions.
         max_materials : int, optional
             Estimated maximum number of materials in any given mesh element.
+        bounding_boxes : bool, optional
+            Whether to compute an axis-aligned bounding box for each
+            (mesh element, material) combination. When enabled, the bounding
+            box encloses the ray-estimator prisms used for the volume
+            estimation.
         **kwargs : dict
             Keyword arguments passed to :func:`openmc.lib.init`
 
@@ -419,7 +497,8 @@ class MeshBase(IDManagerMixin, ABC):
 
             # Compute material volumes
             volumes = mesh.material_volumes(
-                n_samples, max_materials, output=kwargs['output'])
+                n_samples, max_materials, output=kwargs['output'],
+                bounding_boxes=bounding_boxes)
 
         # Restore original tallies
         model.tallies = original_tallies
@@ -557,10 +636,19 @@ class StructuredMesh(MeshBase):
         s0 = (slice(0, -1),)*ndim + (slice(None),)
         s1 = (slice(1, None),)*ndim + (slice(None),)
         return (vertices[s0] + vertices[s1]) / 2
+    
+    @property
+    def n_elements(self):
+        return np.prod(self.dimension)
 
     @property
     def num_mesh_cells(self):
-        return np.prod(self.dimension)
+        warnings.warn(
+            "The 'num_mesh_cells' attribute is deprecated and will be removed in a future version. "
+            "Use 'n_elements' instead.",
+            FutureWarning, stacklevel=2
+        )
+        return self.n_elements
 
     def write_data_to_vtk(self,
                           filename: PathLike,
@@ -822,10 +910,10 @@ class StructuredMesh(MeshBase):
         """
         cv.check_type('data label', label, str)
 
-        if dataset.size != self.num_mesh_cells:
+        if dataset.size != self.n_elements:
             raise ValueError(
                 f"The size of the dataset '{label}' ({dataset.size}) should be"
-                f" equal to the number of mesh cells ({self.num_mesh_cells})"
+                f" equal to the number of mesh cells ({self.n_elements})"
             )
 
         # accept a flat array as-is, assuming it is in the correct order

openmc/model/model.py

@@ -1160,6 +1160,16 @@ 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 and check that mg cross sections path is accessible
+        for mat in self.geometry.get_all_materials().values():
+            if mat._macroscopic is not None:
+                self.settings.energy_mode = 'multi-group'
+                if 'mg_cross_sections' not in openmc.config:
+                    raise RuntimeError("'mg_cross_sections' path must be set in "
+                                       "openmc.config before plotting.")
+                break
+
         # Get ID map from the C API
         id_map = self.id_map(
             origin=origin,
@@ -1179,8 +1189,8 @@ class Model:
 
         # Convert ID map to RGB image
         img = id_map_to_rgb(
-            id_map=id_map, 
-            color_by=color_by, 
+            id_map=id_map,
+            color_by=color_by,
             colors=colors,
             overlap_color=overlap_color
         )
@@ -1217,7 +1227,7 @@ class Model:
                 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)
@@ -1685,6 +1695,87 @@ class Model:
                 self.geometry.get_all_materials().values()
             )
 
+    def _auto_generate_mgxs_lib(
+        self,
+        model: openmc.model.model,
+        groups: openmc.mgxs.EnergyGroups,
+        correction: str | none,
+        directory: pathlike,
+    ) -> openmc.mgxs.Library:
+        """
+        Automatically generate a multi-group cross section libray from a model
+        with the specified group structure.
+
+        Parameters
+        ----------
+        groups : openmc.mgxs.EnergyGroups
+            Energy group structure for the MGXS.
+        nparticles : int
+            Number of particles to simulate per batch when generating MGXS.
+        mgxs_path : str
+            Filename for the MGXS HDF5 file.
+        correction : str
+            Transport correction to apply to the MGXS. Options are None and
+            "P0".
+        directory : str
+            Directory to run the simulation in, so as to contain XML files.
+
+        Returns
+        -------
+        mgxs_lib : openmc.mgxs.Library
+            OpenMC MGXS Library object
+        """
+
+        # Initialize MGXS library with a finished OpenMC geometry object
+        mgxs_lib = openmc.mgxs.Library(model.geometry)
+
+        # Pick energy group structure
+        mgxs_lib.energy_groups = groups
+
+        # Disable transport correction
+        mgxs_lib.correction = correction
+
+        # Specify needed cross sections for random ray
+        if correction == 'P0':
+            mgxs_lib.mgxs_types = [
+                'nu-transport', 'absorption', 'nu-fission', 'fission',
+                'consistent nu-scatter matrix', 'multiplicity matrix', 'chi',
+                'kappa-fission'
+            ]
+        elif correction is None:
+            mgxs_lib.mgxs_types = [
+                'total', 'absorption', 'nu-fission', 'fission',
+                'consistent nu-scatter matrix', 'multiplicity matrix', 'chi',
+                'kappa-fission'
+            ]
+
+        # Specify a "material" domain type for the cross section tally filters
+        mgxs_lib.domain_type = "material"
+
+        # Specify the domains over which to compute multi-group cross sections
+        mgxs_lib.domains = model.geometry.get_all_materials().values()
+
+        # Do not compute cross sections on a nuclide-by-nuclide basis
+        mgxs_lib.by_nuclide = False
+
+        # Check the library - if no errors are raised, then the library is satisfactory.
+        mgxs_lib.check_library_for_openmc_mgxs()
+
+        # Construct all tallies needed for the multi-group cross section library
+        mgxs_lib.build_library()
+
+        # Create a "tallies.xml" file for the MGXS Library
+        mgxs_lib.add_to_tallies(model.tallies, merge=True)
+
+        # Run
+        statepoint_filename = model.run(cwd=directory)
+
+        # Load MGXS
+        with openmc.StatePoint(statepoint_filename) as sp:
+            mgxs_lib.load_from_statepoint(sp)
+
+        return mgxs_lib
+
     def _create_mgxs_sources(
         self,
         groups: openmc.mgxs.EnergyGroups,
@@ -1848,52 +1939,8 @@ class Model:
             model.geometry.root_universe = infinite_universe
 
             # Add MGXS Tallies
-
-            # Initialize MGXS library with a finished OpenMC geometry object
-            mgxs_lib = openmc.mgxs.Library(model.geometry)
-
-            # Pick energy group structure
-            mgxs_lib.energy_groups = groups
-
-            # Disable transport correction
-            mgxs_lib.correction = correction
-
-            # Specify needed cross sections for random ray
-            if correction == 'P0':
-                mgxs_lib.mgxs_types = [
-                    'nu-transport', 'absorption', 'nu-fission', 'fission',
-                    'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'
-                ]
-            elif correction is None:
-                mgxs_lib.mgxs_types = [
-                    'total', 'absorption', 'nu-fission', 'fission',
-                    'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'
-                ]
-
-            # Specify a "cell" domain type for the cross section tally filters
-            mgxs_lib.domain_type = "material"
-
-            # Specify the cell domains over which to compute multi-group cross sections
-            mgxs_lib.domains = model.geometry.get_all_materials().values()
-
-            # Do not compute cross sections on a nuclide-by-nuclide basis
-            mgxs_lib.by_nuclide = False
-
-            # Check the library - if no errors are raised, then the library is satisfactory.
-            mgxs_lib.check_library_for_openmc_mgxs()
-
-            # Construct all tallies needed for the multi-group cross section library
-            mgxs_lib.build_library()
-
-            # Create a "tallies.xml" file for the MGXS Library
-            mgxs_lib.add_to_tallies(model.tallies, merge=True)
-
-            # Run
-            statepoint_filename = model.run(cwd=directory)
-
-            # Load MGXS
-            with openmc.StatePoint(statepoint_filename) as sp:
-                mgxs_lib.load_from_statepoint(sp)
+            mgxs_lib = self._auto_generate_mgxs_lib(
+                model, groups, correction, directory)
 
             # Create a MGXS File which can then be written to disk
             mgxs_set = mgxs_lib.get_xsdata(domain=material, xsdata_name=name)
@@ -2055,48 +2102,8 @@ class Model:
         model.settings.output = {'summary': True, 'tallies': False}
 
         # Add MGXS Tallies
-
-        # Initialize MGXS library with a finished OpenMC geometry object
-        mgxs_lib = openmc.mgxs.Library(model.geometry)
-
-        # Pick energy group structure
-        mgxs_lib.energy_groups = groups
-
-        # Disable transport correction
-        mgxs_lib.correction = correction
-
-       # Specify needed cross sections for random ray
-        if correction == 'P0':
-            mgxs_lib.mgxs_types = ['nu-transport', 'absorption', 'nu-fission', 'fission',
-                                   'consistent nu-scatter matrix', 'multiplicity matrix', 'chi']
-        elif correction is None:
-            mgxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission',
-                                   'consistent nu-scatter matrix', 'multiplicity matrix', 'chi']
-
-        # Specify a "cell" domain type for the cross section tally filters
-        mgxs_lib.domain_type = "material"
-
-        # Specify the cell domains over which to compute multi-group cross sections
-        mgxs_lib.domains = model.geometry.get_all_materials().values()
-
-        # Do not compute cross sections on a nuclide-by-nuclide basis
-        mgxs_lib.by_nuclide = False
-
-        # Check the library - if no errors are raised, then the library is satisfactory.
-        mgxs_lib.check_library_for_openmc_mgxs()
-
-        # Construct all tallies needed for the multi-group cross section library
-        mgxs_lib.build_library()
-
-        # Create a "tallies.xml" file for the MGXS Library
-        mgxs_lib.add_to_tallies(model.tallies, merge=True)
-
-        # Run
-        statepoint_filename = model.run(cwd=directory)
-
-        # Load MGXS
-        with openmc.StatePoint(statepoint_filename) as sp:
-            mgxs_lib.load_from_statepoint(sp)
+        mgxs_lib = self._auto_generate_mgxs_lib(
+                model, groups, correction, directory)
 
         names = [mat.name for mat in mgxs_lib.domains]
 
@@ -2146,52 +2153,8 @@ class Model:
         model.settings.output = {'summary': True, 'tallies': False}
 
         # Add MGXS Tallies
-
-        # Initialize MGXS library with a finished OpenMC geometry object
-        mgxs_lib = openmc.mgxs.Library(model.geometry)
-
-        # Pick energy group structure
-        mgxs_lib.energy_groups = groups
-
-        # Disable transport correction
-        mgxs_lib.correction = correction
-
-        # Specify needed cross sections for random ray
-        if correction == 'P0':
-            mgxs_lib.mgxs_types = [
-                'nu-transport', 'absorption', 'nu-fission', 'fission',
-                'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'
-            ]
-        elif correction is None:
-            mgxs_lib.mgxs_types = [
-                'total', 'absorption', 'nu-fission', 'fission',
-                'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'
-            ]
-
-        # Specify a "cell" domain type for the cross section tally filters
-        mgxs_lib.domain_type = "material"
-
-        # Specify the cell domains over which to compute multi-group cross sections
-        mgxs_lib.domains = model.geometry.get_all_materials().values()
-
-        # Do not compute cross sections on a nuclide-by-nuclide basis
-        mgxs_lib.by_nuclide = False
-
-        # Check the library - if no errors are raised, then the library is satisfactory.
-        mgxs_lib.check_library_for_openmc_mgxs()
-
-        # Construct all tallies needed for the multi-group cross section library
-        mgxs_lib.build_library()
-
-        # Create a "tallies.xml" file for the MGXS Library
-        mgxs_lib.add_to_tallies(model.tallies, merge=True)
-
-        # Run
-        statepoint_filename = model.run(cwd=directory)
-
-        # Load MGXS
-        with openmc.StatePoint(statepoint_filename) as sp:
-            mgxs_lib.load_from_statepoint(sp)
+        mgxs_lib = self._auto_generate_mgxs_lib(
+                model, groups, correction, directory)
 
         names = [mat.name for mat in mgxs_lib.domains]
 
@@ -2617,5 +2580,3 @@ class SearchResult:
     def total_batches(self) -> int:
         """Total number of active batches used across all evaluations."""
         return sum(self.batches)
-
-

openmc/particle_restart.py

@@ -1,6 +1,7 @@
 import h5py
 
 import openmc.checkvalue as cv
+from .particle_type import ParticleType
 
 _VERSION_PARTICLE_RESTART = 2
 
@@ -28,8 +29,8 @@ class Particle:
         Type of simulation (criticality or fixed source)
     id : long
         Identifier of the particle
-    type : int
-        Particle type (1 = neutron, 2 = photon, 3 = electron, 4 = positron)
+    type : openmc.ParticleType
+        Particle type
     weight : float
         Weight of the particle
     energy : float
@@ -52,7 +53,7 @@ class Particle:
             self.energy = f['energy'][()]
             self.generations_per_batch = f['generations_per_batch'][()]
             self.id = f['id'][()]
-            self.type = f['type'][()]
+            self.type = ParticleType(f['type'][()])
             self.n_particles = f['n_particles'][()]
             self.run_mode = f['run_mode'][()].decode()
             self.uvw = f['uvw'][()]

openmc/particle_type.py

@@ -0,0 +1,228 @@
+from numbers import Integral
+
+from openmc.data import gnds_name, zam, ATOMIC_SYMBOL
+
+
+_PDG_NAME = {
+    2112: 'neutron',
+    22: 'photon',
+    11: 'electron',
+    -11: 'positron',
+    2212: 'H1',
+}
+
+_ALIAS_PDG = {
+    'neutron': 2112,
+    'n': 2112,
+    'photon': 22,
+    'gamma': 22,
+    'electron': 11,
+    'positron': -11,
+    'proton': 2212,
+    'p': 2212,
+    'h1': 2212,
+    'deuteron': 1000010020,
+    'd': 1000010020,
+    'h2': 1000010020,
+    'triton': 1000010030,
+    't': 1000010030,
+    'h3': 1000010030,
+    'alpha': 1000020040,
+    'he4': 1000020040,
+}
+
+_LEGACY_PARTICLE_INDEX = {
+    0: 2112,
+    1: 22,
+    2: 11,
+    3: -11,
+}
+
+
+class ParticleType:
+    """Particle type defined by a PDG number.
+
+    ParticleType uses the Particle Data Group (PDG) Monte Carlo numbering scheme
+    to uniquely identify particle types. This includes elementary particles
+    (neutrons, photons, etc.) and nuclear codes for isotopes.
+
+    Parameters
+    ----------
+    value : str, int, or ParticleType
+        The particle identifier. Can be:
+
+        - A string name (e.g., 'neutron', 'photon', 'He4', 'U235')
+        - An integer PDG number (e.g., 2112 for neutron)
+        - A string with PDG prefix (e.g., 'pdg:2112')
+        - An existing ParticleType instance
+
+    Attributes
+    ----------
+    pdg_number : int
+        The PDG number for this particle type
+    zam : tuple of int or None
+        For nuclear particles, the (Z, A, m) tuple where Z is atomic number,
+        A is mass number, and m is metastable state. None for elementary particles.
+    is_nucleus : bool
+        Whether this particle is a nucleus (ion)
+
+    Examples
+    --------
+    >>> neutron = ParticleType('neutron')
+    >>> neutron.pdg_number
+    2112
+    >>> he4 = ParticleType('He4')
+    >>> he4.zam
+    (2, 4, 0)
+    >>> ParticleType(2112) == ParticleType('neutron')
+    True
+
+    """
+
+    __slots__ = ('_pdg_number',)
+
+    def __init__(self, value: 'str | int | ParticleType'):
+        if isinstance(value, ParticleType):
+            pdg = value._pdg_number
+        elif isinstance(value, str):
+            pdg = self._pdg_number_from_string(value)
+        elif isinstance(value, Integral):
+            pdg = int(value)
+            # Handle legacy particle indices (0, 1, 2, 3)
+            if pdg in _LEGACY_PARTICLE_INDEX:
+                pdg = _LEGACY_PARTICLE_INDEX[pdg]
+        else:
+            raise TypeError(f"Cannot create ParticleType from {type(value).__name__}")
+
+        self._pdg_number = pdg
+
+    def __eq__(self, other):
+        if isinstance(other, ParticleType):
+            return self._pdg_number == other._pdg_number
+        if isinstance(other, Integral):
+            return self._pdg_number == int(other)
+        if isinstance(other, str):
+            try:
+                return self._pdg_number == ParticleType(other)._pdg_number
+            except (ValueError, TypeError):
+                return False
+        return NotImplemented
+
+    def __hash__(self) -> int:
+        return hash(self._pdg_number)
+
+    def __int__(self) -> int:
+        return self._pdg_number
+
+    @property
+    def pdg_number(self) -> int:
+        return self._pdg_number
+
+    @staticmethod
+    def _pdg_number_from_string(value: str) -> int:
+        """Parse a string to get a PDG number.
+
+        Parameters
+        ----------
+        value : str
+            Particle identifier string
+
+        Returns
+        -------
+        int
+            PDG number
+
+        Raises
+        ------
+        ValueError
+            If string cannot be parsed as a valid particle identifier
+
+        """
+        s = value.strip()
+        if not s:
+            raise ValueError('Particle identifier cannot be empty.')
+
+        lower = s.lower()
+        if lower.startswith('pdg:'):
+            code_str = lower[4:]
+            try:
+                return int(code_str)
+            except ValueError:
+                raise ValueError(f'Invalid PDG number: {code_str}')
+
+        if lower in _ALIAS_PDG:
+            return _ALIAS_PDG[lower]
+
+        # Assume it is a GNDS nuclide name
+        Z, A, m = zam(s)
+        if Z <= 0 or Z > 999 or A <= 0 or A > 999 or m < 0 or m > 9:
+            raise ValueError('Invalid Z/A/m for nuclear PDG number.')
+        return 1000000000 + Z * 10000 + A * 10 + m
+
+    def __repr__(self) -> str:
+        return f'<ParticleType: {str(self)} (PDG={self._pdg_number})>'
+
+    def __str__(self) -> str:
+        """Return a canonical string representation of the particle type.
+
+        Returns
+        -------
+        str
+            Canonical name (e.g., 'neutron', 'He4', 'pdg:12345')
+
+        """
+        if self._pdg_number in _PDG_NAME:
+            return _PDG_NAME[self._pdg_number]
+
+        if (zam_tuple := self.zam) is not None:
+            Z, A, m = zam_tuple
+            if Z <= 0 or Z > max(ATOMIC_SYMBOL) or A <= 0 or A > 999:
+                raise ValueError(f"Invalid nuclear PDG number: {self._pdg_number}")
+            return gnds_name(Z, A, m)
+
+        return f'pdg:{self._pdg_number}'
+
+    @property
+    def zam(self) -> 'tuple[int, int, int] | None':
+        """Return the (Z, A, m) tuple for nuclear particles.
+
+        Returns
+        -------
+        tuple of int or None
+            For nuclear particles, returns (Z, A, m) where Z is atomic number,
+            A is mass number, and m is metastable state. Returns None for
+            elementary particles.
+
+        """
+        if self._pdg_number < 1000000000:
+            return None
+        Z = (self._pdg_number // 10000) % 1000
+        A = (self._pdg_number // 10) % 1000
+        m = self._pdg_number % 10
+        if Z <= 0 or A <= 0:
+            return None
+        else:
+            return (Z, A, m)
+
+    @property
+    def is_nucleus(self) -> bool:
+        """Return whether this particle is a nucleus.
+
+        Returns
+        -------
+        bool
+            True if the particle is a nucleus (ion), False otherwise
+
+        """
+        return self.zam is not None
+
+
+# Define common particle constants
+ParticleType.NEUTRON = ParticleType(2112)
+ParticleType.PHOTON = ParticleType(22)
+ParticleType.ELECTRON = ParticleType(11)
+ParticleType.POSITRON = ParticleType(-11)
+ParticleType.PROTON = ParticleType(2212)
+ParticleType.DEUTERON = ParticleType(1000010020)
+ParticleType.TRITON = ParticleType(1000010030)
+ParticleType.ALPHA = ParticleType(1000020040)

openmc/settings.py

@@ -4,13 +4,14 @@ import itertools
 from math import ceil
 from numbers import Integral, Real
 from pathlib import Path
+import traceback
 
 import lxml.etree as ET
 import warnings
 import openmc
 import openmc.checkvalue as cv
 from openmc.checkvalue import PathLike
-from openmc.stats.multivariate import MeshSpatial, Box, PolarAzimuthal, Isotropic
+from openmc.stats.multivariate import MeshSpatial
 from ._xml import clean_indentation, get_elem_list, get_text
 from .mesh import _read_meshes, RegularMesh, MeshBase
 from .source import SourceBase, MeshSource, IndependentSource
@@ -467,6 +468,16 @@ class Settings:
         for key, value in kwargs.items():
             setattr(self, key, value)
 
+    def __setattr__(self, name: str, value):
+        if not name.startswith('_'):
+            try:
+                getattr(self, name)
+            except AttributeError as e:
+                msg, = traceback.format_exception_only(e)
+                msg = msg.strip().split(maxsplit=1)[-1]
+                warnings.warn(msg, stacklevel=2)
+        super().__setattr__(name, value)
+
     @property
     def run_mode(self) -> str:
         return self._run_mode.value

openmc/source.py

@@ -1,7 +1,6 @@
 from __future__ import annotations
 from abc import ABC, abstractmethod
 from collections.abc import Iterable, Sequence
-from enum import IntEnum
 from numbers import Real
 from pathlib import Path
 import warnings
@@ -19,6 +18,8 @@ from openmc.stats.multivariate import UnitSphere, Spatial
 from openmc.stats.univariate import Univariate
 from ._xml import get_elem_list, get_text
 from .mesh import MeshBase, StructuredMesh, UnstructuredMesh
+from .particle_type import ParticleType
+from .statepoint import _VERSION_STATEPOINT
 from .utility_funcs import input_path
 
 
@@ -265,8 +266,8 @@ class IndependentSource(SourceBase):
         time distribution of source sites
     strength : float
         Strength of the source
-    particle : {'neutron', 'photon', 'electron', 'positron'}
-        Source particle type
+    particle : str or int or openmc.ParticleType
+        Source particle type (name, PDG number, or type)
     domains : iterable of openmc.Cell, openmc.Material, or openmc.Universe
         Domains to reject based on, i.e., if a sampled spatial location is not
         within one of these domains, it will be rejected.
@@ -302,10 +303,9 @@ class IndependentSource(SourceBase):
     type : str
         Indicator of source type: 'independent'
 
-    .. versionadded:: 0.14.0
-
-    particle : {'neutron', 'photon', 'electron', 'positron'}
-        Source particle type
+        .. versionadded:: 0.14.0
+    particle : str or int or openmc.ParticleType
+        Source particle type (alias, PDG number, or GNDS nuclide name)
     constraints : dict
         Constraints on sampled source particles. Valid keys include
         'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds',
@@ -320,7 +320,7 @@ class IndependentSource(SourceBase):
         energy: openmc.stats.Univariate | None = None,
         time: openmc.stats.Univariate | None = None,
         strength: float = 1.0,
-        particle: str = 'neutron',
+        particle: str | int | ParticleType = 'neutron',
         domains: Sequence[openmc.Cell | openmc.Material |
                           openmc.Universe] | None = None,
         constraints: dict[str, Any] | None = None
@@ -405,14 +405,12 @@ class IndependentSource(SourceBase):
         self._time = time
 
     @property
-    def particle(self):
+    def particle(self) -> ParticleType:
         return self._particle
 
     @particle.setter
     def particle(self, particle):
-        cv.check_value('source particle', particle,
-                       ['neutron', 'photon', 'electron', 'positron'])
-        self._particle = particle
+        self._particle = ParticleType(particle)
 
     def populate_xml_element(self, element):
         """Add necessary source information to an XML element
@@ -423,7 +421,7 @@ class IndependentSource(SourceBase):
             XML element containing source data
 
         """
-        element.set("particle", self.particle)
+        element.set("particle", str(self.particle))
         if self.space is not None:
             element.append(self.space.to_xml_element())
         if self.angle is not None:
@@ -572,10 +570,10 @@ class MeshSource(SourceBase):
         s = np.asarray(s)
 
         if isinstance(self.mesh, StructuredMesh):
-            if s.size != self.mesh.num_mesh_cells:
+            if s.size != self.mesh.n_elements:
                 raise ValueError(
                     f'The length of the source array ({s.size}) does not match '
-                    f'the number of mesh elements ({self.mesh.num_mesh_cells}).')
+                    f'the number of mesh elements ({self.mesh.n_elements}).')
 
             # If user gave a multidimensional array, flatten in the order
             # of the mesh indices
@@ -898,76 +896,6 @@ class FileSource(SourceBase):
         return cls(**kwargs)
 
 
-class ParticleType(IntEnum):
-    """
-    IntEnum class representing a particle type. Type
-    values mirror those found in the C++ class.
-    """
-    NEUTRON = 0
-    PHOTON = 1
-    ELECTRON = 2
-    POSITRON = 3
-
-    @classmethod
-    def from_string(cls, value: str):
-        """
-        Constructs a ParticleType instance from a string.
-
-        Parameters
-        ----------
-        value : str
-            The string representation of the particle type.
-
-        Returns
-        -------
-        The corresponding ParticleType instance.
-        """
-        try:
-            return cls[value.upper()]
-        except KeyError:
-            raise ValueError(
-                f"Invalid string for creation of {cls.__name__}: {value}")
-
-    @classmethod
-    def from_pdg_number(cls, pdg_number: int) -> ParticleType:
-        """Constructs a ParticleType instance from a PDG number.
-
-        The Particle Data Group at LBNL publishes a Monte Carlo particle
-        numbering scheme as part of the `Review of Particle Physics
-        <10.1103/PhysRevD.110.030001>`_. This method maps PDG numbers to the
-        corresponding :class:`ParticleType`.
-
-        Parameters
-        ----------
-        pdg_number : int
-            The PDG number of the particle type.
-
-        Returns
-        -------
-        The corresponding ParticleType instance.
-        """
-        try:
-            return {
-                2112: ParticleType.NEUTRON,
-                22: ParticleType.PHOTON,
-                11: ParticleType.ELECTRON,
-                -11: ParticleType.POSITRON,
-            }[pdg_number]
-        except KeyError:
-            raise ValueError(f"Unrecognized PDG number: {pdg_number}")
-
-    def __repr__(self) -> str:
-        """
-        Returns a string representation of the ParticleType instance.
-
-        Returns:
-            str: The lowercase name of the ParticleType instance.
-        """
-        return self.name.lower()
-
-    # needed for < Python 3.11
-    def __str__(self) -> str:
-        return self.__repr__()
 
 
 class SourceParticle:
@@ -992,8 +920,8 @@ class SourceParticle:
         Delayed group particle was created in (neutrons only)
     surf_id : int
         Surface ID where particle is at, if any.
-    particle : ParticleType
-        Type of the particle
+    particle : ParticleType or str or int
+        Type of the particle (type, name, or PDG number)
 
     """
 
@@ -1006,7 +934,7 @@ class SourceParticle:
         wgt: float = 1.0,
         delayed_group: int = 0,
         surf_id: int = 0,
-        particle: ParticleType = ParticleType.NEUTRON
+        particle: ParticleType | str | int = ParticleType.NEUTRON
     ):
 
         self.r = tuple(r)
@@ -1018,9 +946,16 @@ class SourceParticle:
         self.surf_id = surf_id
         self.particle = particle
 
+    @property
+    def particle(self) -> ParticleType:
+        return self._particle
+
+    @particle.setter
+    def particle(self, particle):
+        self._particle = ParticleType(particle)
+
     def __repr__(self):
-        name = self.particle.name.lower()
-        return f'<SourceParticle: {name} at E={self.E:.6e} eV>'
+        return f'<SourceParticle: {str(self.particle)} at E={self.E:.6e} eV>'
 
     def to_tuple(self) -> tuple:
         """Return source particle attributes as a tuple
@@ -1032,7 +967,7 @@ class SourceParticle:
 
         """
         return (self.r, self.u, self.E, self.time, self.wgt,
-                self.delayed_group, self.surf_id, self.particle.value)
+                self.delayed_group, self.surf_id, self.particle.pdg_number)
 
 
 def write_source_file(
@@ -1116,12 +1051,7 @@ class ParticleList(list):
         particles = []
         with mcpl.MCPLFile(filename) as f:
             for particle in f.particles:
-                # Determine particle type based on the PDG number
-                try:
-                    particle_type = ParticleType.from_pdg_number(
-                        particle.pdgcode)
-                except ValueError:
-                    particle_type = "UNKNOWN"
+                particle_type = ParticleType(particle.pdgcode)
 
                 # Create a source particle instance. Note that MCPL stores
                 # energy in MeV and time in ms.
@@ -1179,7 +1109,7 @@ class ParticleList(list):
         # Extract the attributes of the source particles into a list of tuples
         data = [(sp.r[0], sp.r[1], sp.r[2], sp.u[0], sp.u[1], sp.u[2],
                  sp.E, sp.time, sp.wgt, sp.delayed_group, sp.surf_id,
-                 sp.particle.name.lower()) for sp in self]
+                 str(sp.particle)) for sp in self]
 
         # Define the column names for the DataFrame
         columns = ['x', 'y', 'z', 'u_x', 'u_y', 'u_z', 'E', 'time', 'wgt',
@@ -1226,6 +1156,7 @@ class ParticleList(list):
         kwargs.setdefault('mode', 'w')
         with h5py.File(filename, **kwargs) as fh:
             fh.attrs['filetype'] = np.bytes_("source")
+            fh.attrs['version'] = np.array([_VERSION_STATEPOINT, 2])
             fh.create_dataset('source_bank', data=arr, dtype=source_dtype)
 
 
@@ -1337,7 +1268,7 @@ def read_collision_track_mcpl(file_path):
             data['material_id'].append(int(values_dict.get('material_id', 0)))
             data['universe_id'].append(int(values_dict.get('universe_id', 0)))
             data['n_collision'].append(int(values_dict.get('n_collision', 0)))
-            data['particle'].append(ParticleType.from_pdg_number(p.pdgcode))
+            data['particle'].append(ParticleType(p.pdgcode))
             data['parent_id'].append(int(values_dict.get('parent_id', 0)))
             data['progeny_id'].append(int(values_dict.get('progeny_id', 0)))
 

openmc/stats/univariate.py

@@ -397,7 +397,7 @@ class Discrete(Univariate):
     def merge(
         cls,
         dists: Sequence[Discrete],
-        probs: Sequence[int]
+        probs: Sequence[float]
     ):
         """Merge multiple discrete distributions into a single distribution
 
@@ -1897,7 +1897,7 @@ class Mixture(Univariate):
 
 
 def combine_distributions(
-    dists: Sequence[Univariate],
+    dists: Sequence[Discrete | Tabular],
     probs: Sequence[float]
 ):
     """Combine distributions with specified probabilities
@@ -1912,41 +1912,40 @@ def combine_distributions(
 
     Parameters
     ----------
-    dists : iterable of openmc.stats.Univariate
+    dists : sequence of openmc.stats.Discrete or openmc.stats.Tabular
         Distributions to combine
-    probs : iterable of float
+    probs : sequence of float
         Probability (or intensity) of each distribution
 
     """
-    # Get copy of distribution list so as not to modify the argument
-    dist_list = deepcopy(dists)
+    for i, dist in enumerate(dists):
+        cv.check_type(f'dists[{i}]', dist, (Discrete, Tabular))
+        cv.check_type(f'probs[{i}]', probs[i], Real)
+        cv.check_greater_than(f'probs[{i}]', probs[i], 0.0)
 
     # Get list of discrete/continuous distribution indices
-    discrete_index = [i for i, d in enumerate(dist_list) if isinstance(d, Discrete)]
-    cont_index = [i for i, d in enumerate(dist_list) if isinstance(d, Tabular)]
+    discrete_index = [i for i, d in enumerate(dists) if isinstance(d, Discrete)]
+    cont_index = [i for i, d in enumerate(dists) if isinstance(d, Tabular)]
 
-    # Apply probabilites to continuous distributions
-    for i in cont_index:
-        dist = dist_list[i]
-        dist._p *= probs[i]
+    cont_dists = [dists[i] for i in cont_index]
+    cont_probs = [probs[i] for i in cont_index]
 
     if discrete_index:
         # Create combined discrete distribution
-        dist_discrete = [dist_list[i] for i in discrete_index]
+        dist_discrete = [dists[i] for i in discrete_index]
         discrete_probs = [probs[i] for i in discrete_index]
         combined_dist = Discrete.merge(dist_discrete, discrete_probs)
-
-        # Replace multiple discrete distributions with merged
-        for idx in reversed(discrete_index):
-            dist_list.pop(idx)
-        dist_list.append(combined_dist)
-
-    # Combine discrete and continuous if present
-    if len(dist_list) > 1:
-        probs = [1.0]*len(dist_list)
-        dist_list[:] = [Mixture(probs, dist_list.copy())]
-
-    return dist_list[0]
+        if cont_index:
+            return Mixture(cont_probs + [1.0], cont_dists + [combined_dist])
+        else:
+            return combined_dist
+    else:
+        if len(cont_dists) == 1:
+            dist = cont_dists[0]
+            return Tabular(dist.x, dist.p * cont_probs[0],
+                           dist.interpolation, bias=dist.bias)
+        else:
+            return Mixture(cont_probs, cont_dists)
 
 
 def check_bias_support(parent: Univariate, bias: Univariate | None):

openmc/tallies.py

@@ -1982,7 +1982,7 @@ class Tally(IDManagerMixin):
 
         # Expand the columns into Pandas MultiIndices for readability
         if pd.__version__ >= '0.16':
-            columns = copy.deepcopy(df.columns.values)
+            columns = copy.deepcopy(list(df.columns.values))
 
             # Convert all elements in columns list to tuples
             for i, column in enumerate(columns):

openmc/tracks.py

@@ -4,7 +4,8 @@ from collections.abc import Sequence
 import h5py
 
 from .checkvalue import check_filetype_version
-from .source import SourceParticle, ParticleType
+from .particle_type import ParticleType
+from .source import SourceParticle
 
 from pathlib import Path
 
@@ -25,7 +26,7 @@ states : numpy.ndarray
 
 """
 def _particle_track_repr(self):
-    return f"<ParticleTrack: {self.particle}, {len(self.states)} states>"
+    return f"<ParticleTrack: {str(self.particle)}, {len(self.states)} states>"
 ParticleTrack.__repr__ = _particle_track_repr
 
 
@@ -92,8 +93,8 @@ class Track(Sequence):
 
         Parameters
         ----------
-        particle : {'neutron', 'photon', 'electron', 'positron'}
-            Matching particle type
+        particle : str or int or openmc.ParticleType
+            Matching particle type (name, PDG number, or type)
         state_filter : function
             Function that takes a state (structured datatype) and returns a bool
             depending on some criteria.
@@ -126,7 +127,7 @@ class Track(Sequence):
         for t in self:
             # Check for matching particle
             if particle is not None:
-                if t.particle.name.lower() != particle:
+                if t.particle != ParticleType(particle):
                     continue
 
             # Apply arbitrary state filter
@@ -184,7 +185,7 @@ class Track(Sequence):
     def sources(self):
         sources = []
         for particle_track in self:
-            particle_type = ParticleType(particle_track.particle)
+            particle_type = particle_track.particle
             state = particle_track.states[0]
             sources.append(
                 SourceParticle(

openmc/weight_windows.py

@@ -10,13 +10,12 @@ import numpy as np
 import h5py
 
 import openmc
-from openmc.filter import _PARTICLES
 from openmc.mesh import MeshBase, RectilinearMesh, CylindricalMesh, SphericalMesh, UnstructuredMesh
 import openmc.checkvalue as cv
 from openmc.checkvalue import PathLike
 from ._xml import get_elem_list, get_text, clean_indentation
 from .mixin import IDManagerMixin
-from .utility_funcs import change_directory
+from .particle_type import ParticleType
 
 
 class WeightWindows(IDManagerMixin):
@@ -51,7 +50,7 @@ class WeightWindows(IDManagerMixin):
         A list of values for which each successive pair constitutes a range of
         energies in [eV] for a single bin. If no energy bins are provided, the
         maximum and minimum energy for the data available at runtime.
-    particle_type : {'neutron', 'photon'}
+    particle_type : str or int or openmc.ParticleType
         Particle type the weight windows apply to
     survival_ratio : float
         Ratio of the survival weight to the lower weight window bound for
@@ -116,7 +115,7 @@ class WeightWindows(IDManagerMixin):
         upper_ww_bounds: Iterable[float] | None = None,
         upper_bound_ratio: float | None = None,
         energy_bounds: Iterable[Real] | None = None,
-        particle_type: str = 'neutron',
+        particle_type: str | int | openmc.ParticleType = 'neutron',
         survival_ratio: float = 3.0,
         max_lower_bound_ratio: float | None = None,
         max_split: int = 10,
@@ -213,13 +212,15 @@ class WeightWindows(IDManagerMixin):
         self._mesh = mesh
 
     @property
-    def particle_type(self) -> str:
+    def particle_type(self) -> ParticleType:
         return self._particle_type
 
     @particle_type.setter
-    def particle_type(self, pt: str):
-        cv.check_value('Particle type', pt, _PARTICLES)
-        self._particle_type = pt
+    def particle_type(self, pt):
+        ptype = ParticleType(pt)
+        if ptype not in {ParticleType.NEUTRON, ParticleType.PHOTON}:
+            raise ValueError("Weight windows can only be applied for neutrons or photons")
+        self._particle_type = ptype
 
     @property
     def energy_bounds(self) -> Iterable[Real]:
@@ -329,7 +330,7 @@ class WeightWindows(IDManagerMixin):
         subelement.text = str(self.mesh.id)
 
         subelement = ET.SubElement(element, 'particle_type')
-        subelement.text = self.particle_type
+        subelement.text = str(self.particle_type)
 
         if self.energy_bounds is not None:
             subelement = ET.SubElement(element, 'energy_bounds')
@@ -494,7 +495,7 @@ class WeightWindowGenerator:
         A list of values for which each successive pair constitutes a range of
         energies in [eV] for a single bin. If no energy bins are provided, the
         maximum and minimum energy for the data available at runtime.
-    particle_type : {'neutron', 'photon'}
+    particle_type : str or int or openmc.ParticleType
         Particle type the weight windows apply to
     method : {'magic', 'fw_cadis'}
         The weight window generation methodology applied during an update.
@@ -513,7 +514,7 @@ class WeightWindowGenerator:
     energy_bounds : Iterable of Real
         A list of values for which each successive pair constitutes a range of
         energies in [eV] for a single bin
-    particle_type : {'neutron', 'photon'}
+    particle_type : openmc.ParticleType
         Particle type the weight windows apply to
     method : {'magic', 'fw_cadis'}
         The weight window generation methodology applied during an update.
@@ -534,7 +535,7 @@ class WeightWindowGenerator:
         self,
         mesh: openmc.MeshBase,
         energy_bounds: Sequence[float] | None = None,
-        particle_type: str = 'neutron',
+        particle_type: str | int | openmc.ParticleType = 'neutron',
         method: str = 'magic',
         max_realizations: int = 1,
         update_interval: int = 1,
@@ -555,7 +556,7 @@ class WeightWindowGenerator:
     def __repr__(self):
         string = type(self).__name__ + '\n'
         string += f'\t{"Mesh":<20}=\t{self.mesh.id}\n'
-        string += f'\t{"Particle:":<20}=\t{self.particle_type}\n'
+        string += f'\t{"Particle:":<20}=\t{str(self.particle_type)}\n'
         string += f'\t{"Energy Bounds:":<20}=\t{self.energy_bounds}\n'
         string += f'\t{"Method":<20}=\t{self.method}\n'
         string += f'\t{"Max Realizations:":<20}=\t{self.max_realizations}\n'
@@ -586,13 +587,15 @@ class WeightWindowGenerator:
         self._energy_bounds = eb
 
     @property
-    def particle_type(self) -> str:
+    def particle_type(self) -> ParticleType:
         return self._particle_type
 
     @particle_type.setter
-    def particle_type(self, pt: str):
-        cv.check_value('particle type', pt, ('neutron', 'photon'))
-        self._particle_type = pt
+    def particle_type(self, pt):
+        ptype = ParticleType(pt)
+        if ptype not in {ParticleType.NEUTRON, ParticleType.PHOTON}:
+            raise ValueError("Weight windows can only be applied for neutrons or photons")
+        self._particle_type = ptype
 
     @property
     def method(self) -> str:
@@ -695,7 +698,7 @@ class WeightWindowGenerator:
             subelement = ET.SubElement(element, 'energy_bounds')
             subelement.text = ' '.join(str(e) for e in self.energy_bounds)
         particle_elem = ET.SubElement(element, 'particle_type')
-        particle_elem.text = self.particle_type
+        particle_elem.text = str(self.particle_type)
         realizations_elem = ET.SubElement(element, 'max_realizations')
         realizations_elem.text = str(self.max_realizations)
         update_interval_elem = ET.SubElement(element, 'update_interval')
@@ -730,7 +733,7 @@ class WeightWindowGenerator:
 
         mesh_id = int(get_text(elem, 'mesh'))
         mesh = meshes[mesh_id]
-        
+
         energy_bounds = get_elem_list(elem, "energy_bounds, float")
         particle_type = get_text(elem, 'particle_type')
 

src/bremsstrahlung.cpp

@@ -31,13 +31,13 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
   if (p.material() == MATERIAL_VOID)
     return;
 
-  int photon = static_cast<int>(ParticleType::photon);
+  int photon = ParticleType::photon().transport_index();
   if (p.E() < settings::energy_cutoff[photon])
     return;
 
   // Get bremsstrahlung data for this material and particle type
   BremsstrahlungData* mat;
-  if (p.type() == ParticleType::positron) {
+  if (p.type() == ParticleType::positron()) {
     mat = &model::materials[p.material()]->ttb_->positron;
   } else {
     mat = &model::materials[p.material()]->ttb_->electron;
@@ -119,7 +119,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
       }
 
       // Create secondary photon
-      p.create_secondary(p.wgt(), p.u(), w, ParticleType::photon);
+      p.create_secondary(p.wgt(), p.u(), w, ParticleType::photon());
       *E_lost += w;
     }
   }

src/cell.cpp

@@ -660,7 +660,7 @@ Region::Region(std::string region_spec, int32_t cell_id)
       if (token == OP_UNION) {
         simple_ = false;
         // Ensure intersections have precedence over unions
-        add_precedence();
+        enforce_precedence();
         break;
       }
     }
@@ -703,7 +703,7 @@ void Region::apply_demorgan(
 //! precedence than unions using parentheses.
 //==============================================================================
 
-int64_t Region::add_parentheses(int64_t start)
+void Region::add_parentheses(int64_t start)
 {
   int32_t start_token = expression_[start];
   // Add left parenthesis and set new position to be after parenthesis
@@ -712,14 +712,6 @@ int64_t Region::add_parentheses(int64_t start)
   }
   expression_.insert(expression_.begin() + start - 1, OP_LEFT_PAREN);
 
-  // Keep track of return iterator distance. If we don't encounter a left
-  // parenthesis, we return an iterator corresponding to wherever the right
-  // parenthesis is inserted. If a left parenthesis is encountered, an iterator
-  // corresponding to the left parenthesis is returned. Also note that we keep
-  // track of a *distance* instead of an iterator because the underlying memory
-  // allocation may change.
-  std::size_t return_it_dist = 0;
-
   // Add right parenthesis
   // While the start iterator is within the bounds of infix
   while (start + 1 < expression_.size()) {
@@ -733,7 +725,6 @@ int64_t Region::add_parentheses(int64_t start)
       // in the region, when the operator is an intersection then include the
       // operator and next surface
       if (expression_[start] == OP_LEFT_PAREN) {
-        return_it_dist = start;
         int depth = 1;
         do {
           start++;
@@ -750,54 +741,73 @@ int64_t Region::add_parentheses(int64_t start)
           --start;
         }
         expression_.insert(expression_.begin() + start, OP_RIGHT_PAREN);
-        if (return_it_dist > 0) {
-          return return_it_dist;
-        } else {
-          return start - 1;
-        }
+        return;
       }
     }
   }
-  // If we get here a right parenthesis hasn't been placed,
-  // return iterator
+  // If we get here a right parenthesis hasn't been placed
   expression_.push_back(OP_RIGHT_PAREN);
-  if (return_it_dist > 0) {
-    return return_it_dist;
-  } else {
-    return start - 1;
-  }
 }
 
+//==============================================================================
+//! Add parentheses to enforce operator precedence in region expressions
+//!
+//! This function ensures that intersection operators have higher precedence
+//! than union operators by adding parentheses where needed. For example:
+//!   "1 2 | 3" becomes "(1 2) | 3"
+//!   "1 | 2 3" becomes "1 | (2 3)"
+//!
+//! The algorithm uses stacks to track the current operator type and its
+//! position at each parenthesis depth level. When it encounters a different
+//! operator at the same depth, it adds parentheses to group the
+//! higher-precedence operations.
 //==============================================================================
 
-void Region::add_precedence()
+void Region::enforce_precedence()
 {
-  int32_t current_op = 0;
-  std::size_t current_dist = 0;
+  // Stack tracking the operator type at each depth (0 = no operator seen yet)
+  vector<int32_t> op_stack = {0};
 
-  for (int64_t i = 0; i < expression_.size(); i++) {
+  // Stack tracking where the operator sequence started at each depth
+  vector<std::size_t> pos_stack = {0};
+
+  for (int64_t i = 0; i < expression_.size(); ++i) {
     int32_t token = expression_[i];
 
-    if (token == OP_UNION || token == OP_INTERSECTION) {
-      if (current_op == 0) {
-        // Set the current operator if is hasn't been set
-        current_op = token;
-        current_dist = i;
-      } else if (token != current_op) {
-        // If the current operator doesn't match the token, add parenthesis to
-        // assert precedence
-        if (current_op == OP_INTERSECTION) {
-          i = add_parentheses(current_dist);
-        } else {
-          i = add_parentheses(i);
-        }
-        current_op = 0;
-        current_dist = 0;
+    if (token == OP_LEFT_PAREN) {
+      // Entering a new parenthesis level - push new tracking state
+      op_stack.push_back(0);
+      pos_stack.push_back(0);
+      continue;
+    } else if (token == OP_RIGHT_PAREN) {
+      // Exiting a parenthesis level - pop tracking state (keep at least one)
+      if (op_stack.size() > 1) {
+        op_stack.pop_back();
+        pos_stack.pop_back();
+      }
+      continue;
+    }
+
+    if (token == OP_UNION || token == OP_INTERSECTION) {
+      if (op_stack.back() == 0) {
+        // First operator at this depth - record it and its position
+        op_stack.back() = token;
+        pos_stack.back() = i;
+      } else if (token != op_stack.back()) {
+        // Encountered a different operator at the same depth - need to add
+        // parentheses to enforce precedence. Intersection has higher
+        // precedence, so we parenthesize the intersection terms.
+        if (op_stack.back() == OP_INTERSECTION) {
+          add_parentheses(pos_stack.back());
+        } else {
+          add_parentheses(i);
+        }
+
+        // Restart the scan since we modified the expression
+        i = -1; // Will be incremented to 0 by the for loop
+        op_stack = {0};
+        pos_stack = {0};
       }
-    } else if (token > OP_COMPLEMENT) {
-      // If the token is a parenthesis reset the current operator
-      current_op = 0;
-      current_dist = 0;
     }
   }
 }

src/distribution_spatial.cpp

@@ -426,6 +426,11 @@ SpatialBox::SpatialBox(pugi::xml_node node, bool fission)
   upper_right_ = Position {params[3], params[4], params[5]};
 }
 
+SpatialBox::SpatialBox(Position lower_left, Position upper_right, bool fission)
+  : lower_left_(lower_left), upper_right_(upper_right),
+    only_fissionable_(fission)
+{}
+
 std::pair<Position, double> SpatialBox::sample(uint64_t* seed) const
 {
   Position xi {prn(seed), prn(seed), prn(seed)};

src/material.cpp

@@ -819,9 +819,9 @@ void Material::calculate_xs(Particle& p) const
   p.macro_xs().fission = 0.0;
   p.macro_xs().nu_fission = 0.0;
 
-  if (p.type() == ParticleType::neutron) {
+  if (p.type().is_neutron()) {
     this->calculate_neutron_xs(p);
-  } else if (p.type() == ParticleType::photon) {
+  } else if (p.type().is_photon()) {
     this->calculate_photon_xs(p);
   }
 }
@@ -829,7 +829,7 @@ void Material::calculate_xs(Particle& p) const
 void Material::calculate_neutron_xs(Particle& p) const
 {
   // Find energy index on energy grid
-  int neutron = static_cast<int>(ParticleType::neutron);
+  int neutron = ParticleType::neutron().transport_index();
   int i_grid =
     std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing;
 

src/mcpl_interface.cpp

@@ -321,25 +321,7 @@ inline void ensure_mcpl_ready_or_fatal()
 SourceSite mcpl_particle_to_site(const mcpl_particle_repr_t* particle_repr)
 {
   SourceSite site;
-  switch (particle_repr->pdgcode) {
-  case 2112:
-    site.particle = ParticleType::neutron;
-    break;
-  case 22:
-    site.particle = ParticleType::photon;
-    break;
-  case 11:
-    site.particle = ParticleType::electron;
-    break;
-  case -11:
-    site.particle = ParticleType::positron;
-    break;
-  default:
-    fatal_error(fmt::format(
-      "MCPL: Encountered unexpected PDG code {} when converting to SourceSite.",
-      particle_repr->pdgcode));
-    break;
-  }
+  site.particle = ParticleType {particle_repr->pdgcode};
 
   // Copy position and direction
   site.r.x = particle_repr->position[0];
@@ -368,7 +350,6 @@ vector<SourceSite> mcpl_source_sites(std::string path)
   }
 
   size_t n_particles_in_file = g_mcpl_api->hdr_nparticles(mcpl_file);
-  size_t n_skipped = 0;
   if (n_particles_in_file > 0) {
     sites.reserve(n_particles_in_file);
   }
@@ -381,31 +362,16 @@ vector<SourceSite> mcpl_source_sites(std::string path)
         path, sites.size(), n_particles_in_file));
       break;
     }
-    if (p_repr->pdgcode == 2112 || p_repr->pdgcode == 22 ||
-        p_repr->pdgcode == 11 || p_repr->pdgcode == -11) {
-      sites.push_back(mcpl_particle_to_site(p_repr));
-    } else {
-      n_skipped++;
-    }
+    sites.push_back(mcpl_particle_to_site(p_repr));
   }
 
   g_mcpl_api->close_file(mcpl_file);
 
-  if (n_skipped > 0 && n_particles_in_file > 0) {
-    double percent_skipped =
-      100.0 * static_cast<double>(n_skipped) / n_particles_in_file;
-    warning(fmt::format(
-      "MCPL: Skipped {} of {} total particles ({:.1f}%) in file '{}' because "
-      "their type is not supported by OpenMC.",
-      n_skipped, n_particles_in_file, percent_skipped, path));
-  }
-
   if (sites.empty()) {
     if (n_particles_in_file > 0) {
       fatal_error(fmt::format(
-        "MCPL file '{}' contained {} particles, but none were of the supported "
-        "types (neutron, photon, electron, positron). OpenMC cannot proceed "
-        "without source particles.",
+        "MCPL file '{}' contained {} particles, but no particles could be "
+        "read.",
         path, n_particles_in_file));
     } else {
       fatal_error(fmt::format(
@@ -461,22 +427,7 @@ void write_mcpl_source_bank_internal(mcpl_outfile_t* file_id,
         p_repr.ekin = site.E * 1e-6;
         p_repr.time = site.time * 1e3;
         p_repr.weight = site.wgt;
-        switch (site.particle) {
-        case ParticleType::neutron:
-          p_repr.pdgcode = 2112;
-          break;
-        case ParticleType::photon:
-          p_repr.pdgcode = 22;
-          break;
-        case ParticleType::electron:
-          p_repr.pdgcode = 11;
-          break;
-        case ParticleType::positron:
-          p_repr.pdgcode = -11;
-          break;
-        default:
-          continue;
-        }
+        p_repr.pdgcode = site.particle.pdg_number();
         g_mcpl_api->add_particle(file_id, &p_repr);
       }
     }
@@ -633,22 +584,7 @@ void write_mcpl_collision_track_internal(mcpl_outfile_t* file_id,
       p_repr.ekin = site.E * 1e-6;
       p_repr.time = site.time * 1e3;
       p_repr.weight = site.wgt;
-      switch (site.particle) {
-      case ParticleType::neutron:
-        p_repr.pdgcode = 2112;
-        break;
-      case ParticleType::photon:
-        p_repr.pdgcode = 22;
-        break;
-      case ParticleType::electron:
-        p_repr.pdgcode = 11;
-        break;
-      case ParticleType::positron:
-        p_repr.pdgcode = -11;
-        break;
-      default:
-        continue;
-      }
+      p_repr.pdgcode = site.particle.pdg_number();
       g_mcpl_api->add_particle(file_id, &p_repr);
     }
   } else {

src/mesh.cpp

@@ -1,6 +1,8 @@
 #include "openmc/mesh.h"
 #include <algorithm> // for copy, equal, min, min_element
 #include <cassert>
+#include <cstdint>        // for uint64_t
+#include <cstring>        // for memcpy
 #define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers
 #include <cmath>          // for ceil
 #include <cstddef>        // for size_t
@@ -140,6 +142,63 @@ inline bool atomic_cas_int32(int32_t* ptr, int32_t& expected, int32_t desired)
 #endif
 }
 
+// Helper function equivalent to std::bit_cast in C++20
+template<typename To, typename From>
+inline To bit_cast_value(const From& value)
+{
+  To out;
+  std::memcpy(&out, &value, sizeof(To));
+  return out;
+}
+
+inline void atomic_update_double(double* ptr, double value, bool is_min)
+{
+#if defined(__GNUC__) || defined(__clang__)
+  using may_alias_uint64_t [[gnu::may_alias]] = uint64_t;
+  auto* bits_ptr = reinterpret_cast<may_alias_uint64_t*>(ptr);
+  uint64_t current_bits = __atomic_load_n(bits_ptr, __ATOMIC_SEQ_CST);
+  double current = bit_cast_value<double>(current_bits);
+  while (is_min ? (value < current) : (value > current)) {
+    uint64_t desired_bits = bit_cast_value<uint64_t>(value);
+    uint64_t expected_bits = current_bits;
+    if (__atomic_compare_exchange_n(bits_ptr, &expected_bits, desired_bits,
+          false, __ATOMIC_SEQ_CST, __ATOMIC_SEQ_CST)) {
+      return;
+    }
+    current_bits = expected_bits;
+    current = bit_cast_value<double>(current_bits);
+  }
+
+#elif defined(_MSC_VER)
+  auto* bits_ptr = reinterpret_cast<volatile long long*>(ptr);
+  long long current_bits = *bits_ptr;
+  double current = bit_cast_value<double>(current_bits);
+  while (is_min ? (value < current) : (value > current)) {
+    long long desired_bits = bit_cast_value<long long>(value);
+    long long old_bits =
+      _InterlockedCompareExchange64(bits_ptr, desired_bits, current_bits);
+    if (old_bits == current_bits) {
+      return;
+    }
+    current_bits = old_bits;
+    current = bit_cast_value<double>(current_bits);
+  }
+
+#else
+#error "No compare-and-swap implementation available for this compiler."
+#endif
+}
+
+inline void atomic_max_double(double* ptr, double value)
+{
+  atomic_update_double(ptr, value, false);
+}
+
+inline void atomic_min_double(double* ptr, double value)
+{
+  atomic_update_double(ptr, value, true);
+}
+
 namespace detail {
 
 //==============================================================================
@@ -147,7 +206,7 @@ namespace detail {
 //==============================================================================
 
 void MaterialVolumes::add_volume(
-  int index_elem, int index_material, double volume)
+  int index_elem, int index_material, double volume, const BoundingBox* bbox)
 {
   // This method handles adding elements to the materials hash table,
   // implementing open addressing with linear probing. Consistency across
@@ -166,10 +225,18 @@ void MaterialVolumes::add_volume(
     // Non-atomic read of current material
     int32_t current_val = *slot_ptr;
 
-    // Found the desired material; accumulate volume
+    // Found the desired material; accumulate volume and bbox
     if (current_val == index_material) {
 #pragma omp atomic
       this->volumes(index_elem, slot) += volume;
+      if (bbox) {
+        atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox->min.x);
+        atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox->min.y);
+        atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox->min.z);
+        atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox->max.x);
+        atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox->max.y);
+        atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox->max.z);
+      }
       return;
     }
 
@@ -185,6 +252,14 @@ void MaterialVolumes::add_volume(
       if (claimed_slot || (expected_val == index_material)) {
 #pragma omp atomic
         this->volumes(index_elem, slot) += volume;
+        if (bbox) {
+          atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox->min.x);
+          atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox->min.y);
+          atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox->min.z);
+          atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox->max.x);
+          atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox->max.y);
+          atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox->max.z);
+        }
         return;
       }
     }
@@ -195,7 +270,7 @@ void MaterialVolumes::add_volume(
 }
 
 void MaterialVolumes::add_volume_unsafe(
-  int index_elem, int index_material, double volume)
+  int index_elem, int index_material, double volume, const BoundingBox* bbox)
 {
   // Linear probe
   for (int attempt = 0; attempt < table_size_; ++attempt) {
@@ -207,9 +282,23 @@ void MaterialVolumes::add_volume_unsafe(
     // Read current material
     int32_t current_val = this->materials(index_elem, slot);
 
-    // Found the desired material; accumulate volume
+    // Found the desired material; accumulate volume and bbox
     if (current_val == index_material) {
       this->volumes(index_elem, slot) += volume;
+      if (bbox) {
+        this->bboxes(index_elem, slot, 0) =
+          std::min(this->bboxes(index_elem, slot, 0), bbox->min.x);
+        this->bboxes(index_elem, slot, 1) =
+          std::min(this->bboxes(index_elem, slot, 1), bbox->min.y);
+        this->bboxes(index_elem, slot, 2) =
+          std::min(this->bboxes(index_elem, slot, 2), bbox->min.z);
+        this->bboxes(index_elem, slot, 3) =
+          std::max(this->bboxes(index_elem, slot, 3), bbox->max.x);
+        this->bboxes(index_elem, slot, 4) =
+          std::max(this->bboxes(index_elem, slot, 4), bbox->max.y);
+        this->bboxes(index_elem, slot, 5) =
+          std::max(this->bboxes(index_elem, slot, 5), bbox->max.z);
+      }
       return;
     }
 
@@ -217,6 +306,20 @@ void MaterialVolumes::add_volume_unsafe(
     if (current_val == EMPTY) {
       this->materials(index_elem, slot) = index_material;
       this->volumes(index_elem, slot) += volume;
+      if (bbox) {
+        this->bboxes(index_elem, slot, 0) =
+          std::min(this->bboxes(index_elem, slot, 0), bbox->min.x);
+        this->bboxes(index_elem, slot, 1) =
+          std::min(this->bboxes(index_elem, slot, 1), bbox->min.y);
+        this->bboxes(index_elem, slot, 2) =
+          std::min(this->bboxes(index_elem, slot, 2), bbox->min.z);
+        this->bboxes(index_elem, slot, 3) =
+          std::max(this->bboxes(index_elem, slot, 3), bbox->max.x);
+        this->bboxes(index_elem, slot, 4) =
+          std::max(this->bboxes(index_elem, slot, 4), bbox->max.y);
+        this->bboxes(index_elem, slot, 5) =
+          std::max(this->bboxes(index_elem, slot, 5), bbox->max.z);
+      }
       return;
     }
   }
@@ -333,6 +436,12 @@ vector<double> Mesh::volumes() const
 
 void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
   int32_t* materials, double* volumes) const
+{
+  this->material_volumes(nx, ny, nz, table_size, materials, volumes, nullptr);
+}
+
+void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
+  int32_t* materials, double* volumes, double* bboxes) const
 {
   if (mpi::master) {
     header("MESH MATERIAL VOLUMES CALCULATION", 7);
@@ -351,7 +460,8 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
   timer.start();
 
   // Create object for keeping track of materials/volumes
-  detail::MaterialVolumes result(materials, volumes, table_size);
+  detail::MaterialVolumes result(materials, volumes, bboxes, table_size);
+  bool compute_bboxes = bboxes != nullptr;
 
   // Determine bounding box
   auto bbox = this->bounding_box();
@@ -370,13 +480,13 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
 #pragma omp parallel
   {
     // Preallocate vector for mesh indices and length fractions and particle
-    std::vector<int> bins;
-    std::vector<double> length_fractions;
+    vector<int> bins;
+    vector<double> length_fractions;
     Particle p;
 
     SourceSite site;
     site.E = 1.0;
-    site.particle = ParticleType::neutron;
+    site.particle = ParticleType::neutron();
 
     for (int axis = 0; axis < 3; ++axis) {
       // Set starting position and direction
@@ -453,12 +563,36 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
             if (i_material != C_NONE) {
               i_material = model::materials[i_material]->id();
             }
+            double cumulative_frac = 0.0;
             for (int i_bin = 0; i_bin < bins.size(); i_bin++) {
               int mesh_index = bins[i_bin];
               double length = distance * length_fractions[i_bin];
+              double volume = length * d1 * d2;
 
-              // Add volume to result
-              result.add_volume(mesh_index, i_material, length * d1 * d2);
+              if (compute_bboxes) {
+                double axis_start = r0[axis] + distance * cumulative_frac;
+                double axis_end = axis_start + length;
+                cumulative_frac += length_fractions[i_bin];
+
+                Position contrib_min = site.r;
+                Position contrib_max = site.r;
+
+                contrib_min[ax1] = site.r[ax1] - 0.5 * d1;
+                contrib_max[ax1] = site.r[ax1] + 0.5 * d1;
+                contrib_min[ax2] = site.r[ax2] - 0.5 * d2;
+                contrib_max[ax2] = site.r[ax2] + 0.5 * d2;
+                contrib_min[axis] = std::min(axis_start, axis_end);
+                contrib_max[axis] = std::max(axis_start, axis_end);
+
+                BoundingBox contrib_bbox {contrib_min, contrib_max};
+                contrib_bbox &= bbox;
+
+                result.add_volume(
+                  mesh_index, i_material, volume, &contrib_bbox);
+              } else {
+                // Add volume to result
+                result.add_volume(mesh_index, i_material, volume);
+              }
             }
 
             if (distance == max_distance)
@@ -505,10 +639,15 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
   // Combine results from multiple MPI processes
   if (mpi::n_procs > 1) {
     int total = this->n_bins() * table_size;
+    int total_bbox = total * 6;
     if (mpi::master) {
       // Allocate temporary buffer for receiving data
-      std::vector<int32_t> mats(total);
-      std::vector<double> vols(total);
+      vector<int32_t> mats(total);
+      vector<double> vols(total);
+      vector<double> recv_bboxes;
+      if (compute_bboxes) {
+        recv_bboxes.resize(total_bbox);
+      }
 
       for (int i = 1; i < mpi::n_procs; ++i) {
         // Receive material indices and volumes from process i
@@ -516,6 +655,10 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
           MPI_STATUS_IGNORE);
         MPI_Recv(vols.data(), total, MPI_DOUBLE, i, i, mpi::intracomm,
           MPI_STATUS_IGNORE);
+        if (compute_bboxes) {
+          MPI_Recv(recv_bboxes.data(), total_bbox, MPI_DOUBLE, i, i,
+            mpi::intracomm, MPI_STATUS_IGNORE);
+        }
 
         // Combine with existing results; we can call thread unsafe version of
         // add_volume because each thread is operating on a different element
@@ -524,7 +667,18 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
           for (int k = 0; k < table_size; ++k) {
             int index = index_elem * table_size + k;
             if (mats[index] != EMPTY) {
-              result.add_volume_unsafe(index_elem, mats[index], vols[index]);
+              if (compute_bboxes) {
+                int bbox_index = index * 6;
+                BoundingBox slot_bbox {
+                  {recv_bboxes[bbox_index + 0], recv_bboxes[bbox_index + 1],
+                    recv_bboxes[bbox_index + 2]},
+                  {recv_bboxes[bbox_index + 3], recv_bboxes[bbox_index + 4],
+                    recv_bboxes[bbox_index + 5]}};
+                result.add_volume_unsafe(
+                  index_elem, mats[index], vols[index], &slot_bbox);
+              } else {
+                result.add_volume_unsafe(index_elem, mats[index], vols[index]);
+              }
             }
           }
         }
@@ -533,6 +687,9 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
       // Send material indices and volumes to process 0
       MPI_Send(materials, total, MPI_INT32_T, 0, mpi::rank, mpi::intracomm);
       MPI_Send(volumes, total, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm);
+      if (compute_bboxes) {
+        MPI_Send(bboxes, total_bbox, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm);
+      }
     }
   }
 
@@ -2428,14 +2585,14 @@ extern "C" int openmc_mesh_bounding_box(int32_t index, double* ll, double* ur)
 }
 
 extern "C" int openmc_mesh_material_volumes(int32_t index, int nx, int ny,
-  int nz, int table_size, int32_t* materials, double* volumes)
+  int nz, int table_size, int32_t* materials, double* volumes, double* bboxes)
 {
   if (int err = check_mesh(index))
     return err;
 
   try {
     model::meshes[index]->material_volumes(
-      nx, ny, nz, table_size, materials, volumes);
+      nx, ny, nz, table_size, materials, volumes, bboxes);
   } catch (const std::exception& e) {
     set_errmsg(e.what());
     if (starts_with(e.what(), "Mesh")) {
@@ -3485,9 +3642,6 @@ void LibMesh::initialize()
   // assuming that unstructured meshes used in OpenMC are 3D
   n_dimension_ = 3;
 
-  if (length_multiplier_ > 0.0) {
-    libMesh::MeshTools::Modification::scale(*m_, length_multiplier_);
-  }
   // if OpenMC is managing the libMesh::MeshBase instance, prepare the mesh.
   // Otherwise assume that it is prepared by its owning application
   if (unique_m_) {
@@ -3537,7 +3691,11 @@ Position LibMesh::centroid(int bin) const
 {
   const auto& elem = this->get_element_from_bin(bin);
   auto centroid = elem.vertex_average();
-  return {centroid(0), centroid(1), centroid(2)};
+  if (length_multiplier_ > 0.0) {
+    return length_multiplier_ * Position(centroid(0), centroid(1), centroid(2));
+  } else {
+    return {centroid(0), centroid(1), centroid(2)};
+  }
 }
 
 int LibMesh::n_vertices() const
@@ -3548,7 +3706,11 @@ int LibMesh::n_vertices() const
 Position LibMesh::vertex(int vertex_id) const
 {
   const auto node_ref = m_->node_ref(vertex_id);
-  return {node_ref(0), node_ref(1), node_ref(2)};
+  if (length_multiplier_ > 0.0) {
+    return length_multiplier_ * Position(node_ref(0), node_ref(1), node_ref(2));
+  } else {
+    return {node_ref(0), node_ref(1), node_ref(2)};
+  }
 }
 
 std::vector<int> LibMesh::connectivity(int elem_id) const
@@ -3689,6 +3851,11 @@ int LibMesh::get_bin(Position r) const
   // look-up a tet using the point locator
   libMesh::Point p(r.x, r.y, r.z);
 
+  if (length_multiplier_ > 0.0) {
+    // Scale the point down
+    p /= length_multiplier_;
+  }
+
   // quick rejection check
   if (!bbox_.contains_point(p)) {
     return -1;
@@ -3722,22 +3889,32 @@ const libMesh::Elem& LibMesh::get_element_from_bin(int bin) const
 
 double LibMesh::volume(int bin) const
 {
-  return this->get_element_from_bin(bin).volume();
+  return this->get_element_from_bin(bin).volume() * length_multiplier_ *
+         length_multiplier_ * length_multiplier_;
 }
 
-AdaptiveLibMesh::AdaptiveLibMesh(
-  libMesh::MeshBase& input_mesh, double length_multiplier)
-  : LibMesh(input_mesh, length_multiplier), num_active_(m_->n_active_elem())
+AdaptiveLibMesh::AdaptiveLibMesh(libMesh::MeshBase& input_mesh,
+  double length_multiplier,
+  const std::set<libMesh::subdomain_id_type>& block_ids)
+  : LibMesh(input_mesh, length_multiplier), block_ids_(block_ids),
+    block_restrict_(!block_ids_.empty()),
+    num_active_(
+      block_restrict_
+        ? std::distance(m_->active_subdomain_set_elements_begin(block_ids_),
+            m_->active_subdomain_set_elements_end(block_ids_))
+        : m_->n_active_elem())
 {
   // if the mesh is adaptive elements aren't guaranteed by libMesh to be
   // contiguous in ID space, so we need to map from bin indices (defined over
   // active elements) to global dof ids
   bin_to_elem_map_.reserve(num_active_);
   elem_to_bin_map_.resize(m_->n_elem(), -1);
-  for (auto it = m_->active_elements_begin(); it != m_->active_elements_end();
-       it++) {
-    auto elem = *it;
-
+  auto begin = block_restrict_
+                 ? m_->active_subdomain_set_elements_begin(block_ids_)
+                 : m_->active_elements_begin();
+  auto end = block_restrict_ ? m_->active_subdomain_set_elements_end(block_ids_)
+                             : m_->active_elements_end();
+  for (const auto& elem : libMesh::as_range(begin, end)) {
     bin_to_elem_map_.push_back(elem->id());
     elem_to_bin_map_[elem->id()] = bin_to_elem_map_.size() - 1;
   }
@@ -3770,6 +3947,27 @@ void AdaptiveLibMesh::write(const std::string& filename) const
     this->id_));
 }
 
+int AdaptiveLibMesh::get_bin(Position r) const
+{
+  // look-up a tet using the point locator
+  libMesh::Point p(r.x, r.y, r.z);
+
+  if (length_multiplier_ > 0.0) {
+    // Scale the point down
+    p /= length_multiplier_;
+  }
+
+  // quick rejection check
+  if (!bbox_.contains_point(p)) {
+    return -1;
+  }
+
+  const auto& point_locator = pl_.at(thread_num());
+
+  const auto elem_ptr = (*point_locator)(p, &block_ids_);
+  return elem_ptr ? get_bin_from_element(elem_ptr) : -1;
+}
+
 int AdaptiveLibMesh::get_bin_from_element(const libMesh::Elem* elem) const
 {
   int bin = elem_to_bin_map_[elem->id()];

src/mgxs_interface.cpp

@@ -244,7 +244,7 @@ void MgxsInterface::read_header(const std::string& path_cross_sections)
 void put_mgxs_header_data_to_globals()
 {
   // Get the minimum and maximum energies
-  int neutron = static_cast<int>(ParticleType::neutron);
+  int neutron = ParticleType::neutron().transport_index();
   data::energy_min[neutron] = data::mg.energy_bins_.back();
   data::energy_max[neutron] = data::mg.energy_bins_.front();
 

src/nuclide.cpp

@@ -379,7 +379,7 @@ void Nuclide::create_derived(
       auto pprod = xt::view(xs_[t], xt::range(j, j + n), XS_PHOTON_PROD);
 
       for (const auto& p : rx->products_) {
-        if (p.particle_ == ParticleType::photon) {
+        if (p.particle_.is_photon()) {
           for (int k = 0; k < n; ++k) {
             double E = grid_[t].energy[k + j];
 
@@ -501,7 +501,7 @@ void Nuclide::create_derived(
 
 void Nuclide::init_grid()
 {
-  int neutron = static_cast<int>(ParticleType::neutron);
+  int neutron = ParticleType::neutron().transport_index();
   double E_min = data::energy_min[neutron];
   double E_max = data::energy_max[neutron];
   int M = settings::n_log_bins;
@@ -552,7 +552,7 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const
         for (int i = 1; i < rx->products_.size(); ++i) {
           // Skip any non-neutron products
           const auto& product = rx->products_[i];
-          if (product.particle_ != ParticleType::neutron)
+          if (!product.particle_.is_neutron())
             continue;
 
           // Evaluate yield

src/output.cpp

@@ -155,21 +155,21 @@ std::string time_stamp()
 void print_particle(Particle& p)
 {
   // Display particle type and ID.
-  switch (p.type()) {
-  case ParticleType::neutron:
+  switch (p.type().pdg_number()) {
+  case PDG_NEUTRON:
     fmt::print("Neutron ");
     break;
-  case ParticleType::photon:
+  case PDG_PHOTON:
     fmt::print("Photon ");
     break;
-  case ParticleType::electron:
+  case PDG_ELECTRON:
     fmt::print("Electron ");
     break;
-  case ParticleType::positron:
+  case PDG_POSITRON:
     fmt::print("Positron ");
     break;
   default:
-    fmt::print("Unknown Particle ");
+    fmt::print("Particle {} ", p.type().str());
   }
   fmt::print("{}\n", p.id());
 

src/particle.cpp

@@ -48,17 +48,19 @@ double Particle::speed() const
   if (settings::run_CE) {
     // Determine mass in eV/c^2
     double mass;
-    switch (this->type()) {
-    case ParticleType::neutron:
+    switch (this->type().pdg_number()) {
+    case PDG_NEUTRON:
       mass = MASS_NEUTRON_EV;
       break;
-    case ParticleType::photon:
+    case PDG_PHOTON:
       mass = 0.0;
       break;
-    case ParticleType::electron:
-    case ParticleType::positron:
+    case PDG_ELECTRON:
+    case PDG_POSITRON:
       mass = MASS_ELECTRON_EV;
       break;
+    default:
+      fatal_error("Unsupported particle for speed calculation.");
     }
     // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<<m:
     return C_LIGHT * std::sqrt(this->E() * (this->E() + 2 * mass)) /
@@ -77,7 +79,11 @@ bool Particle::create_secondary(
 {
   // If energy is below cutoff for this particle, don't create secondary
   // particle
-  if (E < settings::energy_cutoff[static_cast<int>(type)]) {
+  int idx = type.transport_index();
+  if (idx == C_NONE) {
+    return false;
+  }
+  if (E < settings::energy_cutoff[idx]) {
     return false;
   }
 
@@ -235,7 +241,8 @@ void Particle::event_advance()
   boundary() = distance_to_boundary(*this);
 
   // Sample a distance to collision
-  if (type() == ParticleType::electron || type() == ParticleType::positron) {
+  if (type() == ParticleType::electron() ||
+      type() == ParticleType::positron()) {
     collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;
   } else if (macro_xs().total == 0.0) {
     collision_distance() = INFINITY;
@@ -244,7 +251,7 @@ void Particle::event_advance()
   }
 
   double speed = this->speed();
-  double time_cutoff = settings::time_cutoff[static_cast<int>(type())];
+  double time_cutoff = settings::time_cutoff[type().transport_index()];
   double distance_cutoff =
     (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
 
@@ -269,8 +276,7 @@ void Particle::event_advance()
   }
 
   // Score track-length estimate of k-eff
-  if (settings::run_mode == RunMode::EIGENVALUE &&
-      type() == ParticleType::neutron) {
+  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
     keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
   }
 
@@ -331,8 +337,7 @@ void Particle::event_cross_surface()
 void Particle::event_collide()
 {
   // Score collision estimate of keff
-  if (settings::run_mode == RunMode::EIGENVALUE &&
-      type() == ParticleType::neutron) {
+  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
     keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
   }
 
@@ -370,8 +375,7 @@ void Particle::event_collide()
     }
   }
 
-  if (!model::active_pulse_height_tallies.empty() &&
-      type() == ParticleType::photon) {
+  if (!model::active_pulse_height_tallies.empty() && type().is_photon()) {
     pht_collision_energy();
   }
 
@@ -442,7 +446,7 @@ void Particle::event_revive_from_secondary()
 
     // Subtract secondary particle energy from interim pulse-height results
     if (!model::active_pulse_height_tallies.empty() &&
-        this->type() == ParticleType::photon) {
+        this->type().is_photon()) {
       // Since the birth cell of the particle has not been set we
       // have to determine it before the energy of the secondary particle can be
       // removed from the pulse-height of this cell.
@@ -525,7 +529,7 @@ void Particle::pht_collision_energy()
 
     // If the energy of the particle is below the cutoff, it will not be sampled
     // so its energy is added to the pulse-height in the cell
-    int photon = static_cast<int>(ParticleType::photon);
+    int photon = ParticleType::photon().transport_index();
     if (E() < settings::energy_cutoff[photon]) {
       pht_storage()[index] += E();
     }
@@ -824,7 +828,7 @@ void Particle::write_restart() const
       break;
     }
     write_dataset(file_id, "id", id());
-    write_dataset(file_id, "type", static_cast<int>(type()));
+    write_dataset(file_id, "type", type().pdg_number());
 
     int64_t i = current_work();
     if (settings::run_mode == RunMode::EIGENVALUE) {
@@ -878,37 +882,6 @@ void Particle::update_neutron_xs(
 //==============================================================================
 // Non-method functions
 //==============================================================================
-
-std::string particle_type_to_str(ParticleType type)
-{
-  switch (type) {
-  case ParticleType::neutron:
-    return "neutron";
-  case ParticleType::photon:
-    return "photon";
-  case ParticleType::electron:
-    return "electron";
-  case ParticleType::positron:
-    return "positron";
-  }
-  UNREACHABLE();
-}
-
-ParticleType str_to_particle_type(std::string str)
-{
-  if (str == "neutron") {
-    return ParticleType::neutron;
-  } else if (str == "photon") {
-    return ParticleType::photon;
-  } else if (str == "electron") {
-    return ParticleType::electron;
-  } else if (str == "positron") {
-    return ParticleType::positron;
-  } else {
-    throw std::invalid_argument {fmt::format("Invalid particle name: {}", str)};
-  }
-}
-
 void add_surf_source_to_bank(Particle& p, const Surface& surf)
 {
   if (simulation::current_batch <= settings::n_inactive ||

src/particle_restart.cpp

@@ -31,6 +31,16 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode)
   hid_t file_id = file_open(settings::path_particle_restart, 'r');
 
   // Read data from file
+  bool legacy_particle_codes = true;
+  if (attribute_exists(file_id, "version")) {
+    array<int, 2> version;
+    read_attribute(file_id, "version", version);
+    if (version[0] > VERSION_PARTICLE_RESTART[0] ||
+        (version[0] == VERSION_PARTICLE_RESTART[0] && version[1] >= 1)) {
+      legacy_particle_codes = false;
+    }
+  }
+
   read_dataset(file_id, "current_batch", simulation::current_batch);
   read_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
   read_dataset(file_id, "current_generation", simulation::current_gen);
@@ -45,7 +55,8 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode)
   read_dataset(file_id, "id", p.id());
   int type;
   read_dataset(file_id, "type", type);
-  p.type() = static_cast<ParticleType>(type);
+  p.type() = legacy_particle_codes ? legacy_particle_index_to_type(type)
+                                   : ParticleType {type};
   read_dataset(file_id, "weight", p.wgt());
   read_dataset(file_id, "energy", p.E());
   read_dataset(file_id, "xyz", p.r());

src/particle_type.cpp

@@ -0,0 +1,246 @@
+#include "openmc/particle_type.h"
+
+#include <algorithm>
+#include <cctype>
+#include <stdexcept>
+
+#include "openmc/string_utils.h"
+
+namespace openmc {
+namespace {
+
+constexpr const char* ATOMIC_SYMBOL[] = {"", "H", "He", "Li", "Be", "B", "C",
+  "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca",
+  "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As",
+  "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd",
+  "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr",
+  "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf",
+  "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At",
+  "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf",
+  "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg",
+  "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og"};
+
+constexpr int MAX_Z =
+  static_cast<int>(sizeof(ATOMIC_SYMBOL) / sizeof(ATOMIC_SYMBOL[0])) - 1;
+
+bool is_integer_string(const std::string& s)
+{
+  if (s.empty())
+    return false;
+  size_t i = 0;
+  if (s[0] == '-' || s[0] == '+') {
+    if (s.size() == 1)
+      return false;
+    i = 1;
+  }
+  for (; i < s.size(); ++i) {
+    if (!std::isdigit(static_cast<unsigned char>(s[i])))
+      return false;
+  }
+  return true;
+}
+
+int atomic_number_from_symbol(std::string_view symbol)
+{
+  for (int z = 1; z <= MAX_Z; ++z) {
+    if (symbol == ATOMIC_SYMBOL[z]) {
+      return z;
+    }
+  }
+  return 0;
+}
+
+bool parse_gnds_nuclide(std::string_view name, int& Z, int& A, int& m)
+{
+  if (name.empty())
+    return false;
+
+  size_t pos = 0;
+  if (!std::isupper(static_cast<unsigned char>(name[pos])))
+    return false;
+
+  std::string symbol;
+  symbol += name[pos++];
+  if (pos < name.size() &&
+      std::islower(static_cast<unsigned char>(name[pos]))) {
+    symbol += name[pos++];
+  }
+
+  if (pos >= name.size() ||
+      !std::isdigit(static_cast<unsigned char>(name[pos]))) {
+    return false;
+  }
+
+  size_t a_start = pos;
+  while (
+    pos < name.size() && std::isdigit(static_cast<unsigned char>(name[pos]))) {
+    ++pos;
+  }
+  A = std::stoi(std::string {name.substr(a_start, pos - a_start)});
+  if (A <= 0 || A > 999)
+    return false;
+
+  m = 0;
+  if (pos < name.size()) {
+    if (name[pos] != '_' || pos + 2 >= name.size() || name[pos + 1] != 'm') {
+      return false;
+    }
+    pos += 2;
+    size_t m_start = pos;
+    while (pos < name.size() &&
+           std::isdigit(static_cast<unsigned char>(name[pos]))) {
+      ++pos;
+    }
+    if (m_start == pos)
+      return false;
+    m = std::stoi(std::string {name.substr(m_start, pos - m_start)});
+    if (m < 0 || m > 9)
+      return false;
+  }
+
+  if (pos != name.size())
+    return false;
+
+  Z = atomic_number_from_symbol(symbol);
+  return Z != 0;
+}
+
+// Helper to convert nuclear PDG number to nuclide name
+std::string nuclide_name_from_pdg(int32_t pdg)
+{
+  int32_t code = pdg;
+  int m = code % 10;
+  int A = (code / 10) % 1000;
+  int Z = (code / 10000) % 1000;
+
+  if (Z <= 0 || Z > MAX_Z || A <= 0 || A > 999) {
+    throw std::invalid_argument {
+      "Invalid nuclear PDG number: " + std::to_string(pdg)};
+  }
+
+  std::string name = ATOMIC_SYMBOL[Z] + std::to_string(A);
+  if (m > 0) {
+    name += "_m" + std::to_string(m);
+  }
+  return name;
+}
+
+} // namespace
+
+//==============================================================================
+// ParticleType member function implementations
+//==============================================================================
+
+ParticleType::ParticleType(std::string_view str)
+{
+  std::string s {str};
+  strtrim(s);
+  if (s.empty()) {
+    throw std::invalid_argument {"Particle string is empty."};
+  }
+
+  std::string lower = s;
+  to_lower(lower);
+
+  // Check for pdg: prefix
+  if (starts_with(lower, "pdg:")) {
+    std::string value_str = lower.substr(4);
+    if (!is_integer_string(value_str)) {
+      throw std::invalid_argument {"Invalid PDG number: " + value_str};
+    }
+    pdg_number_ = std::stoi(value_str);
+    return;
+  }
+
+  // Check for known particle names
+  if (lower == "neutron" || lower == "n") {
+    pdg_number_ = PDG_NEUTRON;
+    return;
+  }
+  if (lower == "photon" || lower == "gamma") {
+    pdg_number_ = PDG_PHOTON;
+    return;
+  }
+  if (lower == "electron") {
+    pdg_number_ = PDG_ELECTRON;
+    return;
+  }
+  if (lower == "positron") {
+    pdg_number_ = PDG_POSITRON;
+    return;
+  }
+  if (lower == "proton" || lower == "p" || lower == "h1") {
+    pdg_number_ = PDG_PROTON;
+    return;
+  }
+  if (lower == "deuteron" || lower == "d" || lower == "h2") {
+    pdg_number_ = PDG_DEUTERON;
+    return;
+  }
+  if (lower == "triton" || lower == "t" || lower == "h3") {
+    pdg_number_ = PDG_TRITON;
+    return;
+  }
+  if (lower == "alpha" || lower == "he4") {
+    pdg_number_ = PDG_ALPHA;
+    return;
+  }
+
+  // Check for integer string
+  if (is_integer_string(s)) {
+    pdg_number_ = std::stoi(s);
+    return;
+  }
+
+  // Try to parse as GNDS nuclide name
+  int Z = 0;
+  int A = 0;
+  int m = 0;
+  if (!parse_gnds_nuclide(s, Z, A, m)) {
+    throw std::invalid_argument {"Invalid nuclide name: " + s};
+  }
+  pdg_number_ = 1000000000 + Z * 10000 + A * 10 + m;
+}
+
+std::string ParticleType::str() const
+{
+  if (pdg_number_ == PDG_NEUTRON)
+    return "neutron";
+  if (pdg_number_ == PDG_PHOTON)
+    return "photon";
+  if (pdg_number_ == PDG_ELECTRON)
+    return "electron";
+  if (pdg_number_ == PDG_POSITRON)
+    return "positron";
+  if (pdg_number_ == PDG_PROTON)
+    return "proton";
+
+  if (is_nucleus()) {
+    return nuclide_name_from_pdg(pdg_number_);
+  }
+
+  return "pdg:" + std::to_string(pdg_number_);
+}
+
+//==============================================================================
+// Free function implementations
+//==============================================================================
+
+ParticleType legacy_particle_index_to_type(int index)
+{
+  switch (index) {
+  case 0:
+    return ParticleType {PDG_NEUTRON};
+  case 1:
+    return ParticleType {PDG_PHOTON};
+  case 2:
+    return ParticleType {PDG_ELECTRON};
+  case 3:
+    return ParticleType {PDG_POSITRON};
+  default:
+    throw std::invalid_argument {
+      "Invalid legacy particle index: " + std::to_string(index)};
+  }
+}
+
+} // namespace openmc

src/photon.cpp

@@ -293,7 +293,7 @@ PhotonInteraction::PhotonInteraction(hid_t group)
     close_group(rgroup);
 
     // Truncate the bremsstrahlung data at the cutoff energy
-    int photon = static_cast<int>(ParticleType::photon);
+    int photon = ParticleType::photon().transport_index();
     const auto& E {electron_energy};
     double cutoff = settings::energy_cutoff[photon];
     if (cutoff > E(0)) {
@@ -805,7 +805,7 @@ void PhotonInteraction::atomic_relaxation(int i_shell, Particle& p) const
     if (shell.transitions.empty()) {
       Direction u = isotropic_direction(p.current_seed());
       double E = shell.binding_energy;
-      p.create_secondary(p.wgt(), u, E, ParticleType::photon);
+      p.create_secondary(p.wgt(), u, E, ParticleType::photon());
       continue;
     }
 
@@ -833,12 +833,13 @@ void PhotonInteraction::atomic_relaxation(int i_shell, Particle& p) const
       holes[n_holes++] = transition.secondary_subshell;
 
       // Create auger electron
-      p.create_secondary(p.wgt(), u, transition.energy, ParticleType::electron);
+      p.create_secondary(
+        p.wgt(), u, transition.energy, ParticleType::electron());
     } else {
       // Radiative transition -- get X-ray energy
 
       // Create fluorescent photon
-      p.create_secondary(p.wgt(), u, transition.energy, ParticleType::photon);
+      p.create_secondary(p.wgt(), u, transition.energy, ParticleType::photon());
     }
   }
 }

src/physics.cpp

@@ -46,27 +46,29 @@ void collision(Particle& p)
   ++(p.n_collision());
 
   // Sample reaction for the material the particle is in
-  switch (p.type()) {
-  case ParticleType::neutron:
+  switch (p.type().pdg_number()) {
+  case PDG_NEUTRON:
     sample_neutron_reaction(p);
     break;
-  case ParticleType::photon:
+  case PDG_PHOTON:
     sample_photon_reaction(p);
     break;
-  case ParticleType::electron:
+  case PDG_ELECTRON:
     sample_electron_reaction(p);
     break;
-  case ParticleType::positron:
+  case PDG_POSITRON:
     sample_positron_reaction(p);
     break;
+  default:
+    fatal_error("Unsupported particle PDG for collision sampling.");
   }
 
   if (settings::weight_window_checkpoint_collision)
     apply_weight_windows(p);
 
   // Kill particle if energy falls below cutoff
-  int type = static_cast<int>(p.type());
-  if (p.E() < settings::energy_cutoff[type]) {
+  int type = p.type().transport_index();
+  if (type != C_NONE && p.E() < settings::energy_cutoff[type]) {
     p.wgt() = 0.0;
   }
 
@@ -75,11 +77,11 @@ void collision(Particle& p)
     std::string msg;
     if (p.event() == TallyEvent::KILL) {
       msg = fmt::format("    Killed. Energy = {} eV.", p.E());
-    } else if (p.type() == ParticleType::neutron) {
+    } else if (p.type().is_neutron()) {
       msg = fmt::format("    {} with {}. Energy = {} eV.",
         reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_,
         p.E());
-    } else if (p.type() == ParticleType::photon) {
+    } else if (p.type().is_photon()) {
       msg = fmt::format("    {} with {}. Energy = {} eV.",
         reaction_name(p.event_mt()),
         to_element(data::nuclides[p.event_nuclide()]->name_), p.E());
@@ -208,7 +210,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
     // Initialize fission site object with particle data
     SourceSite site;
     site.r = p.r();
-    site.particle = ParticleType::neutron;
+    site.particle = ParticleType::neutron();
     site.time = p.time();
     site.wgt = 1. / weight;
     site.surf_id = 0;
@@ -218,7 +220,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
 
     // Reject site if it exceeds time cutoff
     if (site.delayed_group > 0) {
-      double t_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      double t_cutoff = settings::time_cutoff[site.particle.transport_index()];
       if (site.time > t_cutoff) {
         continue;
       }
@@ -287,7 +289,7 @@ void sample_photon_reaction(Particle& p)
   // Kill photon if below energy cutoff -- an extra check is made here because
   // photons with energy below the cutoff may have been produced by neutrons
   // reactions or atomic relaxation
-  int photon = static_cast<int>(ParticleType::photon);
+  int photon = ParticleType::photon().transport_index();
   if (p.E() < settings::energy_cutoff[photon]) {
     p.E() = 0.0;
     p.wgt() = 0.0;
@@ -337,13 +339,13 @@ void sample_photon_reaction(Particle& p)
     // Create Compton electron
     double phi = uniform_distribution(0., 2.0 * PI, p.current_seed());
     double E_electron = (alpha - alpha_out) * MASS_ELECTRON_EV - e_b;
-    int electron = static_cast<int>(ParticleType::electron);
+    int electron = ParticleType::electron().transport_index();
     if (E_electron >= settings::energy_cutoff[electron]) {
       double mu_electron = (alpha - alpha_out * p.mu()) /
                            std::sqrt(alpha * alpha + alpha_out * alpha_out -
                                      2.0 * alpha * alpha_out * p.mu());
       Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed());
-      p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
+      p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron());
     }
 
     // Allow electrons to fill orbital and produce Auger electrons and
@@ -416,7 +418,7 @@ void sample_photon_reaction(Particle& p)
         u.z = std::sqrt(1.0 - mu * mu) * std::sin(phi);
 
         // Create secondary electron
-        p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
+        p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron());
 
         // Allow electrons to fill orbital and produce auger electrons
         // and fluorescent photons
@@ -441,12 +443,11 @@ void sample_photon_reaction(Particle& p)
 
     // Create secondary electron
     Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed());
-    p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
+    p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron());
 
     // Create secondary positron
     u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed());
-    p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron);
-
+    p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron());
     p.event() = TallyEvent::ABSORB;
     p.event_mt() = PAIR_PROD;
     p.wgt() = 0.0;
@@ -481,8 +482,8 @@ void sample_positron_reaction(Particle& p)
   Direction u = isotropic_direction(p.current_seed());
 
   // Create annihilation photon pair traveling in opposite directions
-  p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon);
-  p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon);
+  p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon());
+  p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon());
 
   p.E() = 0.0;
   p.wgt() = 0.0;
@@ -607,7 +608,7 @@ void sample_photon_product(
       continue;
 
     for (int j = 0; j < rx->products_.size(); ++j) {
-      if (rx->products_[j].particle_ == ParticleType::photon) {
+      if (rx->products_[j].particle_.is_photon()) {
         // For fission, artificially increase the photon yield to account
         // for delayed photons
         double f = 1.0;
@@ -1096,7 +1097,7 @@ void sample_fission_neutron(
     rx.products_[site->delayed_group].sample(E_in, site->E, mu, seed);
 
     // resample if energy is greater than maximum neutron energy
-    constexpr int neutron = static_cast<int>(ParticleType::neutron);
+    int neutron = ParticleType::neutron().transport_index();
     if (site->E < data::energy_max[neutron])
       break;
 
@@ -1156,7 +1157,7 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
   if (std::floor(yield) == yield && yield > 0) {
     // If yield is integral, create exactly that many secondary particles
     for (int i = 0; i < static_cast<int>(std::round(yield)) - 1; ++i) {
-      p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron);
+      p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron());
     }
   } else {
     // Otherwise, change weight of particle based on yield
@@ -1212,7 +1213,7 @@ void sample_secondary_photons(Particle& p, int i_nuclide)
     }
 
     // Create the secondary photon
-    bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon);
+    bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon());
 
     // Tag secondary particle with parent nuclide
     if (created_photon && settings::use_decay_photons) {

src/physics_mg.cpp

@@ -136,7 +136,7 @@ void create_fission_sites(Particle& p)
     // Initialize fission site object with particle data
     SourceSite site;
     site.r = p.r();
-    site.particle = ParticleType::neutron;
+    site.particle = ParticleType::neutron();
     site.time = p.time();
     site.wgt = 1. / weight;
 
@@ -171,7 +171,7 @@ void create_fission_sites(Particle& p)
       site.time -= std::log(prn(p.current_seed())) / decay_rate;
 
       // Reject site if it exceeds time cutoff
-      double t_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      double t_cutoff = settings::time_cutoff[site.particle.transport_index()];
       if (site.time > t_cutoff) {
         continue;
       }

src/random_ray/flat_source_domain.cpp

@@ -109,18 +109,21 @@ void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
 
   // Add scattering + fission source
   int material = srh.material();
+  double density_mult = srh.density_mult();
   if (material != MATERIAL_VOID) {
     double inverse_k_eff = 1.0 / k_eff_;
     for (int g_out = 0; g_out < negroups_; g_out++) {
-      double sigma_t = sigma_t_[material * negroups_ + g_out];
+      double sigma_t = sigma_t_[material * negroups_ + g_out] * density_mult;
       double scatter_source = 0.0;
       double fission_source = 0.0;
 
       for (int g_in = 0; g_in < negroups_; g_in++) {
         double scalar_flux = srh.scalar_flux_old(g_in);
-        double sigma_s =
-          sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];
-        double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];
+        double sigma_s = sigma_s_[material * negroups_ * negroups_ +
+                                  g_out * negroups_ + g_in] *
+                         density_mult;
+        double nu_sigma_f =
+          nu_sigma_f_[material * negroups_ + g_in] * density_mult;
         double chi = chi_[material * negroups_ + g_out];
 
         scatter_source += sigma_s * scalar_flux;
@@ -198,7 +201,8 @@ void FlatSourceDomain::set_flux_to_flux_plus_source(
         source_regions_.volume_sq(sr);
     }
   } else {
-    double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g];
+    double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g] *
+                     source_regions_.density_mult(sr);
     source_regions_.scalar_flux_new(sr, g) /= (sigma_t * volume);
     source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g);
   }
@@ -332,7 +336,8 @@ void FlatSourceDomain::compute_k_eff()
     double sr_fission_source_new = 0;
 
     for (int g = 0; g < negroups_; g++) {
-      double nu_sigma_f = nu_sigma_f_[material * negroups_ + g];
+      double nu_sigma_f = nu_sigma_f_[material * negroups_ + g] *
+                          source_regions_.density_mult(sr);
       sr_fission_source_old +=
         nu_sigma_f * source_regions_.scalar_flux_old(sr, g);
       sr_fission_source_new +=
@@ -562,7 +567,8 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const
       // to get the total source strength in the expected units.
       double sigma_t = 1.0;
       if (material != MATERIAL_VOID) {
-        sigma_t = sigma_t_[material * negroups_ + g];
+        sigma_t =
+          sigma_t_[material * negroups_ + g] * source_regions_.density_mult(sr);
       }
       simulation_external_source_strength +=
         source_regions_.external_source(sr, g) * sigma_t * volume;
@@ -618,7 +624,9 @@ void FlatSourceDomain::random_ray_tally()
     // source strength.
     double volume = source_regions_.volume(sr) * simulation_volume_;
 
-    double material = source_regions_.material(sr);
+    int material = source_regions_.material(sr);
+    double density_mult = source_regions_.density_mult(sr);
+
     for (int g = 0; g < negroups_; g++) {
       double flux =
         source_regions_.scalar_flux_new(sr, g) * source_normalization_factor;
@@ -634,19 +642,22 @@ void FlatSourceDomain::random_ray_tally()
 
         case SCORE_TOTAL:
           if (material != MATERIAL_VOID) {
-            score = flux * volume * sigma_t_[material * negroups_ + g];
+            score =
+              flux * volume * sigma_t_[material * negroups_ + g] * density_mult;
           }
           break;
 
         case SCORE_FISSION:
           if (material != MATERIAL_VOID) {
-            score = flux * volume * sigma_f_[material * negroups_ + g];
+            score =
+              flux * volume * sigma_f_[material * negroups_ + g] * density_mult;
           }
           break;
 
         case SCORE_NU_FISSION:
           if (material != MATERIAL_VOID) {
-            score = flux * volume * nu_sigma_f_[material * negroups_ + g];
+            score = flux * volume * nu_sigma_f_[material * negroups_ + g] *
+                    density_mult;
           }
           break;
 
@@ -654,10 +665,15 @@ void FlatSourceDomain::random_ray_tally()
           score = 1.0;
           break;
 
+        case SCORE_KAPPA_FISSION:
+          score = flux * volume * kappa_fission_[material * negroups_ + g] *
+                  density_mult;
+          break;
+
         default:
           fatal_error("Invalid score specified in tallies.xml. Only flux, "
-                      "total, fission, nu-fission, and events are supported in "
-                      "random ray mode.");
+                      "total, fission, nu-fission, kappa-fission, and events "
+                      "are supported in random ray mode.");
           break;
         }
         // Apply score to the appropriate tally bin
@@ -913,7 +929,8 @@ void FlatSourceDomain::output_to_vtk() const
             for (int g = 0; g < negroups_; g++) {
               int64_t source_element = fsr * negroups_ + g;
               float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g);
-              double sigma_f = sigma_f_[mat * negroups_ + g];
+              double sigma_f = sigma_f_[mat * negroups_ + g] *
+                               source_regions_.density_mult(fsr);
               total_fission += sigma_f * flux;
             }
           }
@@ -934,7 +951,8 @@ void FlatSourceDomain::output_to_vtk() const
             // multiply it back to get the true external source.
             double sigma_t = 1.0;
             if (mat != MATERIAL_VOID) {
-              sigma_t = sigma_t_[mat * negroups_ + g];
+              sigma_t = sigma_t_[mat * negroups_ + g] *
+                        source_regions_.density_mult(fsr);
             }
             total_external += source_regions_.external_source(fsr, g) * sigma_t;
           }
@@ -1152,6 +1170,10 @@ void FlatSourceDomain::flatten_xs()
         }
         chi_.push_back(chi);
 
+        double kappa_fission =
+          m.get_xs(MgxsType::KAPPA_FISSION, g_out, NULL, NULL, NULL, t, a);
+        kappa_fission_.push_back(kappa_fission);
+
         for (int g_in = 0; g_in < negroups_; g_in++) {
           double sigma_s =
             m.get_xs(MgxsType::NU_SCATTER, g_in, &g_out, NULL, NULL, t, a);
@@ -1167,6 +1189,7 @@ void FlatSourceDomain::flatten_xs()
         nu_sigma_f_.push_back(0);
         sigma_f_.push_back(0);
         chi_.push_back(0);
+        kappa_fission_.push_back(0);
         for (int g_in = 0; g_in < negroups_; g_in++) {
           sigma_s_.push_back(0);
         }
@@ -1244,7 +1267,8 @@ void FlatSourceDomain::set_adjoint_sources()
       continue;
     }
     for (int g = 0; g < negroups_; g++) {
-      double sigma_t = sigma_t_[material * negroups_ + g];
+      double sigma_t =
+        sigma_t_[material * negroups_ + g] * source_regions_.density_mult(sr);
       source_regions_.external_source(sr, g) /= sigma_t;
     }
   }
@@ -1495,6 +1519,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
 
   handle.material() = material;
 
+  handle.density_mult() = cell.density_mult(gs.cell_instance());
+
   // Store the mesh index (if any) assigned to this source region
   handle.mesh() = mesh_idx;
 
@@ -1523,7 +1549,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
     // Divide external source term by sigma_t
     if (material != C_NONE) {
       for (int g = 0; g < negroups_; g++) {
-        double sigma_t = sigma_t_[material * negroups_ + g];
+        double sigma_t =
+          sigma_t_[material * negroups_ + g] * handle.density_mult();
         handle.external_source(g) /= sigma_t;
       }
     }
@@ -1598,6 +1625,7 @@ void FlatSourceDomain::apply_transport_stabilization()
 #pragma omp parallel for
   for (int64_t sr = 0; sr < n_source_regions(); sr++) {
     int material = source_regions_.material(sr);
+    double density_mult = source_regions_.density_mult(sr);
     if (material == MATERIAL_VOID) {
       continue;
     }
@@ -1605,9 +1633,10 @@ void FlatSourceDomain::apply_transport_stabilization()
       // Only apply stabilization if the diagonal (in-group) scattering XS is
       // negative
       double sigma_s =
-        sigma_s_[material * negroups_ * negroups_ + g * negroups_ + g];
+        sigma_s_[material * negroups_ * negroups_ + g * negroups_ + g] *
+        density_mult;
       if (sigma_s < 0.0) {
-        double sigma_t = sigma_t_[material * negroups_ + g];
+        double sigma_t = sigma_t_[material * negroups_ + g] * density_mult;
         double phi_new = source_regions_.scalar_flux_new(sr, g);
         double phi_old = source_regions_.scalar_flux_old(sr, g);
 

src/random_ray/linear_source_domain.cpp

@@ -43,12 +43,13 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
 
   // Add scattering + fission source
   int material = srh.material();
+  double density_mult = srh.density_mult();
   if (material != MATERIAL_VOID) {
     double inverse_k_eff = 1.0 / k_eff_;
     MomentMatrix invM = srh.mom_matrix().inverse();
 
     for (int g_out = 0; g_out < negroups_; g_out++) {
-      double sigma_t = sigma_t_[material * negroups_ + g_out];
+      double sigma_t = sigma_t_[material * negroups_ + g_out] * density_mult;
 
       double scatter_flat = 0.0f;
       double fission_flat = 0.0f;
@@ -61,9 +62,11 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
         MomentArray flux_linear = srh.flux_moments_old(g_in);
 
         // Handles for cross sections
-        double sigma_s =
-          sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];
-        double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];
+        double sigma_s = sigma_s_[material * negroups_ * negroups_ +
+                                  g_out * negroups_ + g_in] *
+                         density_mult;
+        double nu_sigma_f =
+          nu_sigma_f_[material * negroups_ + g_in] * density_mult;
         double chi = chi_[material * negroups_ + g_out];
 
         // Compute source terms for flat and linear components of the flux

src/random_ray/random_ray.cpp

@@ -435,7 +435,8 @@ void RandomRay::attenuate_flux_flat_source(
 
   // MOC incoming flux attenuation + source contribution/attenuation equation
   for (int g = 0; g < negroups_; g++) {
-    float sigma_t = domain_->sigma_t_[material * negroups_ + g];
+    float sigma_t =
+      domain_->sigma_t_[material * negroups_ + g] * srh.density_mult();
     float tau = sigma_t * distance;
     float exponential = cjosey_exponential(tau); // exponential = 1 - exp(-tau)
     float new_delta_psi = (angular_flux_[g] - srh.source(g)) * exponential;
@@ -558,7 +559,8 @@ void RandomRay::attenuate_flux_linear_source(
   for (int g = 0; g < negroups_; g++) {
 
     // Compute tau, the optical thickness of the ray segment
-    float sigma_t = domain_->sigma_t_[material * negroups_ + g];
+    float sigma_t =
+      domain_->sigma_t_[material * negroups_ + g] * srh.density_mult();
     float tau = sigma_t * distance;
 
     // If tau is very small, set it to zero to avoid numerical issues.

src/random_ray/random_ray_simulation.cpp

@@ -28,17 +28,13 @@ void openmc_run_random_ray()
   // Run forward simulation
   //////////////////////////////////////////////////////////
 
-  // Check if adjoint calculation is needed. If it is, we will run the forward
-  // calculation first and then the adjoint calculation later.
-  bool adjoint_needed = FlatSourceDomain::adjoint_;
-
-  // Configure the domain for forward simulation
-  FlatSourceDomain::adjoint_ = false;
-
-  // If we're going to do an adjoint simulation afterwards, report that this is
-  // the initial forward flux solve.
-  if (adjoint_needed && mpi::master)
-    header("FORWARD FLUX SOLVE", 3);
+  if (mpi::master) {
+    if (FlatSourceDomain::adjoint_) {
+      FlatSourceDomain::adjoint_ = false;
+      openmc::print_adjoint_header();
+      FlatSourceDomain::adjoint_ = true;
+    }
+  }
 
   // Initialize OpenMC general data structures
   openmc_simulation_init();
@@ -53,76 +49,20 @@ void openmc_run_random_ray()
   // Initialize fixed sources, if present
   sim.apply_fixed_sources_and_mesh_domains();
 
-  // Begin main simulation timer
-  simulation::time_total.start();
-
-  // Execute random ray simulation
+  // Run initial random ray simulation
   sim.simulate();
 
-  // End main simulation timer
-  simulation::time_total.stop();
-
-  // Normalize and save the final forward flux
-  double source_normalization_factor =
-    sim.domain()->compute_fixed_source_normalization_factor() /
-    (settings::n_batches - settings::n_inactive);
-
-#pragma omp parallel for
-  for (uint64_t se = 0; se < sim.domain()->n_source_elements(); se++) {
-    sim.domain()->source_regions_.scalar_flux_final(se) *=
-      source_normalization_factor;
-  }
-
-  // Finalize OpenMC
-  openmc_simulation_finalize();
-
-  // Output all simulation results
-  sim.output_simulation_results();
-
   //////////////////////////////////////////////////////////
   // Run adjoint simulation (if enabled)
   //////////////////////////////////////////////////////////
 
-  if (!adjoint_needed) {
-    return;
+  if (sim.adjoint_needed_) {
+    // Setup for adjoint simulation
+    sim.prepare_adjoint_simulation();
+
+    // Run adjoint simulation
+    sim.simulate();
   }
-
-  reset_timers();
-
-  // Configure the domain for adjoint simulation
-  FlatSourceDomain::adjoint_ = true;
-
-  if (mpi::master)
-    header("ADJOINT FLUX SOLVE", 3);
-
-  // Initialize OpenMC general data structures
-  openmc_simulation_init();
-
-  sim.domain()->k_eff_ = 1.0;
-
-  // Initialize adjoint fixed sources, if present
-  sim.prepare_fixed_sources_adjoint();
-
-  // Transpose scattering matrix
-  sim.domain()->transpose_scattering_matrix();
-
-  // Swap nu_sigma_f and chi
-  sim.domain()->nu_sigma_f_.swap(sim.domain()->chi_);
-
-  // Begin main simulation timer
-  simulation::time_total.start();
-
-  // Execute random ray simulation
-  sim.simulate();
-
-  // End main simulation timer
-  simulation::time_total.stop();
-
-  // Finalize OpenMC
-  openmc_simulation_finalize();
-
-  // Output all simulation results
-  sim.output_simulation_results();
 }
 
 // Enforces restrictions on inputs in random ray mode.  While there are
@@ -143,11 +83,12 @@ void validate_random_ray_inputs()
       case SCORE_FISSION:
       case SCORE_NU_FISSION:
       case SCORE_EVENTS:
+      case SCORE_KAPPA_FISSION:
         break;
       default:
         fatal_error(
-          "Invalid score specified. Only flux, total, fission, nu-fission, and "
-          "event scores are supported in random ray mode.");
+          "Invalid score specified. Only flux, total, fission, nu-fission, "
+          "kappa-fission, and event scores are supported in random ray mode.");
       }
     }
 
@@ -353,6 +294,17 @@ void openmc_reset_random_ray()
   RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
 }
 
+void print_adjoint_header()
+{
+  if (!FlatSourceDomain::adjoint_)
+    // If we're going to do an adjoint simulation afterwards, report that this
+    // is the initial forward flux solve.
+    header("FORWARD FLUX SOLVE", 3);
+  else
+    // Otherwise report that we are doing the adjoint simulation
+    header("ADJOINT FLUX SOLVE", 3);
+}
+
 //==============================================================================
 // RandomRaySimulation implementation
 //==============================================================================
@@ -383,6 +335,16 @@ RandomRaySimulation::RandomRaySimulation()
   // Convert OpenMC native MGXS into a more efficient format
   // internal to the random ray solver
   domain_->flatten_xs();
+
+  // Check if adjoint calculation is needed. If it is, we will run the forward
+  // calculation first and then the adjoint calculation later.
+  adjoint_needed_ = FlatSourceDomain::adjoint_;
+
+  // Adjoint is always false for the forward calculation
+  FlatSourceDomain::adjoint_ = false;
+
+  // The first simulation is run after initialization
+  is_first_simulation_ = true;
 }
 
 void RandomRaySimulation::apply_fixed_sources_and_mesh_domains()
@@ -403,8 +365,41 @@ void RandomRaySimulation::prepare_fixed_sources_adjoint()
   }
 }
 
+void RandomRaySimulation::prepare_adjoint_simulation()
+{
+  // Configure the domain for adjoint simulation
+  FlatSourceDomain::adjoint_ = true;
+
+  // Reset k-eff
+  domain_->k_eff_ = 1.0;
+
+  // Initialize adjoint fixed sources, if present
+  prepare_fixed_sources_adjoint();
+
+  // Transpose scattering matrix
+  domain_->transpose_scattering_matrix();
+
+  // Swap nu_sigma_f and chi
+  domain_->nu_sigma_f_.swap(domain_->chi_);
+}
+
 void RandomRaySimulation::simulate()
 {
+  if (!is_first_simulation_) {
+    if (mpi::master && adjoint_needed_)
+      openmc::print_adjoint_header();
+
+    // Reset the timers and reinitialize the general OpenMC datastructures if
+    // this is after the first simulation
+    reset_timers();
+
+    // Initialize OpenMC general data structures
+    openmc_simulation_init();
+  }
+
+  // Begin main simulation timer
+  simulation::time_total.start();
+
   // Random ray power iteration loop
   while (simulation::current_batch < settings::n_batches) {
     // Initialize the current batch
@@ -493,6 +488,31 @@ void RandomRaySimulation::simulate()
   } // End random ray power iteration loop
 
   domain_->count_external_source_regions();
+
+  // End main simulation timer
+  simulation::time_total.stop();
+
+  // Normalize and save the final flux
+  double source_normalization_factor =
+    domain_->compute_fixed_source_normalization_factor() /
+    (settings::n_batches - settings::n_inactive);
+
+#pragma omp parallel for
+  for (uint64_t se = 0; se < domain_->n_source_elements(); se++) {
+    domain_->source_regions_.scalar_flux_final(se) *=
+      source_normalization_factor;
+  }
+
+  // Finalize OpenMC
+  openmc_simulation_finalize();
+
+  // Output all simulation results
+  output_simulation_results();
+
+  // Toggle that the simulation object has been initialized after the first
+  // simulation
+  if (is_first_simulation_)
+    is_first_simulation_ = false;
 }
 
 void RandomRaySimulation::output_simulation_results() const

src/random_ray/source_region.cpp

@@ -11,10 +11,11 @@ namespace openmc {
 //==============================================================================
 SourceRegionHandle::SourceRegionHandle(SourceRegion& sr)
   : negroups_(sr.scalar_flux_old_.size()), material_(&sr.material_),
-    is_small_(&sr.is_small_), n_hits_(&sr.n_hits_),
-    is_linear_(sr.source_gradients_.size() > 0), lock_(&sr.lock_),
-    volume_(&sr.volume_), volume_t_(&sr.volume_t_), volume_sq_(&sr.volume_sq_),
-    volume_sq_t_(&sr.volume_sq_t_), volume_naive_(&sr.volume_naive_),
+    density_mult_(&sr.density_mult_), is_small_(&sr.is_small_),
+    n_hits_(&sr.n_hits_), is_linear_(sr.source_gradients_.size() > 0),
+    lock_(&sr.lock_), volume_(&sr.volume_), volume_t_(&sr.volume_t_),
+    volume_sq_(&sr.volume_sq_), volume_sq_t_(&sr.volume_sq_t_),
+    volume_naive_(&sr.volume_naive_),
     position_recorded_(&sr.position_recorded_),
     external_source_present_(&sr.external_source_present_),
     position_(&sr.position_), centroid_(&sr.centroid_),
@@ -70,6 +71,7 @@ void SourceRegionContainer::push_back(const SourceRegion& sr)
 
   // Scalar fields
   material_.push_back(sr.material_);
+  density_mult_.push_back(sr.density_mult_);
   is_small_.push_back(sr.is_small_);
   n_hits_.push_back(sr.n_hits_);
   lock_.push_back(sr.lock_);
@@ -123,6 +125,7 @@ void SourceRegionContainer::assign(
   // Clear existing data
   n_source_regions_ = 0;
   material_.clear();
+  density_mult_.clear();
   is_small_.clear();
   n_hits_.clear();
   lock_.clear();
@@ -180,6 +183,7 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr)
   SourceRegionHandle handle;
   handle.negroups_ = negroups();
   handle.material_ = &material(sr);
+  handle.density_mult_ = &density_mult(sr);
   handle.is_small_ = &is_small(sr);
   handle.n_hits_ = &n_hits(sr);
   handle.is_linear_ = is_linear();

src/reaction.cpp

@@ -70,9 +70,8 @@ Reaction::Reaction(
 
   if (settings::use_decay_photons) {
     // Remove photon products for D1S method
-    products_.erase(
-      std::remove_if(products_.begin(), products_.end(),
-        [](const auto& p) { return p.particle_ == ParticleType::photon; }),
+    products_.erase(std::remove_if(products_.begin(), products_.end(),
+                      [](const auto& p) { return p.particle_.is_photon(); }),
       products_.end());
 
     // Determine product for D1S method

src/reaction_product.cpp

@@ -26,7 +26,7 @@ ReactionProduct::ReactionProduct(hid_t group)
   // Read particle type
   std::string temp;
   read_attribute(group, "particle", temp);
-  particle_ = str_to_particle_type(temp);
+  particle_ = ParticleType {temp};
 
   // Read emission mode and decay rate
   read_attribute(group, "emission_mode", temp);
@@ -42,7 +42,7 @@ ReactionProduct::ReactionProduct(hid_t group)
   if (emission_mode_ == EmissionMode::delayed) {
     if (attribute_exists(group, "decay_rate")) {
       read_attribute(group, "decay_rate", decay_rate_);
-    } else if (particle_ == ParticleType::neutron) {
+    } else if (particle_.is_neutron()) {
       warning(fmt::format("Decay rate doesn't exist for delayed neutron "
                           "emission ({}).",
         object_name(group)));
@@ -85,7 +85,7 @@ ReactionProduct::ReactionProduct(hid_t group)
 
 ReactionProduct::ReactionProduct(const ChainNuclide::Product& product)
 {
-  particle_ = ParticleType::photon;
+  particle_ = ParticleType::photon();
   emission_mode_ = EmissionMode::delayed;
 
   // Get chain nuclide object for radionuclide

src/simulation.cpp

@@ -701,7 +701,7 @@ void initialize_data()
 
   for (const auto& nuc : data::nuclides) {
     if (nuc->grid_.size() >= 1) {
-      int neutron = static_cast<int>(ParticleType::neutron);
+      int neutron = ParticleType::neutron().transport_index();
       data::energy_min[neutron] =
         std::max(data::energy_min[neutron], nuc->grid_[0].energy.front());
       data::energy_max[neutron] =
@@ -712,7 +712,7 @@ void initialize_data()
   if (settings::photon_transport) {
     for (const auto& elem : data::elements) {
       if (elem->energy_.size() >= 1) {
-        int photon = static_cast<int>(ParticleType::photon);
+        int photon = ParticleType::photon().transport_index();
         int n = elem->energy_.size();
         data::energy_min[photon] =
           std::max(data::energy_min[photon], std::exp(elem->energy_(1)));
@@ -725,9 +725,9 @@ void initialize_data()
       // Determine if minimum/maximum energy for bremsstrahlung is greater/less
       // than the current minimum/maximum
       if (data::ttb_e_grid.size() >= 1) {
-        int photon = static_cast<int>(ParticleType::photon);
-        int electron = static_cast<int>(ParticleType::electron);
-        int positron = static_cast<int>(ParticleType::positron);
+        int photon = ParticleType::photon().transport_index();
+        int electron = ParticleType::electron().transport_index();
+        int positron = ParticleType::positron().transport_index();
         int n_e = data::ttb_e_grid.size();
 
         const std::vector<int> charged = {electron, positron};
@@ -751,7 +751,7 @@ void initialize_data()
     // grid has not been allocated
     if (nuc->grid_.size() > 0) {
       double max_E = nuc->grid_[0].energy.back();
-      int neutron = static_cast<int>(ParticleType::neutron);
+      int neutron = ParticleType::neutron().transport_index();
       if (max_E == data::energy_max[neutron]) {
         write_message(7, "Maximum neutron transport energy: {} eV for {}",
           data::energy_max[neutron], nuc->name_);
@@ -768,7 +768,7 @@ void initialize_data()
   for (auto& nuc : data::nuclides) {
     nuc->init_grid();
   }
-  int neutron = static_cast<int>(ParticleType::neutron);
+  int neutron = ParticleType::neutron().transport_index();
   simulation::log_spacing =
     std::log(data::energy_max[neutron] / data::energy_min[neutron]) /
     settings::n_log_bins;

src/source.cpp

@@ -37,6 +37,20 @@
 
 namespace openmc {
 
+namespace {
+
+void validate_particle_type(ParticleType type, const std::string& context)
+{
+  if (type.is_transportable())
+    return;
+
+  fatal_error(
+    fmt::format("Unsupported source particle type '{}' (PDG {}) in {}.",
+      type.str(), type.pdg_number(), context));
+}
+
+} // namespace
+
 //==============================================================================
 // Global variables
 //==============================================================================
@@ -284,22 +298,15 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node)
 {
   // Check for particle type
   if (check_for_node(node, "particle")) {
-    auto temp_str = get_node_value(node, "particle", true, true);
-    if (temp_str == "neutron") {
-      particle_ = ParticleType::neutron;
-    } else if (temp_str == "photon") {
-      particle_ = ParticleType::photon;
+    auto temp_str = get_node_value(node, "particle", false, true);
+    particle_ = ParticleType(temp_str);
+    if (particle_ == ParticleType::photon() ||
+        particle_ == ParticleType::electron() ||
+        particle_ == ParticleType::positron()) {
       settings::photon_transport = true;
-    } else if (temp_str == "electron") {
-      particle_ = ParticleType::electron;
-      settings::photon_transport = true;
-    } else if (temp_str == "positron") {
-      particle_ = ParticleType::positron;
-      settings::photon_transport = true;
-    } else {
-      fatal_error(std::string("Unknown source particle type: ") + temp_str);
     }
   }
+  validate_particle_type(particle_, "IndependentSource");
 
   // Check for external source file
   if (check_for_node(node, "file")) {
@@ -390,7 +397,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
   // Sample energy and time for neutron and photon sources
   if (settings::solver_type != SolverType::RANDOM_RAY) {
     // Check for monoenergetic source above maximum particle energy
-    auto p = static_cast<int>(particle_);
+    auto p = particle_.transport_index();
     auto energy_ptr = dynamic_cast<Discrete*>(energy_.get());
     if (energy_ptr) {
       auto energies = xt::adapt(energy_ptr->x());
@@ -472,6 +479,11 @@ void FileSource::load_sites_from_file(const std::string& path)
     // Close file
     file_close(file_id);
   }
+
+  // Make sure particles in source file have valid types
+  for (const auto& site : this->sites_) {
+    validate_particle_type(site.particle, "FileSource");
+  }
 }
 
 SourceSite FileSource::sample(uint64_t* seed) const
@@ -585,6 +597,11 @@ MeshSource::MeshSource(pugi::xml_node node) : Source(node)
       std::make_unique<MeshElementSpatial>(mesh_idx, elem_index));
   }
 
+  // Make sure sources use valid particle types
+  for (const auto& src : sources_) {
+    validate_particle_type(src->particle_type(), "MeshSource");
+  }
+
   // the number of source distributions should either be one or equal to the
   // number of mesh elements
   if (sources_.size() > 1 && sources_.size() != mesh->n_bins()) {

src/state_point.cpp

@@ -22,6 +22,7 @@
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/output.h"
+#include "openmc/particle_type.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
 #include "openmc/tallies/derivative.h"
@@ -632,6 +633,7 @@ void write_h5_source_point(const char* filename, span<SourceSite> source_bank,
   if (mpi::master || parallel) {
     file_id = file_open(filename_.c_str(), 'w', true);
     write_attribute(file_id, "filetype", "source");
+    write_attribute(file_id, "version", VERSION_STATEPOINT);
   }
 
   // Get pointer to source bank and write to file
@@ -677,6 +679,16 @@ std::string dtype_member_names(hid_t dtype_id)
 void read_source_bank(
   hid_t group_id, vector<SourceSite>& sites, bool distribute)
 {
+  bool legacy_particle_codes = true;
+  if (attribute_exists(group_id, "version")) {
+    array<int, 2> version;
+    read_attribute(group_id, "version", version);
+    if (version[0] > VERSION_STATEPOINT[0] ||
+        (version[0] == VERSION_STATEPOINT[0] && version[1] >= 2)) {
+      legacy_particle_codes = false;
+    }
+  }
+
   hid_t banktype = h5banktype(true);
 
   // Open the dataset
@@ -738,6 +750,12 @@ void read_source_bank(
     H5Sclose(memspace);
   H5Dclose(dset);
   H5Tclose(banktype);
+
+  if (legacy_particle_codes) {
+    for (auto& site : sites) {
+      site.particle = legacy_particle_index_to_type(site.particle.pdg_number());
+    }
+  }
 }
 
 void write_unstructured_mesh_results()

src/tallies/filter_particle.cpp

@@ -13,7 +13,7 @@ void ParticleFilter::from_xml(pugi::xml_node node)
   // Convert to vector of ParticleType
   vector<ParticleType> types;
   for (auto& p : particles) {
-    types.push_back(str_to_particle_type(p));
+    types.emplace_back(p);
   }
   this->set_particles(types);
 }
@@ -47,7 +47,7 @@ void ParticleFilter::to_statepoint(hid_t filter_group) const
   Filter::to_statepoint(filter_group);
   vector<std::string> particles;
   for (auto p : particles_) {
-    particles.push_back(particle_type_to_str(p));
+    particles.push_back(p.str());
   }
   write_dataset(filter_group, "bins", particles);
 }
@@ -55,10 +55,10 @@ void ParticleFilter::to_statepoint(hid_t filter_group) const
 std::string ParticleFilter::text_label(int bin) const
 {
   const auto& p = particles_.at(bin);
-  return fmt::format("Particle: {}", particle_type_to_str(p));
+  return fmt::format("Particle: {}", p.str());
 }
 
-extern "C" int openmc_particle_filter_get_bins(int32_t idx, int bins[])
+extern "C" int openmc_particle_filter_get_bins(int32_t idx, int32_t bins[])
 {
   if (int err = verify_filter(idx))
     return err;
@@ -68,7 +68,7 @@ extern "C" int openmc_particle_filter_get_bins(int32_t idx, int bins[])
   if (pf) {
     const auto& particles = pf->particles();
     for (int i = 0; i < particles.size(); i++) {
-      bins[i] = static_cast<int>(particles[i]);
+      bins[i] = particles[i].pdg_number();
     }
   } else {
     set_errmsg("The filter at the specified index is not a ParticleFilter");

src/tallies/tally.cpp

@@ -231,12 +231,12 @@ Tally::Tally(pugi::xml_node node)
       const auto& f = model::tally_filters[particle_filter_index].get();
       auto pf = dynamic_cast<ParticleFilter*>(f);
       for (auto p : pf->particles()) {
-        if (p != ParticleType::neutron) {
+        if (!p.is_neutron()) {
           warning(fmt::format(
             "Particle filter other than NEUTRON used with "
             "photon transport turned off. All tallies for particle type {}"
             " will have no scores",
-            static_cast<int>(p)));
+            p.str()));
         }
       }
     }

src/tallies/tally_scoring.cpp

@@ -328,10 +328,10 @@ double score_neutron_heating(const Particle& p, const Tally& tally, double flux,
 //! Helper function to obtain reaction Q value for photons and charged particles
 double get_reaction_q_value(const Particle& p)
 {
-  if (p.type() == ParticleType::photon && p.event_mt() == PAIR_PROD) {
+  if (p.type().is_photon() && p.event_mt() == PAIR_PROD) {
     // pair production
     return -2 * MASS_ELECTRON_EV;
-  } else if (p.type() == ParticleType::positron) {
+  } else if (p.type() == ParticleType::positron()) {
     // positron annihilation
     return 2 * MASS_ELECTRON_EV;
   } else {
@@ -344,7 +344,7 @@ double get_reaction_q_value(const Particle& p)
 double score_particle_heating(const Particle& p, const Tally& tally,
   double flux, int rxn_bin, int i_nuclide, double atom_density)
 {
-  if (p.type() == ParticleType::neutron)
+  if (p.type().is_neutron())
     return score_neutron_heating(
       p, tally, flux, rxn_bin, i_nuclide, atom_density);
   if (i_nuclide == -1 || i_nuclide == p.event_nuclide() ||
@@ -584,8 +584,6 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
   // Get the pre-collision energy of the particle.
   auto E = p.E_last();
 
-  using Type = ParticleType;
-
   for (auto i = 0; i < tally.scores_.size(); ++i) {
     auto score_bin = tally.scores_[i];
     auto score_index = start_index + i;
@@ -598,9 +596,9 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
 
     case SCORE_TOTAL:
       if (i_nuclide >= 0) {
-        if (p.type() == Type::neutron) {
+        if (p.type().is_neutron()) {
           score = p.neutron_xs(i_nuclide).total * atom_density * flux;
-        } else if (p.type() == Type::photon) {
+        } else if (p.type().is_photon()) {
           score = p.photon_xs(i_nuclide).total * atom_density * flux;
         }
       } else {
@@ -609,7 +607,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_INVERSE_VELOCITY:
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Score inverse velocity in units of s/cm.
@@ -617,11 +615,11 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_SCATTER:
-      if (p.type() != Type::neutron && p.type() != Type::photon)
+      if (!p.type().is_neutron() && !p.type().is_photon())
         continue;
 
       if (i_nuclide >= 0) {
-        if (p.type() == Type::neutron) {
+        if (p.type().is_neutron()) {
           const auto& micro = p.neutron_xs(i_nuclide);
           score = (micro.total - micro.absorption) * atom_density * flux;
         } else {
@@ -629,7 +627,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
           score = (micro.coherent + micro.incoherent) * atom_density * flux;
         }
       } else {
-        if (p.type() == Type::neutron) {
+        if (p.type().is_neutron()) {
           score = (p.macro_xs().total - p.macro_xs().absorption) * flux;
         } else {
           score = (p.macro_xs().coherent + p.macro_xs().incoherent) * flux;
@@ -638,11 +636,11 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_ABSORPTION:
-      if (p.type() != Type::neutron && p.type() != Type::photon)
+      if (!p.type().is_neutron() && !p.type().is_photon())
         continue;
 
       if (i_nuclide >= 0) {
-        if (p.type() == Type::neutron) {
+        if (p.type().is_neutron()) {
           score = p.neutron_xs(i_nuclide).absorption * atom_density * flux;
         } else {
           const auto& xs = p.photon_xs(i_nuclide);
@@ -650,7 +648,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
             (xs.total - xs.coherent - xs.incoherent) * atom_density * flux;
         }
       } else {
-        if (p.type() == Type::neutron) {
+        if (p.type().is_neutron()) {
           score = p.macro_xs().absorption * flux;
         } else {
           score =
@@ -806,7 +804,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
           // this loop.
           for (auto d = 1; d < rxn.products_.size(); ++d) {
             const auto& product = rxn.products_[d];
-            if (product.particle_ != Type::neutron)
+            if (!product.particle_.is_neutron())
               continue;
 
             auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
@@ -860,7 +858,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
                 // this loop.
                 for (auto d = 1; d < rxn.products_.size(); ++d) {
                   const auto& product = rxn.products_[d];
-                  if (product.particle_ != Type::neutron)
+                  if (!product.particle_.is_neutron())
                     continue;
 
                   auto yield =
@@ -911,7 +909,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       continue;
 
     case ELASTIC:
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       if (i_nuclide >= 0) {
@@ -942,7 +940,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_IFP_TIME_NUM:
-      if ((p.type() == Type::neutron) && (p.fission())) {
+      if (p.type().is_neutron() && p.fission()) {
         const auto& lifetime =
           simulation::ifp_source_lifetime_bank[p.current_work() - 1][0];
         score = lifetime * p.wgt_last();
@@ -950,7 +948,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_IFP_BETA_NUM:
-      if ((p.type() == Type::neutron) && (p.fission())) {
+      if (p.type().is_neutron() && p.fission()) {
         const auto& delayed_group =
           simulation::ifp_source_delayed_group_bank[p.current_work() - 1][0];
         if (delayed_group > 0) {
@@ -968,7 +966,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_IFP_DENOM:
-      if ((p.type() == Type::neutron) && (p.fission()))
+      if (p.type().is_neutron() && p.fission())
         score = p.wgt_last();
       break;
 
@@ -984,7 +982,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       if (!simulation::need_depletion_rx)
         goto default_case;
 
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       int m;
@@ -1017,7 +1015,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
     case INCOHERENT:
     case PHOTOELECTRIC:
     case PAIR_PROD:
-      if (p.type() != Type::photon)
+      if (!p.type().is_photon())
         continue;
 
       if (i_nuclide >= 0) {
@@ -1047,7 +1045,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
 
       // The default block is really only meant for redundant neutron reactions
       // (e.g. 444, 901)
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Any other cross section has to be calculated on-the-fly
@@ -1101,8 +1099,6 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
                             p.neutron_xs(p.event_nuclide()).total
                         : 0.0;
 
-  using Type = ParticleType;
-
   for (auto i = 0; i < tally.scores_.size(); ++i) {
     auto score_bin = tally.scores_[i];
     auto score_index = start_index + i;
@@ -1113,7 +1109,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       // All events score to a flux bin. We actually use a collision estimator
       // in place of an analog one since there is no way to count 'events'
       // exactly for the flux
-      if (p.type() == Type::neutron || p.type() == Type::photon) {
+      if (p.type().is_neutron() || p.type().is_photon()) {
         score = flux * p.wgt_last() / p.macro_xs().total;
       } else {
         score = 0.;
@@ -1127,7 +1123,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_INVERSE_VELOCITY:
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // All events score to an inverse velocity bin. We actually use a
@@ -1137,7 +1133,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_SCATTER:
-      if (p.type() != Type::neutron && p.type() != Type::photon)
+      if (!p.type().is_neutron() && !p.type().is_photon())
         continue;
 
       // Skip any event where the particle didn't scatter
@@ -1149,7 +1145,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_NU_SCATTER:
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Only analog estimators are available.
@@ -1174,7 +1170,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_ABSORPTION:
-      if (p.type() != Type::neutron && p.type() != Type::photon)
+      if (!p.type().is_neutron() && !p.type().is_photon())
         continue;
 
       if (settings::survival_biasing) {
@@ -1403,7 +1399,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
             // this loop.
             for (auto d = 1; d < rxn.products_.size(); ++d) {
               const auto& product = rxn.products_[d];
-              if (product.particle_ != Type::neutron)
+              if (!product.particle_.is_neutron())
                 continue;
 
               auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
@@ -1495,7 +1491,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       continue;
 
     case ELASTIC:
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Check if event MT matches
@@ -1524,7 +1520,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       if (!simulation::need_depletion_rx)
         goto default_case;
 
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Check if the event MT matches
@@ -1537,7 +1533,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
     case INCOHERENT:
     case PHOTOELECTRIC:
     case PAIR_PROD:
-      if (p.type() != Type::photon)
+      if (!p.type().is_photon())
         continue;
 
       if (score_bin == PHOTOELECTRIC) {
@@ -1564,7 +1560,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
 
       // The default block is really only meant for redundant neutron reactions
       // (e.g. 444, 901)
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Any other score is assumed to be a MT number. Thus, we just need
@@ -2288,10 +2284,7 @@ void score_analog_tally_ce(Particle& p)
   // Since electrons/positrons are not transported, we assign a flux of zero.
   // Note that the heating score does NOT use the flux and will be non-zero for
   // electrons/positrons.
-  double flux =
-    (p.type() == ParticleType::neutron || p.type() == ParticleType::photon)
-      ? 1.0
-      : 0.0;
+  double flux = (p.type().is_neutron() || p.type().is_photon()) ? 1.0 : 0.0;
 
   for (auto i_tally : model::active_analog_tallies) {
     const Tally& tally {*model::tallies[i_tally]};
@@ -2419,7 +2412,7 @@ void score_tracklength_tally_general(
             if (j == C_NONE) {
               // Determine log union grid index
               if (i_log_union == C_NONE) {
-                int neutron = static_cast<int>(ParticleType::neutron);
+                int neutron = ParticleType::neutron().transport_index();
                 i_log_union = std::log(p.E() / data::energy_min[neutron]) /
                               simulation::log_spacing;
               }
@@ -2516,7 +2509,7 @@ void score_collision_tally(Particle& p)
 {
   // Determine the collision estimate of the flux
   double flux = 0.0;
-  if (p.type() == ParticleType::neutron || p.type() == ParticleType::photon) {
+  if (p.type().is_neutron() || p.type().is_photon()) {
     flux = p.wgt_last() / p.macro_xs().total;
   }
 
@@ -2550,7 +2543,7 @@ void score_collision_tally(Particle& p)
           if (j == C_NONE) {
             // Determine log union grid index
             if (i_log_union == C_NONE) {
-              int neutron = static_cast<int>(ParticleType::neutron);
+              int neutron = ParticleType::neutron().transport_index();
               i_log_union = std::log(p.E() / data::energy_min[neutron]) /
                             simulation::log_spacing;
             }

src/track_output.cpp

@@ -122,7 +122,7 @@ void finalize_particle_track(Particle& p)
   int offset = 0;
   for (auto& track_i : p.tracks()) {
     offsets.push_back(offset);
-    particles.push_back(static_cast<int>(track_i.particle));
+    particles.push_back(track_i.particle.pdg_number());
     offset += track_i.states.size();
     tracks.insert(tracks.end(), track_i.states.begin(), track_i.states.end());
   }

src/weight_windows.cpp

@@ -59,7 +59,7 @@ void apply_weight_windows(Particle& p)
     return;
 
   // WW on photon and neutron only
-  if (p.type() != ParticleType::neutron && p.type() != ParticleType::photon)
+  if (!p.type().is_neutron() && !p.type().is_photon())
     return;
 
   // skip dead or no energy
@@ -179,7 +179,7 @@ WeightWindows::WeightWindows(pugi::xml_node node)
 
   // get the particle type
   auto particle_type_str = std::string(get_node_value(node, "particle_type"));
-  particle_type_ = openmc::str_to_particle_type(particle_type_str);
+  particle_type_ = ParticleType {particle_type_str};
 
   // Determine associated mesh
   int32_t mesh_id = std::stoi(get_node_value(node, "mesh"));
@@ -252,7 +252,7 @@ WeightWindows* WeightWindows::from_hdf5(
 
   std::string particle_type;
   read_dataset(ww_group, "particle_type", particle_type);
-  wws->particle_type_ = openmc::str_to_particle_type(particle_type);
+  wws->particle_type_ = ParticleType {particle_type};
 
   read_dataset<double>(ww_group, "energy_bounds", wws->energy_bounds_);
 
@@ -284,7 +284,10 @@ void WeightWindows::set_defaults()
 {
   // set energy bounds to the min/max energy supported by the data
   if (energy_bounds_.size() == 0) {
-    int p_type = static_cast<int>(particle_type_);
+    int p_type = particle_type_.transport_index();
+    if (p_type == C_NONE) {
+      fatal_error("Weight windows particle is not supported for transport.");
+    }
     energy_bounds_.push_back(data::energy_min[p_type]);
     energy_bounds_.push_back(data::energy_max[p_type]);
   }
@@ -345,10 +348,9 @@ void WeightWindows::set_energy_bounds(span<const double> bounds)
 
 void WeightWindows::set_particle_type(ParticleType p_type)
 {
-  if (p_type != ParticleType::neutron && p_type != ParticleType::photon)
-    fatal_error(
-      fmt::format("Particle type '{}' cannot be applied to weight windows.",
-        particle_type_to_str(p_type)));
+  if (!p_type.is_neutron() && !p_type.is_photon())
+    fatal_error(fmt::format(
+      "Particle type '{}' cannot be applied to weight windows.", p_type.str()));
   particle_type_ = p_type;
 }
 
@@ -608,8 +610,7 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value,
     if (p_it == particles.end()) {
       auto msg = fmt::format("Particle type '{}' not present on Filter {} for "
                              "Tally {} used to update WeightWindows {}",
-        particle_type_to_str(this->particle_type_), pf->id(), tally->id(),
-        this->id());
+        this->particle_type_.str(), pf->id(), tally->id(), this->id());
       fatal_error(msg);
     }
 
@@ -818,8 +819,7 @@ void WeightWindows::to_hdf5(hid_t group) const
   hid_t ww_group = create_group(group, fmt::format("weight_windows_{}", id()));
 
   write_dataset(ww_group, "mesh", this->mesh()->id());
-  write_dataset(
-    ww_group, "particle_type", openmc::particle_type_to_str(particle_type_));
+  write_dataset(ww_group, "particle_type", particle_type_.str());
   write_dataset(ww_group, "energy_bounds", energy_bounds_);
   write_dataset(ww_group, "lower_ww_bounds", lower_ww_);
   write_dataset(ww_group, "upper_ww_bounds", upper_ww_);
@@ -846,8 +846,8 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
         max_realizations_, active_batches);
     warning(msg);
   }
-  auto tmp_str = get_node_value(node, "particle_type", true, true);
-  auto particle_type = str_to_particle_type(tmp_str);
+  auto tmp_str = get_node_value(node, "particle_type", false, true);
+  auto particle_type = ParticleType {tmp_str};
 
   update_interval_ = std::stoi(get_node_value(node, "update_interval"));
   on_the_fly_ = get_node_value_bool(node, "on_the_fly");
@@ -856,7 +856,10 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
   if (check_for_node(node, "energy_bounds")) {
     e_bounds = get_node_array<double>(node, "energy_bounds");
   } else {
-    int p_type = static_cast<int>(particle_type);
+    int p_type = particle_type.transport_index();
+    if (p_type == C_NONE) {
+      fatal_error("Weight windows particle is not supported for transport.");
+    }
     e_bounds.push_back(data::energy_min[p_type]);
     e_bounds.push_back(data::energy_max[p_type]);
   }
@@ -1110,23 +1113,25 @@ extern "C" int openmc_weight_windows_get_energy_bounds(
   return 0;
 }
 
-extern "C" int openmc_weight_windows_set_particle(int32_t index, int particle)
+extern "C" int openmc_weight_windows_set_particle(
+  int32_t index, int32_t particle)
 {
   if (int err = verify_ww_index(index))
     return err;
 
   const auto& wws = variance_reduction::weight_windows.at(index);
-  wws->set_particle_type(static_cast<ParticleType>(particle));
+  wws->set_particle_type(ParticleType {particle});
   return 0;
 }
 
-extern "C" int openmc_weight_windows_get_particle(int32_t index, int* particle)
+extern "C" int openmc_weight_windows_get_particle(
+  int32_t index, int32_t* particle)
 {
   if (int err = verify_ww_index(index))
     return err;
 
   const auto& wws = variance_reduction::weight_windows.at(index);
-  *particle = static_cast<int>(wws->particle_type());
+  *particle = wws->particle_type().pdg_number();
   return 0;
 }
 

tests/cpp_unit_tests/CMakeLists.txt

@@ -6,6 +6,7 @@ set(TEST_NAMES
   test_math
   test_mcpl_stat_sum
   test_mesh
+  test_region
   # Add additional unit test files here
 )
 

tests/cpp_unit_tests/test_distribution.cpp

@@ -1,4 +1,6 @@
 #include "openmc/distribution.h"
+#include "openmc/distribution_spatial.h"
+#include "openmc/position.h"
 #include "openmc/random_lcg.h"
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
@@ -79,3 +81,14 @@ TEST_CASE("Test alias sampling method for pugixml constructor")
     REQUIRE(dist.alias()[i] == correct_alias[i]);
   }
 }
+
+TEST_CASE("Test construction of SpatialBox with parameters")
+{
+  openmc::Position ll {-1, -2, -3};
+  openmc::Position ur {30, 15, 5};
+  openmc::SpatialBox box(ll, ur);
+
+  REQUIRE(box.lower_left() == openmc::Position {-1, -2, -3});
+  REQUIRE(box.upper_right() == openmc::Position {30, 15, 5});
+  REQUIRE_FALSE(box.only_fissionable());
+}

tests/cpp_unit_tests/test_mcpl_stat_sum.cpp

@@ -25,7 +25,7 @@ TEST_CASE("MCPL stat:sum field")
 
     // Initialize test particles
     for (int i = 0; i < 100; ++i) {
-      source_bank[i].particle = openmc::ParticleType::neutron;
+      source_bank[i].particle = openmc::ParticleType::neutron();
       source_bank[i].r = {i * 0.1, i * 0.2, i * 0.3};
       source_bank[i].u = {0.0, 0.0, 1.0};
       source_bank[i].E = 2.0e6; // 2 MeV
@@ -63,7 +63,7 @@ TEST_CASE("MCPL stat:sum field")
 
       // Initialize particles
       for (int i = 0; i < count; ++i) {
-        source_bank[i].particle = openmc::ParticleType::neutron;
+        source_bank[i].particle = openmc::ParticleType::neutron();
         source_bank[i].r = {0.0, 0.0, 0.0};
         source_bank[i].u = {0.0, 0.0, 1.0};
         source_bank[i].E = 1.0e6;

tests/cpp_unit_tests/test_region.cpp

@@ -0,0 +1,101 @@
+#include <catch2/catch_test_macros.hpp>
+
+#include "openmc/cell.h"
+#include "openmc/surface.h"
+
+#include <pugixml.hpp>
+
+namespace {
+
+// Helper class to set up and tear down test surfaces
+class SurfaceFixture {
+public:
+  SurfaceFixture()
+  {
+    pugi::xml_document doc;
+    pugi::xml_node surf_node = doc.append_child("surface");
+    surf_node.set_name("surface");
+    surf_node.append_attribute("id") = "0";
+    surf_node.append_attribute("type") = "x-plane";
+    surf_node.append_attribute("coeffs") = "1";
+
+    for (int i = 1; i < 10; ++i) {
+      surf_node.attribute("id") = i;
+      openmc::model::surfaces.push_back(
+        std::make_unique<openmc::SurfaceXPlane>(surf_node));
+      openmc::model::surface_map[i] = i - 1;
+    }
+  }
+
+  ~SurfaceFixture()
+  {
+    openmc::model::surfaces.clear();
+    openmc::model::surface_map.clear();
+  }
+};
+
+} // anonymous namespace
+
+TEST_CASE("Test region simplification")
+{
+  SurfaceFixture fixture;
+
+  SECTION("Original bug case from issue #3685")
+  {
+    // Input: "-1 2 (-3 4) | (-5 6)" was being incorrectly interpreted
+    auto region = openmc::Region("(-1 2 (-3 4) | (-5 6))", 0);
+    REQUIRE(region.str() == " ( ( -1 2 ( -3 4 ) ) | ( -5 6 ) )");
+  }
+
+  SECTION("Simple union - no extra parentheses needed")
+  {
+    auto region = openmc::Region("1 | 2", 0);
+    REQUIRE(region.str() == " 1 | 2");
+  }
+
+  SECTION("Intersection then union")
+  {
+    // Intersection should have higher precedence, so (1 2) grouped
+    auto region = openmc::Region("1 2 | 3", 0);
+    REQUIRE(region.str() == " ( 1 2 ) | 3");
+  }
+
+  SECTION("Union then intersection")
+  {
+    // The (2 3) intersection should be grouped
+    auto region = openmc::Region("1 | 2 3", 0);
+    REQUIRE(region.str() == " 1 | ( 2 3 )");
+  }
+
+  SECTION("Nested parentheses preserved")
+  {
+    // These parentheses are meaningful and should be preserved
+    auto region = openmc::Region("(1 | 2) (3 | 4)", 0);
+    REQUIRE(region.str() == " ( 1 | 2 ) ( 3 | 4 )");
+  }
+
+  SECTION("Deep nesting")
+  {
+    auto region = openmc::Region("((1 2) | (3 4)) 5", 0);
+    REQUIRE(region.str() == " ( ( 1 2 ) | ( 3 4 ) ) 5");
+  }
+
+  SECTION("Multiple unions")
+  {
+    auto region = openmc::Region("1 | 2 | 3", 0);
+    REQUIRE(region.str() == " 1 | 2 | 3");
+  }
+
+  SECTION("Multiple intersections")
+  {
+    auto region = openmc::Region("1 2 3", 0);
+    // Simple cell - no operators in output
+    REQUIRE(region.str() == " 1 2 3");
+  }
+
+  SECTION("Complex mixed expression")
+  {
+    auto region = openmc::Region("1 2 | 3 4 | 5 6", 0);
+    REQUIRE(region.str() == " ( 1 2 ) | ( 3 4 ) | ( 5 6 )");
+  }
+}

tests/dummy_operator.py

@@ -243,4 +243,4 @@ class DummyOperator(TransportOperator):
             Maps cell name to index in global geometry.
 
         """
-        return self.volume, self.nuc_list, self.local_mats, self.burnable_mats
+        return self.volume, self.nuc_list, self.local_mats, self.burnable_mats, {"1": ""}

tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat

@@ -0,0 +1,73 @@
+<?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>
+  <tallies>
+    <filter id="61" type="material">
+      <bins>1</bins>
+    </filter>
+    <tally id="193" name="KF Tally">
+      <filters>61</filters>
+      <scores>kappa-fission</scores>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/results_true.dat

@@ -0,0 +1,5 @@
+k-combined:
+7.479770E-01 1.624548E-02
+tally 1:
+2.965503E+08
+1.762157E+16

tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat

@@ -0,0 +1,73 @@
+<?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>
+  <tallies>
+    <filter id="97" type="material">
+      <bins>1</bins>
+    </filter>
+    <tally id="203" name="KF Tally">
+      <filters>97</filters>
+      <scores>kappa-fission</scores>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/results_true.dat

@@ -0,0 +1,5 @@
+k-combined:
+7.372542E-01 6.967831E-03
+tally 1:
+2.909255E+08
+1.693395E+16

tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat

@@ -0,0 +1,73 @@
+<?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>
+  <tallies>
+    <filter id="97" type="material">
+      <bins>1</bins>
+    </filter>
+    <tally id="203" name="KF Tally">
+      <filters>97</filters>
+      <scores>kappa-fission</scores>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/results_true.dat

@@ -0,0 +1,5 @@
+k-combined:
+6.413334E-01 2.083132E-02
+tally 1:
+2.541463E+08
+1.297259E+16

tests/regression_tests/random_ray_auto_convert_kappa_fission/test.py

@@ -0,0 +1,60 @@
+import os
+
+import openmc
+from openmc.examples import pwr_pin_cell
+from openmc import RegularMesh
+from openmc.utility_funcs import change_directory
+import pytest
+
+from tests.testing_harness import TolerantPyAPITestHarness
+
+
+class MGXSTestHarness(TolerantPyAPITestHarness):
+    def _cleanup(self):
+        super()._cleanup()
+        f = 'mgxs.h5'
+        if os.path.exists(f):
+            os.remove(f)
+
+
+@pytest.mark.parametrize("method", ["material_wise", "stochastic_slab", "infinite_medium"])
+def test_random_ray_auto_convert(method):
+    with change_directory(method):
+        openmc.reset_auto_ids()
+
+        # Start with a normal continuous energy model
+        model = pwr_pin_cell()
+
+        # Convert to a multi-group model
+        model.convert_to_multigroup(
+            method=method, groups='CASMO-2', nparticles=100,
+            overwrite_mgxs_library=False, mgxs_path="mgxs.h5"
+        )
+
+        # Convert to a random ray model
+        model.convert_to_random_ray()
+
+        # Set the number of particles
+        model.settings.particles = 100
+
+        # Overlay a basic 2x2 mesh
+        n = 2
+        mesh = RegularMesh()
+        mesh.dimension = (n, n)
+        bbox = model.geometry.bounding_box
+        mesh.lower_left = (bbox.lower_left[0], bbox.lower_left[1])
+        mesh.upper_right = (bbox.upper_right[0], bbox.upper_right[1])
+        model.settings.random_ray['source_region_meshes'] = [
+            (mesh, [model.geometry.root_universe])]
+
+        # Set the source shape to linear
+        model.settings.random_ray['source_shape'] = 'linear'
+
+        # Set a material tally
+        t = openmc.Tally(name = 'KF Tally')
+        t.filters = [openmc.MaterialFilter(bins=1)]
+        t.scores = ['kappa-fission']
+        model.tallies.append(t)
+
+        harness = MGXSTestHarness('statepoint.10.h5', model)
+        harness.main()

tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat

@@ -0,0 +1,109 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="UO2 fuel">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="UO2"/>
+    </material>
+    <material id="2" name="Water">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="LWTR"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" name="fuel inner a" material="1" region="-2" density="2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0" universe="1"/>
+    <cell id="2" name="fuel inner b" material="1" region="2 -3" density="2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0" universe="1"/>
+    <cell id="3" name="fuel inner c" material="1" region="3 -1" density="2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0" universe="1"/>
+    <cell id="4" name="moderator inner a" material="2" region="1 -4" universe="1"/>
+    <cell id="5" name="moderator outer b" material="2" region="4 -5" universe="1"/>
+    <cell id="6" name="moderator outer c" material="2" region="5" universe="1"/>
+    <cell id="7" name="azimuthal_cell_0" fill="1" region="6 -7" universe="2"/>
+    <cell id="8" name="azimuthal_cell_1" fill="1" region="7 -8" universe="2"/>
+    <cell id="9" name="azimuthal_cell_2" fill="1" region="8 -9" universe="2"/>
+    <cell id="10" name="azimuthal_cell_3" fill="1" region="9 -10" universe="2"/>
+    <cell id="11" name="azimuthal_cell_4" fill="1" region="10 -11" universe="2"/>
+    <cell id="12" name="azimuthal_cell_5" fill="1" region="11 -12" universe="2"/>
+    <cell id="13" name="azimuthal_cell_6" fill="1" region="12 -13" universe="2"/>
+    <cell id="14" name="azimuthal_cell_7" fill="1" region="13 -6" universe="2"/>
+    <cell id="15" name="moderator infinite" material="2" universe="3"/>
+    <cell id="16" fill="4" universe="5"/>
+    <cell id="17" name="assembly" fill="6" region="14 -15 16 -17" universe="7"/>
+    <lattice id="4">
+      <pitch>0.126 0.126</pitch>
+      <dimension>10 10</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <universes>
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 </universes>
+    </lattice>
+    <lattice id="6">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+2 2 
+2 5 </universes>
+    </lattice>
+    <surface id="1" name="Fuel OR" type="z-cylinder" coeffs="0.0 0.0 0.54"/>
+    <surface id="2" name="inner ring a" type="z-cylinder" coeffs="0.0 0.0 0.33"/>
+    <surface id="3" name="inner ring b" type="z-cylinder" coeffs="0.0 0.0 0.45"/>
+    <surface id="4" name="outer ring a" type="z-cylinder" coeffs="0.0 0.0 0.6"/>
+    <surface id="5" name="outer ring b" type="z-cylinder" coeffs="0.0 0.0 0.69"/>
+    <surface id="6" type="plane" coeffs="-0.0 1.0 0 0"/>
+    <surface id="7" type="plane" coeffs="-0.7071067811865475 0.7071067811865476 0 0"/>
+    <surface id="8" type="plane" coeffs="-1.0 6.123233995736766e-17 0 0"/>
+    <surface id="9" type="plane" coeffs="-0.7071067811865476 -0.7071067811865475 0 0"/>
+    <surface id="10" type="plane" coeffs="-1.2246467991473532e-16 -1.0 0 0"/>
+    <surface id="11" type="plane" coeffs="0.7071067811865475 -0.7071067811865477 0 0"/>
+    <surface id="12" type="plane" coeffs="1.0 -1.8369701987210297e-16 0 0"/>
+    <surface id="13" type="plane" coeffs="0.7071067811865477 0.7071067811865474 0 0"/>
+    <surface id="14" name="minimum x" type="x-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="15" name="maximum x" type="x-plane" boundary="reflective" coeffs="1.26"/>
+    <surface id="16" name="minimum y" type="y-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="17" name="maximum y" type="y-plane" boundary="reflective" coeffs="1.26"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>100.0</distance_active>
+      <distance_inactive>20.0</distance_inactive>
+      <source type="independent" strength="1.0" particle="neutron">
+        <space type="box">
+          <parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
+        </space>
+      </source>
+      <volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
+    </random_ray>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="energy">
+      <bins>1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0</bins>
+    </filter>
+    <tally id="1" name="Mesh tally">
+      <filters>1 2</filters>
+      <scores>flux fission nu-fission</scores>
+      <estimator>analog</estimator>
+    </tally>
+  </tallies>
+</model>