move some global vars from settings.cpp to simulation.cpp

- is_initial_condition
- current_timetstep

Author: yardasol

include/openmc/settings.h

@@ -203,11 +203,6 @@ extern double weight_survive;      //!< Survival weight after Russian roulette
 // Timestep variables for kinetic simulation
 extern int n_timesteps;      //!< number of timesteps
 extern double dt;            //!< fixed timestep size
-extern int current_timestep; //!< current timestep index
-extern bool
-  is_initial_condition; //!< if eigenvalue/fixed source sim is an initial
-                        // condition for a time-dependent simulation
-
 } // namespace settings
 
 //==============================================================================

include/openmc/simulation.h

@@ -51,7 +51,11 @@ extern vector<int64_t> work_index;
 
 //-----------------------------------------------------------------------------
 // Global variables for kinetic simulations
-extern double current_time;
+extern bool
+  is_initial_condition;      //!< if eigenvalue/fixed source sim is an initial
+                             // condition for a kinetic simulation
+extern int current_timestep; // !< current time step in kinetic simulation
+extern double current_time;  // !< current time in kinetic simulation
 
 } // namespace simulation
 

src/random_ray/flat_source_domain.cpp

@@ -100,7 +100,7 @@ void FlatSourceDomain::batch_reset()
 #pragma omp parallel for
     for (int64_t se = 0; se < n_source_elements(); se++) {
       source_regions_.precursors_new(se) = 0.0;
-      if (!settings::is_initial_condition)
+      if (!simulation::is_initial_condition)
         source_regions_.scalar_flux_td_new(se) = 0.0;
     }
   }
@@ -156,7 +156,7 @@ void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
     }
   }
   // Set souce_td to source for IC calculation
-  if (settings::is_initial_condition) {
+  if (simulation::is_initial_condition) {
 #pragma omp parallel for
     for (int g = 0; g < negroups_; g++)
       srh.source_td(g) += srh.source(g);
@@ -191,7 +191,7 @@ void FlatSourceDomain::normalize_scalar_flux_and_volumes(
 #pragma omp parallel for
   for (int64_t se = 0; se < n_source_elements(); se++) {
     source_regions_.scalar_flux_new(se) *= normalization_factor;
-    if (settings::kinetic_simulation && !settings::is_initial_condition)
+    if (settings::kinetic_simulation && !simulation::is_initial_condition)
       source_regions_.scalar_flux_td_new(se) *= normalization_factor;
   }
 
@@ -226,7 +226,7 @@ void FlatSourceDomain::set_flux_to_flux_plus_source(
     double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g];
     source_regions_.scalar_flux_new(sr, g) /= (sigma_t * volume);
     source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g);
-    if (settings::kinetic_simulation && !settings::is_initial_condition) {
+    if (settings::kinetic_simulation && !simulation::is_initial_condition) {
       double sigma_t_td =
         sigma_t_td_[source_regions_.material(sr) * negroups_ + g];
       source_regions_.scalar_flux_td_new(sr, g) /= (sigma_t_td * volume);
@@ -252,7 +252,7 @@ void FlatSourceDomain::set_flux_to_old_flux(int64_t sr, int g)
 {
   source_regions_.scalar_flux_new(sr, g) =
     source_regions_.scalar_flux_old(sr, g);
-  if (settings::kinetic_simulation && !settings::is_initial_condition) {
+  if (settings::kinetic_simulation && !simulation::is_initial_condition) {
     source_regions_.scalar_flux_td_new(sr, g) =
       source_regions_.scalar_flux_td_old(sr, g);
   }
@@ -261,7 +261,7 @@ void FlatSourceDomain::set_flux_to_old_flux(int64_t sr, int g)
 void FlatSourceDomain::set_flux_to_source(int64_t sr, int g)
 {
   source_regions_.scalar_flux_new(sr, g) = source_regions_.source(sr, g);
-  if (settings::kinetic_simulation && !settings::is_initial_condition) {
+  if (settings::kinetic_simulation && !simulation::is_initial_condition) {
     source_regions_.scalar_flux_td_new(sr, g) =
       source_regions_.source_td(sr, g);
   }
@@ -344,7 +344,7 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux()
       }
       // Set time-dependent flux to unperturbed flux during the initial
       // condition calculation
-      if (settings::is_initial_condition) {
+      if (simulation::is_initial_condition) {
         source_regions_.scalar_flux_td_new(sr, g) =
           source_regions_.scalar_flux_new(sr, g);
       }
@@ -751,7 +751,7 @@ void FlatSourceDomain::random_ray_tally()
     double material = source_regions_.material(sr);
     for (int g = 0; g < negroups_; g++) {
       double flux;
-      if (settings::kinetic_simulation && !settings::is_initial_condition) {
+      if (settings::kinetic_simulation && !simulation::is_initial_condition) {
         flux = source_regions_.scalar_flux_td_new(sr, g) *
                source_normalization_factor;
       } else {
@@ -771,7 +771,8 @@ void FlatSourceDomain::random_ray_tally()
         case SCORE_TOTAL:
           if (material != MATERIAL_VOID) {
             double sigma_t;
-            if (settings::kinetic_simulation && !settings::is_initial_condition)
+            if (settings::kinetic_simulation &&
+                !simulation::is_initial_condition)
               sigma_t = sigma_t_td_[material * negroups_ + g];
             else
               sigma_t = sigma_t_[material * negroups_ + g];
@@ -782,7 +783,8 @@ void FlatSourceDomain::random_ray_tally()
         case SCORE_FISSION:
           if (material != MATERIAL_VOID) {
             double sigma_f;
-            if (settings::kinetic_simulation && !settings::is_initial_condition)
+            if (settings::kinetic_simulation &&
+                !simulation::is_initial_condition)
               sigma_f = sigma_f_td_[material * negroups_ + g];
             else
               sigma_f = sigma_f_[material * negroups_ + g];
@@ -793,7 +795,8 @@ void FlatSourceDomain::random_ray_tally()
         case SCORE_NU_FISSION:
           if (material != MATERIAL_VOID) {
             double nu_sigma_f;
-            if (settings::kinetic_simulation && !settings::is_initial_condition)
+            if (settings::kinetic_simulation &&
+                !simulation::is_initial_condition)
               nu_sigma_f = nu_sigma_f_td_[material * negroups_ + g];
             else
               nu_sigma_f = nu_sigma_f_[material * negroups_ + g];
@@ -1804,7 +1807,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
 
   // Compute the combined source term
   update_single_neutron_source(handle);
-  if (settings::kinetic_simulation && !settings::is_initial_condition) {
+  if (settings::kinetic_simulation && !simulation::is_initial_condition) {
     update_single_neutron_source_td(handle);
     if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
       compute_single_neutron_source_time_derivative(handle);
@@ -1907,7 +1910,7 @@ void FlatSourceDomain::apply_transport_stabilization()
       }
       // TODO: test this and td with mesh
       // Duplicated stabilization for time-dependent simulations
-      if (settings::kinetic_simulation && !settings::is_initial_condition) {
+      if (settings::kinetic_simulation && !simulation::is_initial_condition) {
         double sigma_s_td =
           sigma_s_td_[material * negroups_ * negroups_ + g * negroups_ + g];
         if (sigma_s_td < 0.0) {
@@ -2110,7 +2113,7 @@ void FlatSourceDomain::compute_single_delayed_fission_source(
       if (lambda != 0.0) {
         for (int g = 0; g < negroups_; g++) {
           double scalar_flux;
-          if (settings::is_initial_condition) {
+          if (simulation::is_initial_condition) {
             scalar_flux = srh.scalar_flux_new(g);
           } else {
             scalar_flux = srh.scalar_flux_td_new(g);
@@ -2138,7 +2141,7 @@ void FlatSourceDomain::compute_single_precursors(SourceRegionHandle& srh)
       double lambda = lambda_[material * ndgroups_ + dg];
       if (lambda != 0.0) {
         double delayed_fission_source = srh.delayed_fission_source(dg);
-        if (settings::is_initial_condition) {
+        if (simulation::is_initial_condition) {
           srh.precursors_new(dg) = delayed_fission_source / lambda;
         } else {
           double precursor_rhs_bd = srh.precursors_rhs_bd(dg);
@@ -2220,7 +2223,7 @@ void FlatSourceDomain::accumulate_iteration_quantities()
         source_regions_.scalar_flux_td_final(sr, g) +=
           source_regions_.scalar_flux_td_new(sr, g);
         if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
-          if (settings::is_initial_condition)
+          if (simulation::is_initial_condition)
             source_regions_.source_final(sr, g) +=
               source_regions_.source(sr, g);
           else
@@ -2244,7 +2247,7 @@ void FlatSourceDomain::normalize_final_quantities()
   double source_normalization_factor;
   if (!settings::kinetic_simulation ||
       settings::kinetic_simulation &&
-        settings::current_timestep == settings::n_timesteps)
+        simulation::current_timestep == settings::n_timesteps)
     source_normalization_factor =
       compute_fixed_source_normalization_factor() * normalization_factor;
   else
@@ -2316,7 +2319,7 @@ void FlatSourceDomain::store_time_step_quantities(bool increment_not_initialize)
       if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
         // Multiply out sigma_t to store the base source
         double sigma_t;
-        if (settings::is_initial_condition) {
+        if (simulation::is_initial_condition) {
           sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g];
         } else {
           sigma_t = sigma_t_td_[source_regions_.material(sr) * negroups_ + g];

src/random_ray/random_ray.cpp

@@ -456,7 +456,7 @@ void RandomRay::attenuate_flux_flat_source(
     float new_delta_psi = (angular_flux_[g] - srh.source(g)) * exponential;
     delta_psi_[g] = new_delta_psi;
     angular_flux_[g] -= new_delta_psi;
-    if (settings::kinetic_simulation && !settings::is_initial_condition) {
+    if (settings::kinetic_simulation && !simulation::is_initial_condition) {
       float sigma_t_td = domain_->sigma_t_td_[material * negroups_ + g];
       float tau_td = sigma_t_td * distance;
       float exponential_td =
@@ -499,7 +499,7 @@ void RandomRay::attenuate_flux_flat_source(
     // this iteration
     for (int g = 0; g < negroups_; g++) {
       srh.scalar_flux_new(g) += delta_psi_[g];
-      if (settings::kinetic_simulation && !settings::is_initial_condition)
+      if (settings::kinetic_simulation && !simulation::is_initial_condition)
         srh.scalar_flux_td_new(g) += delta_psi_td_[g];
     }
 
@@ -860,7 +860,7 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
     for (int g = 0; g < negroups_; g++) {
       angular_flux_[g] = srh.source(g);
     }
-    if (settings::kinetic_simulation && !settings::is_initial_condition) {
+    if (settings::kinetic_simulation && !simulation::is_initial_condition) {
       for (int g = 0; g < negroups_; g++) {
         angular_flux_td_[g] = srh.source_td(g);
       }

src/random_ray/random_ray_simulation.cpp

@@ -35,11 +35,16 @@ void openmc_run_random_ray()
   // TODO: no initial condition needed fof td fixed source simulations
   // Check if this simulation is to establish an initial condition
   if (settings::kinetic_simulation)
-    settings::is_initial_condition = true;
+    simulation::is_initial_condition = true;
 
   // Configure the domain for forward simulation
   FlatSourceDomain::adjoint_ = false;
 
+  // If we're going to do a kinetic simulation, report that this is
+  // the initial condition.
+  if (settings::kinetic_simulation && mpi::master)
+    header("KINETIC SIMULATION INITIAL CONDITION", 3);
+
   // If we're going to do an adjoint simulation afterwards, report that this is
   // the initial forward flux solve.
   if (adjoint_needed && mpi::master)
@@ -91,16 +96,21 @@ void openmc_run_random_ray()
     // Timestepping loop
     // TODO: Add support for time-dependent restart
     for (int i = 0; i < settings::n_timesteps; i++) {
-      settings::current_timestep = i + 1;
+      simulation::current_timestep = i + 1;
 
       // Print simulation information
       if (mpi::master) {
         std::string message = fmt::format(
-          "KINETIC SIMULATION TIME STEP {0}", settings::current_timestep);
+          "KINETIC SIMULATION TIME STEP {0}", simulation::current_timestep);
         const char* msg = message.c_str();
         header(msg, 3);
       }
 
+      // If we're going to do an adjoint simulation afterwards, report that this
+      // is the initial forward flux solve.
+      if (adjoint_needed && mpi::master)
+        header("FORWARD FLUX SOLVE", 3);
+
       reset_timers();
 
       // Initialize OpenMC general data structures
@@ -133,9 +143,9 @@ void openmc_run_random_ray()
       sim.output_simulation_results();
 
       // Rename statepoint and tallies file
-      rename_statepoint_file(settings::current_timestep);
+      rename_statepoint_file(simulation::current_timestep);
       if (settings::output_tallies) {
-        rename_tallies_file(settings::current_timestep);
+        rename_tallies_file(simulation::current_timestep);
       }
 
       // Normalize and store final quantities for next time step
@@ -493,7 +503,7 @@ void write_random_ray_hdf5(hid_t group)
     RandomRay::total_geometric_intersections_ * data::mg.num_energy_groups_;
   write_dataset(random_ray_group, "n_integrations", n_integrations);
 
-  if (settings::kinetic_simulation && !settings::is_initial_condition) {
+  if (settings::kinetic_simulation && !simulation::is_initial_condition) {
     write_dataset(random_ray_group, "bd_order", RandomRay::bd_order_);
     switch (RandomRay::time_method_) {
     case RandomRayTimeMethod::ISOTROPIC:
@@ -515,7 +525,7 @@ void write_random_ray_hdf5(hid_t group)
 void set_time_dependent_settings()
 {
   // Reset flags
-  settings::is_initial_condition = false;
+  simulation::is_initial_condition = false;
 
   // Set current time
   simulation::current_time = settings::dt;
@@ -619,7 +629,7 @@ void RandomRaySimulation::simulate()
       // domain_->compute_neutron_source()
       // Update source term (scattering + fission)
       domain_->update_all_neutron_sources();
-      if (settings::kinetic_simulation && !settings::is_initial_condition) {
+      if (settings::kinetic_simulation && !simulation::is_initial_condition) {
         domain_->update_all_neutron_sources_td();
       }
 
@@ -688,7 +698,7 @@ void RandomRaySimulation::simulate()
 
       // Set phi_old = phi_new
       domain_->flux_swap();
-      if (settings::kinetic_simulation && !settings::is_initial_condition) {
+      if (settings::kinetic_simulation && !simulation::is_initial_condition) {
         domain_->flux_td_swap();
         domain_->precursors_swap();
       }
@@ -766,7 +776,7 @@ void RandomRaySimulation::print_results_random_ray() const
                        time_transport.elapsed() - time_tallies.elapsed() -
                        time_bank_sendrecv.elapsed();
 
-    if (settings::kinetic_simulation && !settings::is_initial_condition) {
+    if (settings::kinetic_simulation && !simulation::is_initial_condition) {
       misc_time -= time_update_bd_vectors_td.elapsed();
     }
     header("Simulation Statistics", 4);
@@ -846,7 +856,7 @@ void RandomRaySimulation::print_results_random_ray() const
     } else {
       fmt::print(" Transport XS Stabilization Used   = NO\n");
     }
-    if (settings::kinetic_simulation && !settings::is_initial_condition) {
+    if (settings::kinetic_simulation && !simulation::is_initial_condition) {
       std::string time_method =
         (RandomRay::time_method_ == RandomRayTimeMethod::ISOTROPIC)
           ? "ISOTROPIC"
@@ -867,7 +877,7 @@ void RandomRaySimulation::print_results_random_ray() const
         "Precursor computation only", time_compute_precursors.elapsed(), 1);
       misc_time -= time_compute_precursors.elapsed();
 
-      if (!settings::is_initial_condition) {
+      if (!simulation::is_initial_condition) {
         show_time(
           "Time-dependent source update only", time_update_src_td.elapsed(), 1);
         misc_time -= time_update_src_td.elapsed();

src/random_ray/source_region.cpp

@@ -412,7 +412,7 @@ void SourceRegionContainer::adjoint_reset()
   std::fill(flux_moments_t_.begin(), flux_moments_t_.end(),
     MomentArray {0.0, 0.0, 0.0});
   // Reset arrays for kinetic adjoint simulations
-  if (settings::kinetic_simulation && !settings::is_initial_condition) {
+  if (settings::kinetic_simulation && !simulation::is_initial_condition) {
     std::fill(scalar_flux_td_old_.begin(), scalar_flux_td_old_.end(), 0.0);
     std::fill(scalar_flux_td_new_.begin(), scalar_flux_td_new_.end(), 0.0);
     std::fill(

src/settings.cpp

@@ -153,8 +153,6 @@ double weight_survive {1.0};
 // Timestep variables for kinetic simulation
 int n_timesteps;
 double dt;
-int current_timestep;
-bool is_initial_condition {false};
 
 } // namespace settings
 

src/simulation.cpp

@@ -323,6 +323,8 @@ const RegularMesh* ufs_mesh {nullptr};
 vector<double> k_generation;
 vector<int64_t> work_index;
 
+bool is_initial_condition {false};
+int current_timestep;
 double current_time;
 
 } // namespace simulation

src/state_point.cpp

@@ -132,7 +132,7 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
       hid_t timestep_group = create_group(file_id, "timestep_data");
       write_dataset(timestep_group, "dt", settings::dt);
       write_dataset(
-        timestep_group, "current_timestep", settings::current_timestep);
+        timestep_group, "current_timestep", simulation::current_timestep);
       write_dataset(timestep_group, "current_time", simulation::current_time);
       close_group(timestep_group);
     }