Store atomic mass in ParticleType. (#3765)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

CMakeLists.txt

@@ -332,6 +332,7 @@ endif()
 #===============================================================================
 
 list(APPEND libopenmc_SOURCES
+  src/atomic_mass.cpp
   src/bank.cpp
   src/boundary_condition.cpp
   src/bremsstrahlung.cpp

include/openmc/atomic_mass.h

@@ -0,0 +1,28 @@
+//==============================================================================
+// atomic masses definitions
+//==============================================================================
+
+#ifndef OPENMC_ATOMIC_MASS_H
+#define OPENMC_ATOMIC_MASS_H
+
+#include <cstdint>
+#include <unordered_map>
+
+namespace openmc {
+
+// Values here are from the Committee on Data for Science and Technology
+// (CODATA) 2018 recommendation (https://physics.nist.gov/cuu/Constants/).
+
+// Physical constants
+constexpr double MASS_ELECTRON {5.48579909065e-4}; // mass of an electron in amu
+constexpr double MASS_NEUTRON {1.00866491595};     // mass of a neutron in amu
+constexpr double MASS_PROTON {1.007276466621};     // mass of a proton in amu
+constexpr double MASS_DEUTRON {2.013553212745};    // mass of a deutron in amu
+constexpr double MASS_HELION {3.014932247175};     // mass of a helion in amu
+constexpr double MASS_ALPHA {4.001506179127};      // mass of an alpha in amu
+
+extern std::unordered_map<int32_t, double> ATOMIC_MASS;
+
+} // namespace openmc
+
+#endif //  OPENMC_ATOMIC_MASS_H

include/openmc/constants.h

@@ -9,6 +9,7 @@
 #include <limits>
 
 #include "openmc/array.h"
+#include "openmc/atomic_mass.h"
 #include "openmc/vector.h"
 #include "openmc/version.h"
 
@@ -86,10 +87,10 @@ constexpr double INFTY {std::numeric_limits<double>::max()};
 // (CODATA) 2018 recommendation (https://physics.nist.gov/cuu/Constants/).
 
 // Physical constants
-constexpr double MASS_NEUTRON {1.00866491595}; // mass of a neutron in amu
+constexpr double AMU_EV {
+  9.3149410242e8}; // atomic mass unit energy equivalent in eV/c^2
 constexpr double MASS_NEUTRON_EV {
-  939.56542052e6};                             // mass of a neutron in eV/c^2
-constexpr double MASS_PROTON {1.007276466621}; // mass of a proton in amu
+  939.56542052e6}; // neutron mass energy equivalent in eV/c^2
 constexpr double MASS_ELECTRON_EV {
   0.51099895000e6}; // electron mass energy equivalent in eV/c^2
 constexpr double FINE_STRUCTURE {

include/openmc/particle_type.h

@@ -12,6 +12,7 @@
 #include <type_traits>
 
 #include "openmc/constants.h"
+#include "openmc/error.h"
 
 namespace openmc {
 
@@ -61,6 +62,17 @@ public:
   //----------------------------------------------------------------------------
   // Methods
 
+  // Get particle mass in [u]
+  double mass() const
+  {
+    int32_t p = std::abs(pdg_number_);
+    if (ATOMIC_MASS.count(p)) {
+      return ATOMIC_MASS[p];
+    } else {
+      fatal_error("Unknown mass for particle " + str());
+    }
+  }
+
   // Convert to string representation
   std::string str() const;
 

src/particle.cpp

@@ -48,20 +48,16 @@ double Particle::speed() const
   if (settings::run_CE) {
     // Determine mass in eV/c^2
     double mass;
-    switch (this->type().pdg_number()) {
+    switch (type().pdg_number()) {
     case PDG_NEUTRON:
       mass = MASS_NEUTRON_EV;
-      break;
-    case PDG_PHOTON:
-      mass = 0.0;
-      break;
     case PDG_ELECTRON:
     case PDG_POSITRON:
       mass = MASS_ELECTRON_EV;
-      break;
     default:
-      fatal_error("Unsupported particle for speed calculation.");
+      mass = this->type().mass() * AMU_EV;
     }
+
     // 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)) /
            (this->E() + mass);