Approximate multigroup velocity

include/openmc/mgxs.h

@@ -95,10 +95,10 @@ class Mgxs {
   //!
   //! @param xs_id HDF5 group id for the cross section data.
   //! @param temperature Temperatures to read.
-  //! @param num_group number of energy groups
+  //! @param energy_bins energy bins
   //! @param num_delay number of delayed groups
-  Mgxs(hid_t xs_id, const vector<double>& temperature, int num_group,
-    int num_delay);
+  Mgxs(hid_t xs_id, const vector<double>& temperature,
+    const vector<double>& energy_bins, int num_delay);
 
   //! \brief Constructor that initializes and populates all data to build a
   //!   macroscopic cross section from microscopic cross sections.

include/openmc/mgxs_interface.h

@@ -61,6 +61,7 @@ class MgxsInterface {
   vector<double> energy_bin_avg_;
   vector<double> rev_energy_bins_;
   vector<vector<double>> nuc_temps_; // all available temperatures
+  vector<double> void_velocities_;   // velocity of particles in void regions
 };
 
 namespace data {

include/openmc/xsdata.h

@@ -120,10 +120,11 @@ class XsData {
   //!   the incoming particle.
   //! @param n_pol Number of polar angles.
   //! @param n_azi Number of azimuthal angles.
+  //! @param energy_bins energy bins of XsData
   void from_hdf5(hid_t xsdata_grp, bool fissionable,
     AngleDistributionType scatter_format,
     AngleDistributionType final_scatter_format, int order_data,
-    bool is_isotropic, int n_pol, int n_azi);
+    bool is_isotropic, int n_pol, int n_azi, const vector<double>& energy_bins);
 
   //! \brief Combines the microscopic data to a macroscopic object.
   //!

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

@@ -271,9 +271,9 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector<double>& temperature,
 
 //==============================================================================
 
-Mgxs::Mgxs(
-  hid_t xs_id, const vector<double>& temperature, int num_group, int num_delay)
-  : num_groups(num_group), num_delayed_groups(num_delay)
+Mgxs::Mgxs(hid_t xs_id, const vector<double>& temperature,
+  const vector<double>& energy_bins, int num_delay)
+  : num_groups(energy_bins.size() - 1), num_delayed_groups(num_delay)
 {
   // Call generic data gathering routine (will populate the metadata)
   int order_data;
@@ -296,7 +296,8 @@ Mgxs::Mgxs(
     hid_t xsdata_grp = open_group(xs_id, temp_str.c_str());
 
     xs[t].from_hdf5(xsdata_grp, fissionable, scatter_format,
-      final_scatter_format, order_data, is_isotropic, n_pol, n_azi);
+      final_scatter_format, order_data, is_isotropic, n_pol, n_azi,
+      energy_bins);
     close_group(xsdata_grp);
 
   } // end temperature loop

src/mgxs_interface.cpp

@@ -109,7 +109,7 @@ void MgxsInterface::add_mgxs(
   }
 
   nuclides_.emplace_back(
-    xs_grp, temperature, num_energy_groups_, num_delayed_groups_);
+    xs_grp, temperature, energy_bins_, num_delayed_groups_);
   close_group(xs_grp);
 }
 
@@ -237,6 +237,48 @@ void MgxsInterface::read_header(const std::string& path_cross_sections)
                 "library file!");
   }
 
+  // Read void velocities
+  bool void_velocities_exist = false;
+  if (object_exists(file_id, "void")) {
+    auto void_grp = open_group(file_id, "void");
+    // Determine the available temperatures
+    hid_t kT_group = open_group(void_grp, "kTs");
+    size_t num_temps = get_num_datasets(kT_group);
+    char** dset_names = new char*[num_temps];
+    for (int i = 0; i < num_temps; i++) {
+      dset_names[i] = new char[151];
+    }
+    get_datasets(kT_group, dset_names);
+    close_group(kT_group);
+    if (num_temps > 1)
+      fatal_error("Multigroup void data must have only one dummy temperature");
+    auto xsdata_grp = open_group(void_grp, dset_names[0]);
+    if (object_exists(xsdata_grp, "inverse-velocity")) {
+      vector<size_t> shape {1, static_cast<size_t>(num_energy_groups_)};
+      xt::xtensor<double, 2> inverse_velocity = xt::zeros<double>(shape);
+      read_nd_vector(xsdata_grp, "inverse-velocity", inverse_velocity);
+      auto velocity = 1.0 / inverse_velocity;
+      for (double v : velocity) {
+        void_velocities_.push_back(v);
+      }
+      void_velocities_exist = true;
+    }
+    close_group(xsdata_grp);
+    close_group(void_grp);
+  }
+  if (!void_velocities_exist) {
+    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 v = C_LIGHT * std::log(e_max / e_min) /
+                 (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));
+      void_velocities_.push_back(v);
+    }
+  }
+
   // Close MGXS HDF5 file
   file_close(file_id);
 }

src/particle.cpp

@@ -66,11 +66,14 @@ double Particle::speed() const
     return C_LIGHT * std::sqrt(this->E() * (this->E() + 2 * mass)) /
            (this->E() + mass);
   } else {
-    auto& macro_xs = data::mg.macro_xs_[this->material()];
+    auto mat = this->material();
+    if (mat == MATERIAL_VOID)
+      return data::mg.void_velocities_[this->g()];
+    auto& macro_xs = data::mg.macro_xs_[mat];
     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 / macro_xs.get_xs(
+                   MgxsType::INVERSE_VELOCITY, this->g(), macro_t, macro_a);
   }
 }
 

src/xsdata.cpp

@@ -81,7 +81,7 @@ XsData::XsData(bool fissionable, AngleDistributionType scatter_format,
 void XsData::from_hdf5(hid_t xsdata_grp, bool fissionable,
   AngleDistributionType scatter_format,
   AngleDistributionType final_scatter_format, int order_data, bool is_isotropic,
-  int n_pol, int n_azi)
+  int n_pol, int n_azi, const vector<double>& energy_bins)
 {
   // Reconstruct the dimension information so it doesn't need to be passed
   size_t n_ang = n_pol * n_azi;
@@ -96,6 +96,22 @@ void XsData::from_hdf5(hid_t xsdata_grp, bool fissionable,
   read_nd_vector(xsdata_grp, "absorption", absorption, true);
   read_nd_vector(xsdata_grp, "inverse-velocity", inverse_velocity);
 
+  if (!object_exists(xsdata_grp, "inverse-velocity")) {
+    xt::xarray<double> inv_vel = xt::zeros<double>({energy_groups});
+    for (int i = 0; i < energy_groups; ++i) {
+      double e_min = std::max(energy_bins[i], 1e-5);
+      double e_max = energy_bins[i + 1];
+      inv_vel[i] = (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)) /
+                   std::log(e_max / e_min) / C_LIGHT;
+    }
+    for (int j = 0; j < n_ang; ++j) {
+      xt::view(inverse_velocity, j, xt::all()) = inv_vel;
+    }
+  }
+
   // Get scattering data
   scatter_from_hdf5(
     xsdata_grp, n_ang, scatter_format, final_scatter_format, order_data);

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