Approximate multigroup velocity

include/openmc/mgxs_interface.h

@@ -49,6 +49,9 @@ class MgxsInterface {
   // Get the group index corresponding to a continuous energy
   int get_group_index(double E);
 
+  // Get the inverse velocity corresponding to particle p
+  double get_inverse_velocity(const Particle& p);
+
   int num_energy_groups_;
   int num_delayed_groups_;
   vector<std::string> xs_names_;            // available names in HDF5 file
@@ -61,6 +64,8 @@ class MgxsInterface {
   vector<double> energy_bin_avg_;
   vector<double> rev_energy_bins_;
   vector<vector<double>> nuc_temps_; // all available temperatures
+  vector<double>
+    default_inverse_velocity_; // approximate default inverse-velocity data
 };
 
 namespace data {

include/openmc/particle.h

@@ -39,6 +39,8 @@ class Particle : public ParticleData {
 
   double speed() const;
 
+  double mass() const;
+
   //! create a secondary particle
   //
   //! stores the current phase space attributes of the particle in the

openmc/mgxs/groups.py

@@ -78,6 +78,7 @@ def group_edges(self):
     @group_edges.setter
     def group_edges(self, edges):
         cv.check_type('group edges', edges, Iterable, Real)
+        cv.check_increasing('group edges', edges)
         cv.check_greater_than('number of group edges', len(edges), 1)
         self._group_edges = np.array(edges)
 

src/mgxs_interface.cpp

@@ -193,6 +193,23 @@ int MgxsInterface::get_group_index(double E)
 
 //==============================================================================
 
+double MgxsInterface::get_inverse_velocity(const Particle& p)
+{
+  auto mat = p.material();
+  if (mat != MATERIAL_VOID) {
+    auto& macro_xs = macro_xs_[mat];
+    int macro_t = p.mg_xs_cache().t;
+    int macro_a = macro_xs.get_angle_index(p.u());
+    double inv_v =
+      macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p.g(), macro_t, macro_a);
+    if (inv_v > 0.0)
+      return inv_v;
+  }
+  return default_inverse_velocity_[p.g()];
+}
+
+//==============================================================================
+
 void MgxsInterface::read_header(const std::string& path_cross_sections)
 {
   // Save name of HDF5 file to be read to struct data
@@ -237,6 +254,18 @@ void MgxsInterface::read_header(const std::string& path_cross_sections)
                 "library file!");
   }
 
+  // Calculate approximate default inverse velocity data
+  for (int i = 0; i < energy_bins_.size() - 1; ++i) {
+    double e_min = std::max(energy_bins_[i], 1e-5);
+    double e_max = energy_bins_[i + 1];
+    double inv_v = (std::acosh(1 + e_max / MASS_NEUTRON_EV) -
+                     std::sqrt(1 + 2 * MASS_NEUTRON_EV / e_max) -
+                     std::acosh(1 + e_min / MASS_NEUTRON_EV) +
+                     std::sqrt(1 + 2 * MASS_NEUTRON_EV / e_min)) /
+                   (C_LIGHT * std::log(e_max / e_min));
+    default_inverse_velocity_.push_back(inv_v);
+  }
+
   // Close MGXS HDF5 file
   file_close(file_id);
 }

src/particle.cpp

@@ -47,26 +47,26 @@ double Particle::speed() const
 {
   if (settings::run_CE) {
     // Determine mass in eV/c^2
-    double mass;
-    switch (type().pdg_number()) {
-    case PDG_NEUTRON:
-      mass = MASS_NEUTRON_EV;
-    case PDG_ELECTRON:
-    case PDG_POSITRON:
-      mass = MASS_ELECTRON_EV;
-    default:
-      mass = this->type().mass() * AMU_EV;
-    }
+    double mass = this->mass();
 
     // 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);
   } else {
-    auto& macro_xs = data::mg.macro_xs_[this->material()];
-    int macro_t = this->mg_xs_cache().t;
-    int macro_a = macro_xs.get_angle_index(this->u());
-    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, this->g(), nullptr,
-                   nullptr, nullptr, macro_t, macro_a);
+    return 1.0 / data::mg.get_inverse_velocity(*this);
+  }
+}
+
+double Particle::mass() const
+{
+  switch (type().pdg_number()) {
+  case PDG_NEUTRON:
+    return MASS_NEUTRON_EV;
+  case PDG_ELECTRON:
+  case PDG_POSITRON:
+    return MASS_ELECTRON_EV;
+  default:
+    return this->type().mass() * AMU_EV;
   }
 }
 

tests/regression_tests/mg_basic_delayed/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-1.017078E+00 1.181139E-02
+1.012441E+00 7.095889E-03