fast forward to use updated functions; add some missing flow control

Author: yardasol

include/openmc/random_ray/random_ray_simulation.h

@@ -21,9 +21,10 @@ class RandomRaySimulation {
   // Methods
   void apply_fixed_sources_and_mesh_domains();
   void prepare_fixed_sources_adjoint();
-  void kinetic_single_time_step(int i);
   void prepare_adjoint_simulation();
   void simulate();
+  void initialize_time_step(int i);
+  void finalize_time_step();
   void output_simulation_results() const;
   void instability_check(
     int64_t n_hits, double k_eff, double& avg_miss_rate) const;

include/openmc/random_ray/source_region.h

@@ -996,6 +996,7 @@ class SourceRegionContainer {
   bool& is_linear() { return is_linear_; }
   const bool is_linear() const { return is_linear_; }
   SourceRegionHandle get_source_region_handle(int64_t sr);
+  void simulation_reset();
   void adjoint_reset();
 
   //---------------------------------------------------------------------------

src/random_ray/random_ray_simulation.cpp

@@ -438,6 +438,9 @@ void RandomRaySimulation::prepare_fixed_sources_adjoint()
 
 void RandomRaySimulation::prepare_adjoint_simulation()
 {
+  // Reset all simulation values
+  domain_->source_regions_.simulation_reset();
+
   // Configure the domain for adjoint simulation
   FlatSourceDomain::adjoint_ = true;
 
@@ -454,59 +457,6 @@ void RandomRaySimulation::prepare_adjoint_simulation()
   domain_->nu_sigma_f_.swap(domain_->chi_);
 }
 
-// TODO: Add support for time-dependent restart
-void RandomRaySimulation::kinetic_single_time_step(int i)
-{
-  // Increment time step
-  simulation::current_timestep = i + 1;
-  if (i >= 0)
-    // Increment the current time
-    simulation::current_time += settings::dt;
-
-  // 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();
-
-  if (i >= 0) {
-    // Compute RHS backward differences
-    domain_->compute_rhs_bd_quantities();
-
-    // Update time dependent cross section based on the density
-    domain_->update_material_density(i);
-  }
-
-  // Run the initial condition
-  simulate();
-
-  if (i == -1) {
-    // Initialize the BD arrays if initial condition
-    domain_->store_time_step_quantities(false);
-    // 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)
-    rename_time_step_file("tallies", ".out", simulation::current_timestep);
-}
-
 void RandomRaySimulation::simulate()
 {
   if (!is_first_simulation_) {
@@ -580,11 +530,11 @@ void RandomRaySimulation::simulate()
       domain_->apply_transport_stabilization();
 
       if (settings::run_mode == RunMode::EIGENVALUE) {
-        // Compute random ray k-eff
-        if (!settings::kinetic_simulation ||
-            settings::kinetic_simulation && simulation::is_initial_condition) {
+        // Compute random ray k-eff for initial condition.
+        // This keff will be preserved
+        if (simulation::is_initial_condition) {
           domain_->compute_k_eff();
-          if (simulation::k_eff_correction) {
+          if (settings::kinetic_simulation && simulation::k_eff_correction) {
             static_fission_rate_.push_back(domain_->fission_rate_);
             static_k_eff_.push_back(domain_->k_eff_);
           }
@@ -655,6 +605,51 @@ void RandomRaySimulation::simulate()
     is_first_simulation_ = false;
 }
 
+void RandomRaySimulation::initialize_time_step(int i)
+{
+  if (simulation::k_eff_correction)
+    static_avg_k_eff_ = simulation::keff;
+  domain_->k_eff_ = static_avg_k_eff_;
+
+  // Increment current timestep and simuation time
+  simulation::current_timestep = i + 1;
+
+  // Propagate previous converted solution for kinetic simulation
+  domain_->source_regions_.simulation_reset();
+  domain_->propagate_final_quantities();
+  domain_->source_regions_.time_step_reset();
+
+  if (!simulation::is_initial_condition) {
+    // Compute RHS backward differences
+    domain_->compute_rhs_bd_quantities();
+
+    // Update time dependent cross section based on the density
+    domain_->update_material_density(i);
+
+    simulation::current_time += settings::dt;
+  }
+}
+
+void RandomRaySimulation::finalize_time_step()
+{
+  if (simulation::is_initial_condition) {
+    // Initialize the BD arrays if initial condition
+    domain_->store_time_step_quantities(false);
+    // Toggle off initial condition and source correction
+    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)
+    rename_time_step_file("tallies", ".out", simulation::current_timestep);
+}
+
 void RandomRaySimulation::output_simulation_results() const
 {
   // Print random ray results
@@ -861,7 +856,6 @@ void openmc_run_random_ray()
   //////////////////////////////////////////////////////////
   // Run forward simulation
   //////////////////////////////////////////////////////////
-
   if (openmc::mpi::master) {
     if (openmc::FlatSourceDomain::adjoint_) {
       openmc::FlatSourceDomain::adjoint_ = false;
@@ -883,14 +877,22 @@ void openmc_run_random_ray()
   // Initialize fixed sources, if present
   sim.apply_fixed_sources_and_mesh_domains();
 
-  // Run initial random ray simulation
+  // Simulate single random ray simulation (static case)
+  // OR get an initial estimate for scattering and fission
+  // distributions (if fissile material exist),
+  // and k-eff (if kinetic eigenvalue simulation)
   sim.simulate();
 
   if (openmc::settings::kinetic_simulation) {
-    // Timestepping loop, including k-eff correction initial
-    // condition (i = -1)
-    for (int i = -1; i < openmc::settings::n_timesteps; i++)
-      sim.kinetic_single_time_step(i);
+    // Toggle initial condition source correction
+    openmc::simulation::k_eff_correction = true;
+    int i_start = -1;
+    // Timestepping loop,
+    for (int i = i_start; i < openmc::settings::n_timesteps; i++) {
+      sim.initialize_time_step(i);
+      sim.simulate();
+      sim.finalize_time_step();
+    }
   }
 
   //////////////////////////////////////////////////////////

src/random_ray/source_region.cpp

@@ -364,17 +364,14 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr)
   return handle;
 }
 
-void SourceRegionContainer::adjoint_reset()
+void SourceRegionContainer::simulation_reset()
 {
   std::fill(n_hits_.begin(), n_hits_.end(), 0);
   std::fill(volume_.begin(), volume_.end(), 0.0);
   std::fill(volume_t_.begin(), volume_t_.end(), 0.0);
   std::fill(volume_sq_.begin(), volume_sq_.end(), 0.0);
   std::fill(volume_sq_t_.begin(), volume_sq_t_.end(), 0.0);
   std::fill(volume_naive_.begin(), volume_naive_.end(), 0.0);
-  std::fill(
-    external_source_present_.begin(), external_source_present_.end(), 0);
-  std::fill(external_source_.begin(), external_source_.end(), 0.0);
   std::fill(centroid_.begin(), centroid_.end(), Position {0.0, 0.0, 0.0});
   std::fill(centroid_iteration_.begin(), centroid_iteration_.end(),
     Position {0.0, 0.0, 0.0});
@@ -390,7 +387,6 @@ void SourceRegionContainer::adjoint_reset()
   }
   std::fill(scalar_flux_new_.begin(), scalar_flux_new_.end(), 0.0);
   std::fill(source_.begin(), source_.end(), 0.0f);
-  std::fill(external_source_.begin(), external_source_.end(), 0.0f);
   std::fill(source_gradients_.begin(), source_gradients_.end(),
     MomentArray {0.0, 0.0, 0.0});
   std::fill(flux_moments_old_.begin(), flux_moments_old_.end(),
@@ -424,6 +420,13 @@ void SourceRegionContainer::adjoint_reset()
   }
 }
 
+void SourceRegionContainer::adjoint_reset()
+{
+  std::fill(
+    external_source_present_.begin(), external_source_present_.end(), 0);
+  std::fill(external_source_.begin(), external_source_.end(), 0.0f);
+}
+
 //-----------------------------------------------------------------------------
 // Methods for kinetic simulations
 

src/settings.cpp

@@ -381,6 +381,9 @@ void get_run_parameters(pugi::xml_node node_base)
     if (check_for_node(random_ray_node, "adjoint")) {
       FlatSourceDomain::adjoint_ =
         get_node_value_bool(random_ray_node, "adjoint");
+      if (FlatSourceDomain::adjoint_ && kinetic_simulation) {
+        fatal_error("Adjoint kinetic simulations are currently unsupported .");
+      }
     }
     if (check_for_node(random_ray_node, "sample_method")) {
       std::string temp_str =
@@ -1343,6 +1346,10 @@ void read_settings_xml(pugi::xml_node root)
 
   // Create weight window generator objects
   if (check_for_node(root, "weight_window_generators")) {
+    if (kinetic_simulation) {
+      fatal_error("Weight window generation is currently unsupported in kinetic"
+                  " random ray solver mode.");
+    }
     auto wwgs_node = root.child("weight_window_generators");
     for (pugi::xml_node node_wwg :
       wwgs_node.children("weight_windows_generator")) {