add kinetic simulation machinery

- add some missing precursor tally machinery in reaction.cpp and
  tally.cpp, and output.cpp
- add testing harness for kinetic simulations
- remove some unuesd timers
- implement time stepping loop in random_ray_simulation.cpp and
  supporting funtions.

Author: yardasol

include/openmc/random_ray/random_ray.h

@@ -50,6 +50,16 @@ class RandomRay : public Particle {
   static RandomRaySourceShape source_shape_; // Flag for linear source
   static RandomRaySampleMethod sample_method_; // Flag for sampling method
 
+  static double avg_miss_rate_;     // Average ray miss rate per
+                                    // iteration for reporting
+  static int64_t n_source_regions_; // Total number of source regions
+  static int64_t
+    n_external_source_regions_; // Total number of source regions with
+                                // non-zero external source terms
+  static uint64_t total_geometric_intersections_; // Tracks the total number of
+                                                  // geometric intersections by
+                                                  // all rays for reporting
+
   // Kinetic simulation variables
   static int bd_order_; // Order of backwards difference approximation for
                         // time-derivatives
@@ -60,9 +70,13 @@ class RandomRay : public Particle {
   //----------------------------------------------------------------------------
   // Public data members
   vector<float> angular_flux_;
-
   bool ray_trace_only_ {false}; // If true, only perform geometry operations
 
+  //---------------------------------------------------------------------------
+  // Public data members for kinetic simulations
+  vector<float> angular_flux_td_;
+  vector<float> angular_flux_td_prime_;
+
 private:
   //----------------------------------------------------------------------------
   // Private data members
@@ -77,6 +91,12 @@ class RandomRay : public Particle {
   double distance_travelled_ {0};
   bool is_active_ {false};
   bool is_alive_ {true};
+
+  //---------------------------------------------------------------------------
+  // Private data members for kinetic simulations
+  vector<double> delta_psi_td_;
+  vector<double> delta_psi_td_prime_;
+
 }; // class RandomRay
 
 } // namespace openmc

include/openmc/random_ray/random_ray_simulation.h

@@ -26,17 +26,15 @@ class RandomRaySimulation {
   void output_simulation_results() const;
   void instability_check(
     int64_t n_hits, double k_eff, double& avg_miss_rate) const;
-  void print_results_random_ray(uint64_t total_geometric_intersections,
-    double avg_miss_rate, int negroups, int64_t n_source_regions,
-    int64_t n_external_source_regions) const;
+  void print_results_random_ray() const;
 
   //----------------------------------------------------------------------------
   // Accessors
   FlatSourceDomain* domain() const { return domain_.get(); }
 
 private:
   //----------------------------------------------------------------------------
-  // Data members
+  // Data Members
 
   // Contains all flat source region data
   unique_ptr<FlatSourceDomain> domain_;
@@ -51,6 +49,9 @@ class RandomRaySimulation {
   // Number of energy groups
   int negroups_;
 
+  // Number of delay groups
+  int ndgroups_;
+
 }; // class RandomRaySimulation
 
 //============================================================================
@@ -61,6 +62,15 @@ void openmc_run_random_ray();
 void validate_random_ray_inputs();
 void openmc_reset_random_ray();
 
+//! Write data related to randaom ray to statepoint
+//! \param[in] group HDF5 group
+void write_random_ray_hdf5(hid_t group);
+
+// Functions for kinetic simulations
+void set_time_dependent_settings();
+void rename_statepoint_file(int i);
+void rename_tallies_file(int i);
+
 } // namespace openmc
 
-#endif // OPENMC_RANDOM_RAY_SIMULATION_H
+#endif // OPENMC_RANDOM_RAY_SIMUgATION_H

include/openmc/simulation.h

@@ -49,6 +49,10 @@ extern const RegularMesh* ufs_mesh;
 extern vector<double> k_generation;
 extern vector<int64_t> work_index;
 
+//-----------------------------------------------------------------------------
+// Global variables for kinetic simulations
+extern double current_time;
+
 } // namespace simulation
 
 //==============================================================================

include/openmc/timer.h

@@ -33,7 +33,6 @@ extern Timer time_event_collision;
 extern Timer time_event_death;
 extern Timer time_update_src;
 
-extern Timer time_initialize_td;
 extern Timer time_update_bd_vectors_td;
 extern Timer time_update_src_td;
 extern Timer time_compute_precursors;

openmc/examples.py

@@ -656,14 +656,14 @@ def slab_mg(num_regions=1, mat_names=None, mgxslib_name='2g.h5') -> openmc.Model
 
     return model
 
-def _generate_c5g7_materials(time_dependent) -> openmc.Materials:
+def _generate_c5g7_materials(kinetic) -> openmc.Materials:
     """Generate materials utilizing multi-group cross sections based on the
     the C5G7 Benchmark.
 
     Parameters
     ----------
-    time_dependent : bool
-        Flag to generate cross sections for a time-dependent model or not.
+    kinetic : bool
+        Flag to generate cross sections for kinetic simulations.
 
     Returns
     -------
@@ -692,7 +692,7 @@ def _generate_c5g7_materials(time_dependent) -> openmc.Materials:
     # come from Hou et al., "OECD/NEA benchmark for time-dependnet neutron
     # transport calculations without homogeniztaion"
     # DOI: 10.1016/j.nucengdes.2017.02.008
-    n_dg = C5G7_N_DG if time_dependent else 0
+    n_dg = C5G7_N_DG if kinetic else 0
 
     # Instantiate the 7-group (C5G7) cross section data
     uo2_xsdata = openmc.XSdata('UO2', groups, num_delayed_groups=n_dg)
@@ -728,7 +728,7 @@ def _generate_c5g7_materials(time_dependent) -> openmc.Materials:
                         0.0000e+00, 0.0000e+00])
 
     # Delayed and prompt cross sections for time-dependent simulation
-    if time_dependent:
+    if kinetic:
 
         # Table A2 in Hou et. al
         beta = np.array([[2.13333e-04, 2.13333e-04, 2.13333e-04, 2.13333e-04, 2.13333e-04, 2.13333e-04, 2.13333e-04],
@@ -804,7 +804,7 @@ def _generate_c5g7_materials(time_dependent) -> openmc.Materials:
     scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
     h2o_xsdata.set_scatter_matrix(scatter_matrix)
 
-    if time_dependent:
+    if kinetic:
         # Table A4 in Hou et al.
         velocities = np.array([2.23517E+09, 4.98880E+08, 3.84974E+07,
                               5.12639E+06, 1.67542E+06, 7.26031E+05, 2.81629E+05])
@@ -822,7 +822,7 @@ def _generate_c5g7_materials(time_dependent) -> openmc.Materials:
     uo2.set_density('macro', 1.0)
     uo2.add_macroscopic('UO2')
 
-    if time_dependent:
+    if kinetic:
         densities = np.linspace(1, 0.95, 100)
     else:
         densities = None
@@ -905,14 +905,14 @@ def _generate_random_ray_pin_cell(uo2, water) -> openmc.Universe:
 
     return pincell
 
-def random_ray_pin_cell(time_dependent=False) -> openmc.Model:
+def random_ray_pin_cell(kinetic=False) -> openmc.Model:
     """Create a PWR pin cell example using C5G7 cross section data.
     cross section data.
 
     Parameters
     ----------
-    time_dependent : bool
-        Flag to generate a time-dependent model or not.
+    kinetic : bool
+        Flag to generate kinetic simulation model or not.
 
     Returns
     -------
@@ -924,7 +924,7 @@ def random_ray_pin_cell(time_dependent=False) -> openmc.Model:
 
     ###########################################################################
     # Create Materials for the problem
-    materials = _generate_c5g7_materials(time_dependent)
+    materials = _generate_c5g7_materials(kinetic)
     uo2 = materials[0]
     water = materials[1]
 
@@ -962,9 +962,10 @@ def random_ray_pin_cell(time_dependent=False) -> openmc.Model:
     settings.random_ray['distance_inactive'] = 20.0
     settings.random_ray['ray_source'] = rr_source
     settings.random_ray['volume_normalized_flux_tallies'] = True
-    if time_dependent:
-        settings.run_mode = "time dependent"
-        settings.time_dependent = {
+    if kinetic:
+        settings.random_ray['bd_order'] = 3
+        settings.kinetic_simulation = True
+        settings.timestep_parameters = {
             "dt": 0.01,
             "n_timesteps": 20,
             "timestep_units": "s",
@@ -980,7 +981,7 @@ def random_ray_pin_cell(time_dependent=False) -> openmc.Model:
     # Instantiate a Tallies collection and export to XML
     tallies = openmc.Tallies([tally])
 
-    if time_dependent:
+    if kinetic:
         delay_filter = openmc.DelayedGroupFilter(np.arange(1, C5G7_N_DG+1, 1))
         tally = openmc.Tally(name="Delayed tally")
         tally.filters = [delay_filter]
@@ -997,7 +998,7 @@ def random_ray_pin_cell(time_dependent=False) -> openmc.Model:
     model.tallies = tallies
     return model
 
-def random_ray_lattice(time_dependent=False) -> openmc.Model:
+def random_ray_lattice(kinetic=False) -> openmc.Model:
     """Create a 2x2 PWR pin cell asymmetrical lattice example.
 
     This model is a 2x2 reflective lattice of fuel pins with one of the lattice
@@ -1006,8 +1007,8 @@ def random_ray_lattice(time_dependent=False) -> openmc.Model:
 
     Parameters
     ----------
-    time_dependent : bool
-        Flag to generate a time-dependent model or not.
+    kinetic : bool
+        Flag to generate a kinetic simulation model or not.
 
     Returns
     -------
@@ -1019,7 +1020,7 @@ def random_ray_lattice(time_dependent=False) -> openmc.Model:
 
     ###########################################################################
     # Create Materials for the problem
-    materials = _generate_c5g7_materials(time_dependent)
+    materials = _generate_c5g7_materials(kinetic)
     uo2 = materials[0]
     water = materials[1]
 
@@ -1031,7 +1032,7 @@ def random_ray_lattice(time_dependent=False) -> openmc.Model:
     # Define a moderator lattice universe
 
     moderator_infinite = openmc.Cell(name='moderator infinite')
-    if time_dependent:
+    if kinetic:
         water_reflector = water.clone()
         water_reflector.name='Water Reflector'
         water_reflector.set_density('macro', 1.0)
@@ -1094,9 +1095,10 @@ def random_ray_lattice(time_dependent=False) -> openmc.Model:
     settings.random_ray['distance_inactive'] = 20.0
     settings.random_ray['ray_source'] = rr_source
     settings.random_ray['volume_normalized_flux_tallies'] = True
-    if time_dependent:
-        settings.run_mode = "time dependent"
-        settings.time_dependent = {
+    if kinetic:
+        settings.random_ray['bd_order'] = 3
+        settings.kinetic_simulation = True
+        settings.timestep_parameters = {
             "dt": 0.01,
             "n_timesteps": 2,
             "timestep_units": "s",
@@ -1127,7 +1129,7 @@ def random_ray_lattice(time_dependent=False) -> openmc.Model:
     # Instantiate a Tallies collection and export to XML
     tallies = openmc.Tallies([tally])
 
-    if time_dependent:
+    if kinetic:
         delay_filter = openmc.DelayedGroupFilter(np.arange(1, C5G7_N_DG+1, 1))
         tally = openmc.Tally(name="Mesh delayed tally")
         tally.filters = [mesh_filter, delay_filter]

src/output.cpp

@@ -598,6 +598,7 @@ const std::unordered_map<int, const char*> score_names = {
   {SCORE_IFP_TIME_NUM, "IFP lifetime numerator"},
   {SCORE_IFP_BETA_NUM, "IFP delayed fraction numerator"},
   {SCORE_IFP_DENOM, "IFP common denominator"},
+  {SCORE_PRECURSORS, "Precursor Population"},
 };
 
 //! Create an ASCII output file showing all tally results.

src/random_ray/flat_source_domain.cpp

@@ -96,10 +96,12 @@ void FlatSourceDomain::batch_reset()
     source_regions_.scalar_flux_new(se) = 0.0;
   }
 
-  if (settings::kinetic_simulation && !settings::is_initial_condition) {
+  if (settings::kinetic_simulation) {
 #pragma omp parallel for
     for (int64_t se = 0; se < n_source_elements(); se++) {
-      source_regions_.scalar_flux_td_new(se) = 0.0;
+      source_regions_.precursors_new(se) = 0.0;
+      if (!settings::is_initial_condition)
+        source_regions_.scalar_flux_td_new(se) = 0.0;
     }
   }
 }