Refactor ParticleType to use PDG Monte Carlo numbering scheme (#3756)

Co-authored-by: GuySten <62616591+GuySten@users.noreply.github.com>
Co-authored-by: Amanda Lund <alund1187@gmail.com>

include/openmc/capi.h

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