Approximate multigroup velocity

docs/source/io_formats/mgxs_library.rst

@@ -133,6 +133,10 @@ Temperature-dependent data, provided for temperature <TTT>K.
              This dataset is optional. This is a 1-D vector if `representation`
              is "isotropic", or a 3-D vector if `representation` is "angle"
              with dimensions of [polar][azimuthal][groups].
+	     When this data is not available, an approximation using the
+	     group energy boundaries is used. For more information see
+	     the particle speed subsection in the multigroup-data section
+	     of the theory manual.
 
 **/<library name>/<TTT>K/scatter_data/**
 

docs/source/methods/cross_sections.rst

@@ -289,6 +289,37 @@ sections. This allows flexibility for the model to use highly anisotropic
 scattering information in the water while the fuel can be simulated with linear
 or even isotropic scattering.
 
+^^^^^^^^^^^^^^
+Particle Speed
+^^^^^^^^^^^^^^
+
+When using a multigroup representation of cross sections the particle speed
+has meaning only in an average sense. The particle speed is important when modeling dynamic behavior. OpenMC calculates the particle speed using the inverse velocity multigroup data if it is available. If it is not OpenMC uses an approximate velocity using the group energy bounds in the following way:
+
+.. math::
+
+   \frac{1}{v_g} = \int_{E_{\text{min}}^g}^{E_{\text{max}}^g} \frac{1}{v(E)} \frac{k}{E} dE 
+
+Where :math:`E_{\text{min}}^g` and :math:`E_{\text{max}}^g` are the group energy boundaries for group :math:`g`.
+:math:`v(E)` is the neutron velocity calculated using relativistic kinematics,
+:math:`k` is a normalization constant for the :math:`\frac{1}{E}` spectrum.
+
+This equation is valid when inside the group boundaries the neutron spectrum follows a typical :math:`\frac{1}{E}` slowing down spectrum. This assumption is widely used when generating fine group neutron cross section data libraries from continuous energy data.
+
+The solution to this equation is:
+
+.. math::
+
+   v_g = \frac{1.0}{c * log(\frac{E_{\text{max}}^g}{E_{\text{min}}^g})} * ( 2*(atan(k_max) - atan(k_min)) + (k_max-k_min))
+
+Where :math:`k_{\text{max}}`, :math:`k_{\text{min}}` are defined by a change of variables:
+
+.. math::
+
+   k = \sqrt{1+\frac{2 M_n}{E}}
+
+Where :math:`E` is the particle kinetic energy in eV and :math:`M_n` is the neutron mass in units of eV.
+
 .. _logarithmic mapping technique:
    https://mcnp.lanl.gov/pdf_files/TechReport_2014_LANL_LA-UR-14-24530_Brown.pdf
 .. _Hwang: https://doi.org/10.13182/NSE87-A16381

include/openmc/mgxs_interface.h

@@ -61,6 +61,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.cpp

@@ -510,6 +510,8 @@ double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
     break;
   case MgxsType::INVERSE_VELOCITY:
     val = xs_t->inverse_velocity(a, gin);
+    if (!(val > 0))
+      val = data::mg.default_inverse_velocity_[gin];
     break;
   case MgxsType::DECAY_RATE:
     if (dg != nullptr) {

src/mgxs_interface.cpp

@@ -237,6 +237,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 + 1], 1e-5);
+    double e_max = energy_bins_[i];
+    double f = 1.0 / (C_LIGHT * std::log(e_max / e_min));
+    double k_max = std::sqrt(1 + 2.0 * MASS_NEUTRON_EV / e_max);
+    double k_min = std::sqrt(1 + 2.0 * MASS_NEUTRON_EV / e_min);
+    double inv_v =
+      f * (2.0 * (std::atan(k_max) - std::atan(k_min)) - (k_max - k_min));
+    default_inverse_velocity_.push_back(inv_v);
+  }
+
   // Close MGXS HDF5 file
   file_close(file_id);
 }

src/particle.cpp

@@ -47,26 +47,33 @@ 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()];
+    auto mat = this->material();
+    if (mat == MATERIAL_VOID)
+      return 1.0 / data::mg.default_inverse_velocity_[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);
+  }
+}
+
+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/unit_tests/test_mg_inverse_velocity.py

@@ -0,0 +1,59 @@
+import openmc
+import numpy as np
+import pytest
+
+@pytest.fixture
+def one_group_lib():
+    groups = openmc.mgxs.EnergyGroups([0.0, 20.0e6])
+    xsdata = openmc.XSdata('slab_mat', groups)
+    xsdata.order = 0
+    xsdata.set_total([0.0])               
+    xsdata.set_absorption([0.0])
+    xsdata.set_scatter_matrix([[[0.0]]])
+
+    mg_library = openmc.MGXSLibrary(groups)
+    mg_library.add_xsdata(xsdata)
+    name = 'mgxs.h5'
+    mg_library.export_to_hdf5(name)
+    yield name
+
+@pytest.fixture
+def slab_model(one_group_lib):
+    model = openmc.Model()
+    mat = openmc.Material(name='slab_material')
+    mat.set_density('macro', 1.0)
+    mat.add_macroscopic('slab_mat')
+
+    model.materials = openmc.Materials([mat])
+    model.materials.cross_sections = one_group_lib
+
+    x_min = openmc.XPlane(x0=0.0, boundary_type='vacuum')
+    x_max = openmc.XPlane(x0=10.0, boundary_type='vacuum')
+
+    y_min = openmc.YPlane(y0=-10.0, boundary_type='vacuum')
+    y_max = openmc.YPlane(y0=10.0, boundary_type='vacuum')
+    z_min = openmc.ZPlane(z0=-10.0, boundary_type='vacuum')
+    z_max = openmc.ZPlane(z0=19.0, boundary_type='vacuum')
+
+    cell = openmc.Cell(fill=mat, region=+z_min & -x_max & +y_min & -y_max & +z_min & -z_max)
+    model.geometry = openmc.Geometry([cell])
+
+    model.settings = openmc.Settings()
+    model.settings.energy_mode = 'multi-group'
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 3
+    model.settings.particles = 10
+
+    source = openmc.IndependentSource()
+    source.space = openmc.stats.Point((5.0, 0.0, 0.0))
+    model.settings.source = source
+    return model
+
+def test_inverse_velocity(run_in_tmpdir, slab_model):
+    tally = openmc.Tally()
+    tally.scores = ['flux','inverse-velocity']
+    slab_model.tallies = [tally]
+    slab_model.run(apply_tally_results=True)
+    inverse_velocity = tally.mean.squeeze()[1]/tally.mean.squeeze()[0]
+
+    assert inverse_velocity == pytest.approx(1.6144e-5, rel=1e-4)