Approximate multigroup velocity (#3766)

Co-authored-by: Adam Nelson <1037107+nelsonag@users.noreply.github.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

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,48 @@ 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 such data is not available,
+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{\alpha}{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:`\alpha` 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::
+
+   \frac{1}{v_g} = \frac{1}{c \log\left(\frac{E_{\text{max}}^g}{E_{\text{min}}^g}\right)}
+   \left[ 2(\operatorname{arctanh}(k_{\text{max}}^{-1}) - \operatorname{arctanh}(k_{\text{min}}^{-1}))
+   - (k_{\text{max}}-k_{\text{min}}) \right]
+
+where :math:`c` is the speed of light and :math:`k_{\text{max}}`,
+:math:`k_{\text{min}}` are defined by a change of variables:
+
+.. math::
+
+   k = \sqrt{1+\frac{2 m_n c^2}{E}}
+
+where :math:`E` is the particle kinetic energy and :math:`m_n` is the neutron
+rest mass.
+
 .. _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 @@ public:
   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 @@ public:
 
   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 @@ class EnergyGroups:
     @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,19 @@ 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 alpha = 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 =
+      alpha * (2.0 * (std::atanh(1.0 / k_max) - std::atanh(1.0 / k_min)) -
+                (k_max - k_min));
+    default_inverse_velocity_.push_back(inv_v);
+  }
+
   // Close MGXS HDF5 file
   file_close(file_id);
 }

src/particle.cpp

@@ -48,26 +48,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)