cleanup

Author: yardasol

include/openmc/random_ray/flat_source_domain.h

@@ -191,7 +191,7 @@ class FlatSourceDomain {
   // 3D arrays stored in 1D representing values for all materials x energy
   // groups x delay groups
   vector<double> nu_d_sigma_f_;
-  vector<double> chi_d_;
+  vector<double> chi_d_lambda_; // chi-delayed * lambda in each delay group
 
   // 2D arrays stored in 1D representing values for all materials x energy
   // groups

include/openmc/random_ray/source_region.h

@@ -209,10 +209,10 @@ class SourceRegionHandle {
 
   // Energy group-wise 2D BD arrays (g x time step)
   std::deque<double>* scalar_flux_bd_;
-  std::deque<double>* precursors_bd_;
+  std::deque<float>* source_bd_;
 
   // Delay group-wise 2D BD arrays (dg x time step)
-  std::deque<float>* source_bd_;
+  std::deque<double>* precursors_bd_;
 
   // Energy group-wise 1D RHS BD arrays
   double* scalar_flux_rhs_bd_;
@@ -374,15 +374,15 @@ class SourceRegionHandle {
     return scalar_flux_bd_[g];
   }
 
+  std::deque<float>& source_bd(int g) { return source_bd_[g]; }
+  const std::deque<float>& source_bd(int g) const { return source_bd_[g]; }
+
   std::deque<double>& precursors_bd(int dg) { return precursors_bd_[dg]; }
   const std::deque<double>& precursors_bd(int dg) const
   {
     return precursors_bd_[dg];
   }
 
-  std::deque<float>& source_bd(int g) { return source_bd_[g]; }
-  const std::deque<float>& source_bd(int g) const { return source_bd_[g]; }
-
   double& scalar_flux_rhs_bd(int g) { return scalar_flux_rhs_bd_[g]; }
   const double scalar_flux_rhs_bd(int g) const
   {
@@ -522,17 +522,17 @@ class SourceRegion {
                      //!< finite number of previous time steps (used for
                      //!< computing the first order (and second order, for SDP)
                      //!< scalar flux time derivative)
-  vector<std::deque<double>>
-    precursors_bd_; //!< The final scalar flux in each energy group from a
-                    //!< finite number of previous time steps (used for
-                    //!< computing the first order source time derivative for
-                    //!< SDP)
+  vector<std::deque<float>>
+    source_bd_; //!< The final scalar flux in each energy group from a
+                //!< finite number of previous time steps (used for
+                //!< computing the first order source time derivative for
+                //!< SDP)
 
   // Delay group-wise 2D BD arrays (dg x time step)
-  vector<std::deque<float>>
-    source_bd_; //!< The final precursor population in each energy group from a
-                //!< finite number of previous time steps (used for computing
-                //!< the first order precursor time derivative)
+  vector<std::deque<double>>
+    precursors_bd_; //!< The final precursor population in each energy group
+                    //!< from a finite number of previous time steps (used for
+                    //!< computing the first order precursor time derivative)
 
   // Energy group-wise 1D RHS BD arrays
   vector<double>
@@ -877,20 +877,6 @@ class SourceRegionContainer {
     return scalar_flux_bd_[se];
   }
 
-  std::deque<double>& precursors_bd(int64_t sr, int dg)
-  {
-    return precursors_bd_[dindex(sr, dg)];
-  }
-  const std::deque<double>& precursors_bd(int64_t sr, int dg) const
-  {
-    return precursors_bd_[dindex(sr, dg)];
-  }
-  std::deque<double>& precursors_bd(int64_t de) { return precursors_bd_[de]; }
-  const std::deque<double>& precursors_bd(int64_t de) const
-  {
-    return precursors_bd_[de];
-  }
-
   std::deque<float>& source_bd(int64_t sr, int g)
   {
     return source_bd_[index(sr, g)];
@@ -905,6 +891,20 @@ class SourceRegionContainer {
     return source_bd_[se];
   }
 
+  std::deque<double>& precursors_bd(int64_t sr, int dg)
+  {
+    return precursors_bd_[dindex(sr, dg)];
+  }
+  const std::deque<double>& precursors_bd(int64_t sr, int dg) const
+  {
+    return precursors_bd_[dindex(sr, dg)];
+  }
+  std::deque<double>& precursors_bd(int64_t de) { return precursors_bd_[de]; }
+  const std::deque<double>& precursors_bd(int64_t de) const
+  {
+    return precursors_bd_[de];
+  }
+
   double& scalar_flux_rhs_bd(int64_t sr, int g)
   {
     return scalar_flux_rhs_bd_[index(sr, g)];
@@ -1067,11 +1067,13 @@ class SourceRegionContainer {
   vector<double> precursors_new_;
   vector<double> precursors_final_;
 
-  // SoA energy group-wise 3D BD arrays (sr x g/dg X timestep) flattened to 2D
+  // SoA energy group-wise 3D BD arrays (sr x g X timestep) flattened to 2D
   vector<std::deque<double>> scalar_flux_bd_;
-  vector<std::deque<double>> precursors_bd_;
   vector<std::deque<float>> source_bd_;
 
+  // SoA delay group-wise 3D BD arrays (sr x dg X timestep) flattened to 2D
+  vector<std::deque<double>> precursors_bd_;
+
   // SoA energy group-wise 2D RHS BD arrays flattened to 1D
   vector<double> scalar_flux_rhs_bd_;
   vector<float> source_rhs_bd_;

src/random_ray/flat_source_domain.cpp

@@ -158,11 +158,10 @@ void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
           settings::create_delayed_neutrons) {
         double delayed_source = 0.0;
         for (int dg = 0; dg < ndgroups_; dg++) {
-          double chi_d =
-            chi_d_[material * negroups_ * ndgroups_ + dg * negroups_ + g_out];
-          double lambda = lambda_[material * ndgroups_ + dg];
+          double chi_d_lambda = chi_d_lambda_[material * negroups_ * ndgroups_ +
+                                              dg * negroups_ + g_out];
           double precursors = srh.precursors_old(dg);
-          delayed_source += chi_d * precursors * lambda;
+          delayed_source += chi_d_lambda * precursors;
         }
         total_source += delayed_source;
       }
@@ -1400,13 +1399,13 @@ void FlatSourceDomain::flatten_xs()
               // material is fissionable but has very small sigma_f.
               chi_d = 0.0;
             }
-            chi_d_.push_back(chi_d);
+            chi_d_lambda_.push_back(chi_d * lambda);
           }
         } else {
           lambda_.push_back(0);
           for (int g_out = 0; g_out < negroups_; g_out++) {
             nu_d_sigma_f_.push_back(0);
-            chi_d_.push_back(0);
+            chi_d_lambda_.push_back(0);
           }
         }
       }

src/random_ray/random_ray_simulation.cpp

@@ -510,24 +510,23 @@ void RandomRaySimulation::prepare_adjoint_simulation()
 // TODO: Add support for time-dependent restart
 void RandomRaySimulation::kinetic_single_time_step(int i)
 {
-  if (i == -1) {
-    // Set flag for k_eff correction if initial condition
-    simulation::k_eff_correction = true;
-
-    // Store average keff from initial simulation
-    static_avg_k_eff_ = simulation::keff;
-  }
-
   // Increment time step
   simulation::current_timestep = i + 1;
-  if (i == -1)
-    // Current time is zero for initial condition
-    simulation::current_time = 0;
-  else
-    // Else, increment the current time
+  if (i >= 0)
+    // Increment the current time
     simulation::current_time += settings::dt;
 
-  domain_->k_eff_ = static_avg_k_eff_;
+  // Set eigenvalue if needed
+  if (settings::run_mode == RunMode::EIGENVALUE) {
+    if (i == -1) {
+      // Set flag for k_eff correction if initial condition
+      simulation::k_eff_correction = true;
+
+      // Store average keff from initial simulation
+      static_avg_k_eff_ = simulation::keff;
+    }
+    domain_->k_eff_ = static_avg_k_eff_;
+  }
   domain_->source_regions_.adjoint_reset();
   domain_->propagate_final_quantities();
   domain_->source_regions_.time_step_reset();
@@ -543,25 +542,22 @@ void RandomRaySimulation::kinetic_single_time_step(int i)
   // Run the initial condition
   simulate();
 
-  if (i == -1)
+  if (i == -1) {
     // Initialize the BD arrays if initial condition
     domain_->store_time_step_quantities(false);
-  else
+    // Reset flags for kinetic simulation if initial condition
+    simulation::is_initial_condition = false;
+    simulation::k_eff_correction = false;
+  } else {
     // Else, store final quantities for the current time step
     domain_->store_time_step_quantities();
+  }
 
   // Rename statepoint and tallies file for the current time step
   rename_time_step_file(fmt::format("statepoint.{0}", settings::n_batches),
     ".h5", simulation::current_timestep);
-  if (settings::output_tallies) {
+  if (settings::output_tallies)
     rename_time_step_file("tallies", ".out", simulation::current_timestep);
-  }
-
-  if (i == -1) {
-    // Reset flags for kinetic simulation if initial condition
-    simulation::is_initial_condition = false;
-    simulation::k_eff_correction = false;
-  }
 }
 
 void RandomRaySimulation::simulate()

src/simulation.cpp

@@ -338,7 +338,7 @@ vector<int64_t> work_index;
 
 bool is_initial_condition {true};
 int current_timestep;
-double current_time;
+double current_time {0.0};
 bool k_eff_correction {false};
 
 } // namespace simulation