Add MG particle speed getter and delayed fission neutron emission

docs/source/methods/neutron_physics.rst

@@ -290,7 +290,10 @@ create and store fission sites for the following generation. First, the average
 number of prompt and delayed neutrons must be determined to decide whether the
 secondary neutrons will be prompt or delayed. This is important because delayed
 neutrons have a markedly different spectrum from prompt neutrons, one that has a
-lower average energy of emission. The total number of neutrons emitted
+lower average energy of emission. Furthermore, in simulations where tracking
+time of neutrons is important, we need to consider the emission time delay of
+the secondary neutrons, which is dependent on the decay constant of the
+delayed neutron precursor. The total number of neutrons emitted
 :math:`\nu_t` is given as a function of incident energy in the ENDF format. Two
 representations exist for :math:`\nu_t`. The first is a polynomial of order
 :math:`N` with coefficients :math:`c_0,c_1,\dots,c_N`. If :math:`\nu_t` has this
@@ -306,8 +309,8 @@ interpolation law. The number of prompt neutrons released per fission event
 :math:`\nu_p` is also given as a function of incident energy and can be
 specified in a polynomial or tabular format. The number of delayed neutrons
 released per fission event :math:`\nu_d` can only be specified in a tabular
-format. In practice, we only need to determine :math:`nu_t` and
-:math:`nu_d`. Once these have been determined, we can calculated the delayed
+format. In practice, we only need to determine :math:`\nu_t` and
+:math:`\nu_d`. Once these have been determined, we can calculated the delayed
 neutron fraction
 
 .. math::
@@ -335,8 +338,14 @@ neutrons. Otherwise, we produce :math:`\lfloor \nu \rfloor + 1` neutrons. Then,
 for each fission site produced, we sample the outgoing angle and energy
 according to the algorithms given in :ref:`sample-angle` and
 :ref:`sample-energy` respectively. If the neutron is to be born delayed, then
-there is an extra step of sampling a delayed neutron precursor group since they
-each have an associated secondary energy distribution.
+there is an extra step of sampling a delayed neutron precursor group to get the
+associated secondary energy distribution and the decay constant
+:math:`\lambda`, which is needed to sample the emission delay time :math:`t_d`:
+
+.. math::
+    :label: sample-delay-time
+
+    t_d = -\frac{\ln \xi}{\lambda}.
 
 The sampled outgoing angle and energy of fission neutrons along with the
 position of the collision site are stored in an array called the fission

include/openmc/settings.h

@@ -123,6 +123,7 @@ extern std::unordered_set<int>
 extern int max_splits; //!< maximum number of particle splits for weight windows
 extern int64_t max_surface_particles; //!< maximum number of particles to be
                                       //!< banked on surfaces per process
+extern int max_secondaries; //!< maximum number of secondaries in the bank
 extern TemperatureMethod
   temperature_method; //!< method for choosing temperatures
 extern double

openmc/settings.py

@@ -114,6 +114,10 @@ class Settings:
     max_splits : int
         Maximum number of times a particle can split during a history
 
+        .. versionadded:: 0.15
+    max_secondaries : int
+        Maximum secondary bank size
+
         .. versionadded:: 0.13
     max_tracks : int
         Maximum number of tracks written to a track file (per MPI process).
@@ -341,6 +345,7 @@ def __init__(self, **kwargs):
         self._weight_windows_file = None
         self._weight_window_checkpoints = {}
         self._max_splits = None
+        self._max_secondaries = None
         self._max_tracks = None
 
         for key, value in kwargs.items():
@@ -986,6 +991,16 @@ def max_splits(self, value: int):
         cv.check_greater_than('max particle splits', value, 0)
         self._max_splits = value
 
+    @property
+    def max_secondaries(self) -> int:
+        return self._max_secondaries
+
+    @max_secondaries.setter
+    def max_secondaries(self, value: int):
+        cv.check_type('maximum secondary bank size', value, Integral)
+        cv.check_greater_than('max secondary bank size', value, 0)
+        self._max_secondaries = value
+
     @property
     def max_tracks(self) -> int:
         return self._max_tracks
@@ -1417,6 +1432,11 @@ def _create_max_splits_subelement(self, root):
             elem = ET.SubElement(root, "max_splits")
             elem.text = str(self._max_splits)
 
+    def _create_max_secondaries_subelement(self, root):
+        if self._max_secondaries is not None:
+            elem = ET.SubElement(root, "max_secondaries")
+            elem.text = str(self._max_secondaries)
+
     def _create_max_tracks_subelement(self, root):
         if self._max_tracks is not None:
             elem = ET.SubElement(root, "max_tracks")
@@ -1763,6 +1783,11 @@ def _max_splits_from_xml_element(self, root):
         if text is not None:
             self.max_splits = int(text)
 
+    def _max_secondaries_from_xml_element(self, root):
+        text = get_text(root, 'max_secondaries')
+        if text is not None:
+            self.max_secondaries = int(text)
+
     def _max_tracks_from_xml_element(self, root):
         text = get_text(root, 'max_tracks')
         if text is not None:
@@ -1828,6 +1853,7 @@ def to_xml_element(self, mesh_memo=None):
         self._create_weight_windows_file_element(element)
         self._create_weight_window_checkpoints_subelement(element)
         self._create_max_splits_subelement(element)
+        self._create_max_secondaries_subelement(element)
         self._create_max_tracks_subelement(element)
 
         # Clean the indentation in the file to be user-readable
@@ -1930,6 +1956,7 @@ def from_xml_element(cls, elem, meshes=None):
         settings._weight_window_generators_from_xml_element(elem, meshes)
         settings._weight_window_checkpoints_from_xml_element(elem)
         settings._max_splits_from_xml_element(elem)
+        settings._max_secondaries_from_xml_element(elem)
         settings._max_tracks_from_xml_element(elem)
 
         # TODO: Get volume calculations

src/finalize.cpp

@@ -91,6 +91,7 @@ int openmc_finalize()
   settings::max_lost_particles = 10;
   settings::max_order = 0;
   settings::max_particles_in_flight = 100000;
+  settings::max_secondaries = 10000;
   settings::max_splits = 1000;
   settings::max_tracks = 1000;
   settings::max_write_lost_particles = -1;

src/particle.cpp

@@ -44,6 +44,15 @@ namespace openmc {
 
 double Particle::speed() const
 {
+  // Multigroup speed?
+  if (!settings::run_CE) {
+    auto& macro_xs = data::mg.macro_xs_[material()];
+    int macro_t = mg_xs_cache().t;
+    int macro_a = macro_xs.get_angle_index(u());
+    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, g(), nullptr,
+                   nullptr, nullptr, macro_t, macro_a);
+  }
+
   // Determine mass in eV/c^2
   double mass;
   switch (this->type()) {
@@ -237,6 +246,9 @@ void Particle::event_advance()
     double push_back_distance = speed() * dt;
     this->move_distance(-push_back_distance);
     hit_time_boundary = true;
+
+    // Reduce the distance traveled for tallying
+    distance -= push_back_distance;
   }
 
   // Score track-length tallies

src/physics.cpp

@@ -113,7 +113,7 @@ void sample_neutron_reaction(Particle& p)
 
       // Make sure particle population doesn't grow out of control for
       // subcritical multiplication problems.
-      if (p.secondary_bank().size() >= 10000) {
+      if (p.secondary_bank().size() >= settings::max_secondaries) {
         fatal_error(
           "The secondary particle bank appears to be growing without "
           "bound. You are likely running a subcritical multiplication problem "
@@ -209,6 +209,20 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
     // Sample delayed group and angle/energy for fission reaction
     sample_fission_neutron(i_nuclide, rx, &site, p);
 
+    // If delayed product production, sample time of emission
+    if (site.delayed_group > 0) {
+      double decay_rate = rx.products_[site.delayed_group].decay_rate_;
+      site.time -= std::log(prn(p.current_seed())) / decay_rate;
+
+      // Reject site if it exceeds time cutoff
+      double time_cutoff =
+        settings::time_cutoff[static_cast<int>(site.particle)];
+      if (site.time > time_cutoff) {
+        p.n_progeny()--;
+        continue;
+      }
+    }
+
     // Store fission site in bank
     if (use_fission_bank) {
       int64_t idx = simulation::fission_bank.thread_safe_append(site);

src/physics_mg.cpp

@@ -127,6 +127,7 @@ void create_fission_sites(Particle& p)
     SourceSite site;
     site.r = p.r();
     site.particle = ParticleType::neutron;
+    site.time = p.time();
     site.wgt = 1. / weight;
     site.parent_id = p.id();
     site.progeny_id = p.n_progeny()++;
@@ -154,6 +155,22 @@ void create_fission_sites(Particle& p)
     // of the code, 0 is prompt.
     site.delayed_group = dg + 1;
 
+    // If delayed product production, sample time of emission
+    if (dg != -1) {
+      auto& macro_xs = data::mg.macro_xs_[p.material()];
+      double decay_rate =
+        macro_xs.get_xs(MgxsType::DECAY_RATE, 0, nullptr, nullptr, &dg, 0, 0);
+      site.time -= std::log(prn(p.current_seed())) / decay_rate;
+
+      // Reject site if it exceeds time cutoff
+      double time_cutoff =
+        settings::time_cutoff[static_cast<int>(site.particle)];
+      if (site.time > time_cutoff) {
+        p.n_progeny()--;
+        continue;
+      }
+    }
+
     // Store fission site in bank
     if (use_fission_bank) {
       int64_t idx = simulation::fission_bank.thread_safe_append(site);

src/settings.cpp

@@ -105,6 +105,7 @@ int max_order {0};
 int n_log_bins {8000};
 int n_batches;
 int n_max_batches;
+int max_secondaries {10000};
 int max_splits {1000};
 int max_tracks {1000};
 ResScatMethod res_scat_method {ResScatMethod::rvs};
@@ -930,6 +931,10 @@ void read_settings_xml(pugi::xml_node root)
     weight_windows_on = get_node_value_bool(root, "weight_windows_on");
   }
 
+  if (check_for_node(root, "max_secondaries")) {
+    settings::max_secondaries = std::stoi(get_node_value(root, "max_secondaries"));
+  }
+
   if (check_for_node(root, "max_splits")) {
     settings::max_splits = std::stoi(get_node_value(root, "max_splits"));
   }