Merge pull request #16 from yardasol/kin-rr-sim-simplification

Consistency changes, explicitly prohibit void regions in kinetic rr until it is fully supported.

Author: yardasol

include/openmc/random_ray/flat_source_domain.h

@@ -76,6 +76,7 @@ class FlatSourceDomain {
 
   //----------------------------------------------------------------------------
   // Methods for kinetic simulations
+  void compute_single_phi_prime(SourceRegionHandle& srh);
   void compute_single_T1(SourceRegionHandle& srh);
 
   void compute_single_delayed_fission_source(SourceRegionHandle& srh);

include/openmc/random_ray/source_region.h

@@ -197,6 +197,7 @@ class SourceRegionHandle {
   // Public Data Members for kinetic simulations
 
   // Energy group-wise 1D time-derivative arrays
+  double* phi_prime_;
   double* T1_;
 
   // Delay group-wise 1D arrays
@@ -342,6 +343,9 @@ class SourceRegionHandle {
 
   //---------------------------------------------------------------------------
   // Public Accessors for kinetic simulations
+  double& phi_prime(int g) { return phi_prime_[g]; }
+  const double phi_prime(int g) const { return phi_prime_[g]; }
+
   double& T1(int g) { return T1_[g]; }
   const double T1(int g) const { return T1_[g]; }
 
@@ -490,7 +494,9 @@ class SourceRegion {
                                //!< active iterations (used for SDP)
 
   // Energy group-wise 1D derivative arrays
-  vector<double> T1_; //!< The combined sourcetime derivative and 2nd order
+  vector<double>
+    phi_prime_; //!< The 1st order scalar flux time derivative (used for TI)
+  vector<double> T1_; //!< The combined source time derivative and 2nd order
                       //!< scalar flux time derivative (used for SDP)
   // Delay group-wise 1D arrays
   vector<double> delayed_fission_source_; //!< The delayed fission source binned
@@ -782,6 +788,14 @@ class SourceRegionContainer {
 
   //---------------------------------------
   // For kinetic simulations
+  double& phi_prime(int64_t sr, int g) { return phi_prime_[index(sr, g)]; }
+  const double& phi_prime(int64_t sr, int g) const
+  {
+    return phi_prime_[index(sr, g)];
+  }
+  double& phi_prime(int64_t se) { return phi_prime_[se]; }
+  const double& phi_prime(int64_t se) const { return phi_prime_[se]; }
+
   double& T1(int64_t sr, int g) { return T1_[index(sr, g)]; }
   const double& T1(int64_t sr, int g) const { return T1_[index(sr, g)]; }
   double& T1(int64_t se) { return T1_[se]; }
@@ -1034,6 +1048,7 @@ class SourceRegionContainer {
   // Private Data Members for kinetic simulations
 
   // SoA energy group-wise 2D derivative arrays flattened to 1D
+  vector<double> phi_prime_;
   vector<double> T1_;
 
   // SoA delay group-wise 2D arrays flattened to 1D

src/random_ray/flat_source_domain.cpp

@@ -149,39 +149,24 @@ void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
       }
       total_source = (scatter_source + fission_source * inverse_k_eff);
 
-      if (settings::kinetic_simulation && !simulation::is_initial_condition) {
-        // Add delayed source for kinetic simulation if delayed neutrons are
-        // turned on
-        if (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 precursors = srh.precursors_old(dg);
-            delayed_source += chi_d * precursors * lambda;
-          }
-          total_source += delayed_source;
-        }
-        // Add derivative of scalar flux to source (only works for isotropic
-        // method)
-        if (RandomRay::time_method_ == RandomRayTimeMethod::ISOTROPIC) {
-          double inverse_vbar = inverse_vbar_[material * negroups_ + g_out];
-          double scalar_flux_rhs_bd = srh.scalar_flux_rhs_bd(g_out);
-          double A0 =
-            (bd_coefficients_first_order_.at(RandomRay::bd_order_))[0] /
-            settings::dt;
-          double scalar_flux = srh.scalar_flux_old(g_out);
-          double scalar_flux_time_derivative =
-            A0 * scalar_flux + scalar_flux_rhs_bd;
-          total_source -= scalar_flux_time_derivative * inverse_vbar;
+      // Add delayed source for kinetic simulation if delayed neutrons are
+      // turned on
+      if (settings::kinetic_simulation && !simulation::is_initial_condition &&
+          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 precursors = srh.precursors_old(dg);
+          delayed_source += chi_d * precursors * lambda;
         }
+        total_source += delayed_source;
       }
       srh.source(g_out) = total_source / sigma_t;
     }
   }
 
-  // TODO: Add control flow for k-eigenvalue forward-weighted adjoint
   // Add external source if in fixed source mode
   if (settings::run_mode == RunMode::FIXED_SOURCE) {
     for (int g = 0; g < negroups_; g++) {
@@ -200,9 +185,11 @@ void FlatSourceDomain::update_all_neutron_sources()
   for (int64_t sr = 0; sr < n_source_regions(); sr++) {
     SourceRegionHandle srh = source_regions_.get_source_region_handle(sr);
     update_single_neutron_source(srh);
-    if (settings::kinetic_simulation && !simulation::is_initial_condition &&
-        RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
-      compute_single_T1(srh);
+    if (settings::kinetic_simulation && !simulation::is_initial_condition) {
+      if (RandomRay::time_method_ == RandomRayTimeMethod::ISOTROPIC)
+        compute_single_phi_prime(srh);
+      else if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION)
+        compute_single_T1(srh);
     }
   }
 
@@ -243,7 +230,6 @@ void FlatSourceDomain::set_flux_to_flux_plus_source(
   int64_t sr, double volume, int g)
 {
   int material = source_regions_.material(sr);
-  // TODO: Implement support for time-dependent void transport
   if (material == MATERIAL_VOID) {
     source_regions_.scalar_flux_new(sr, g) /= volume;
     if (settings::run_mode == RunMode::FIXED_SOURCE) {
@@ -252,20 +238,21 @@ void FlatSourceDomain::set_flux_to_flux_plus_source(
         source_regions_.volume_sq(sr);
     }
   } else {
-    double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g];
+    double sigma_t = sigma_t_[material * 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 && !simulation::is_initial_condition &&
-        RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
-      double inverse_vbar =
-        inverse_vbar_[source_regions_.material(sr) * negroups_ + g];
-      double scalar_flux_rhs_bd = source_regions_.scalar_flux_rhs_bd(sr, g);
-      double A0 = (bd_coefficients_first_order_.at(RandomRay::bd_order_))[0] /
-                  settings::dt;
-      source_regions_.scalar_flux_new(sr, g) -=
-        scalar_flux_rhs_bd * inverse_vbar / sigma_t;
-      source_regions_.scalar_flux_new(sr, g) /= 1 + A0 * inverse_vbar / sigma_t;
-    }
+  }
+  if (settings::kinetic_simulation && !simulation::is_initial_condition) {
+    double inverse_vbar =
+      inverse_vbar_[source_regions_.material(sr) * negroups_ + g];
+    double scalar_flux_rhs_bd = source_regions_.scalar_flux_rhs_bd(sr, g);
+    double A0 =
+      (bd_coefficients_first_order_.at(RandomRay::bd_order_))[0] / settings::dt;
+    // TODO: Add support for expicit void regions
+    double sigma_t = sigma_t_[material * negroups_ + g];
+    source_regions_.scalar_flux_new(sr, g) -=
+      scalar_flux_rhs_bd * inverse_vbar / sigma_t;
+    source_regions_.scalar_flux_new(sr, g) /= 1 + A0 * inverse_vbar / sigma_t;
   }
 }
 
@@ -1736,6 +1723,11 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
     }
   }
 
+  if (settings::kinetic_simulation && material == MATERIAL_VOID) {
+    fatal_error("Explicit void treatment for kinetic simulations "
+                " is not currently supported.");
+  }
+
   handle.material() = material;
 
   // Store the mesh index (if any) assigned to this source region
@@ -1774,9 +1766,11 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
 
   // Compute the combined source term
   update_single_neutron_source(handle);
-  if (settings::kinetic_simulation && !simulation::is_initial_condition &&
-      RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
-    compute_single_T1(handle);
+  if (settings::kinetic_simulation && !simulation::is_initial_condition) {
+    if (RandomRay::time_method_ == RandomRayTimeMethod::ISOTROPIC)
+      compute_single_phi_prime(handle);
+    else if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION)
+      compute_single_T1(handle);
   }
 
   // Unlock the parallel map. Note: we may be tempted to release
@@ -1936,6 +1930,22 @@ int64_t FlatSourceDomain::lookup_mesh_bin(int64_t sr, Position r) const
 // kinetic simulations) sources in each source region based on the flux
 // estimate from the previous iteration.
 
+void FlatSourceDomain::compute_single_phi_prime(SourceRegionHandle& srh)
+{
+  double A0 =
+    (bd_coefficients_first_order_.at(RandomRay::bd_order_))[0] / settings::dt;
+  int material = srh.material();
+  for (int g = 0; g < negroups_; g++) {
+    double inverse_vbar = inverse_vbar_[material * negroups_ + g];
+    // TODO: add support for explicit void
+    double sigma_t = sigma_t_[material * negroups_ + g];
+
+    double scalar_flux_time_derivative =
+      A0 * srh.scalar_flux_old(g) + srh.scalar_flux_rhs_bd(g);
+    srh.phi_prime(g) = scalar_flux_time_derivative * inverse_vbar / sigma_t;
+  }
+}
+
 // T1 calculation
 void FlatSourceDomain::compute_single_T1(SourceRegionHandle& srh)
 {
@@ -1946,9 +1956,8 @@ void FlatSourceDomain::compute_single_T1(SourceRegionHandle& srh)
   int material = srh.material();
   for (int g = 0; g < negroups_; g++) {
     double inverse_vbar = inverse_vbar_[material * negroups_ + g];
-    double sigma_t = 1.0;
-    if (material != MATERIAL_VOID)
-      sigma_t = sigma_t_[material * negroups_ + g];
+    // TODO: add support for explicit void
+    double sigma_t = sigma_t_[material * negroups_ + g];
 
     // Multiply out sigma_t to correctly compute the derivative term
     float source_time_derivative =
@@ -2131,10 +2140,10 @@ void FlatSourceDomain::store_time_step_quantities(bool increment_not_initialize)
         source_regions_.scalar_flux_final(sr, g), increment_not_initialize,
         RandomRay::bd_order_ + j);
       if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
-        // TODO: add support for void regions
         // Multiply out sigma_t to store the base source
-        double sigma_t = sigma_t =
-          sigma_t_[source_regions_.material(sr) * negroups_ + g];
+        int material = source_regions_.material(sr);
+        // TODO: add support for explicit void regions
+        double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g];
         float source = source_regions_.source_final(sr, g) * sigma_t;
         add_value_to_bd_vector(source_regions_.source_bd(sr, g), source,
           increment_not_initialize, RandomRay::bd_order_);

src/random_ray/random_ray.cpp

@@ -451,18 +451,21 @@ void RandomRay::attenuate_flux_flat_source(
     float tau = sigma_t * distance;
     float exponential = cjosey_exponential(tau); // exponential = 1 - exp(-tau)
     float new_delta_psi = (angular_flux_[g] - srh.source(g)) * exponential;
-    if (settings::kinetic_simulation && !simulation::is_initial_condition &&
-        RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
-      float T1 = srh.T1(g);
-
-      // Source Derivative Propogation terms for Characteristic Equation
-      float inverse_vbar = domain_->inverse_vbar_[material * negroups_ + g];
-      new_delta_psi += T1 * inverse_vbar * exponential / sigma_t;
-      new_delta_psi += distance * inverse_vbar * (angular_flux_prime_[g] - T1) *
-                       (1 - exponential);
-
-      // Time Derivative Characteristic Equation
-      angular_flux_prime_[g] -= (angular_flux_prime_[g] - T1) * exponential;
+    if (settings::kinetic_simulation && !simulation::is_initial_condition) {
+      if (RandomRay::time_method_ == RandomRayTimeMethod::ISOTROPIC) {
+        new_delta_psi += srh.phi_prime(g) * exponential;
+      } else if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
+        // Source Derivative Propogation terms for Characteristic Equation
+        float inverse_vbar = domain_->inverse_vbar_[material * negroups_ + g];
+        float T1 = srh.T1(g);
+        float new_delta_psi_prime = (angular_flux_prime_[g] - T1);
+        new_delta_psi += T1 * inverse_vbar * exponential / sigma_t;
+        new_delta_psi +=
+          distance * inverse_vbar * new_delta_psi_prime * (1 - exponential);
+
+        // Time Derivative Characteristic Equation
+        angular_flux_prime_[g] -= new_delta_psi_prime * exponential;
+      }
     }
     delta_psi_[g] = new_delta_psi;
     angular_flux_[g] -= new_delta_psi;
@@ -501,7 +504,7 @@ void RandomRay::attenuate_flux_flat_source(
 }
 
 // Alternative flux attenuation function for true void regions.
-// TODO: Implement support for time dependent voids
+// TODO: Implement support for time-dependent voids
 void RandomRay::attenuate_flux_flat_source_void(
   SourceRegionHandle& srh, double distance, bool is_active, Position r)
 {
@@ -841,7 +844,6 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
     if (settings::kinetic_simulation && !simulation::is_initial_condition &&
         RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
       for (int g = 0; g < negroups_; g++) {
-        // T1
         angular_flux_prime_[g] = srh.T1(g);
       }
     }

src/random_ray/random_ray_simulation.cpp

@@ -53,10 +53,6 @@ void openmc_run_random_ray()
     // Second steady state simulation to correct the batchwise k-eff
     sim.kinetic_initial_condition();
 
-    warning(
-      "Time-dependent explicit void treatment has not yet been "
-      "implemented. Use caution when interpreting results from models with "
-      "voids, as they may contain large inaccuracies.");
     // Timestepping loop
     for (int i = 0; i < settings::n_timesteps; i++)
       sim.kinetic_single_time_step(i);

src/random_ray/source_region.cpp

@@ -31,7 +31,8 @@ SourceRegionHandle::SourceRegionHandle(SourceRegion& sr)
     flux_moments_old_(sr.flux_moments_old_.data()),
     flux_moments_new_(sr.flux_moments_new_.data()),
     flux_moments_t_(sr.flux_moments_t_.data()),
-    tally_task_(sr.tally_task_.data()), T1_(sr.T1_.data()),
+    tally_task_(sr.tally_task_.data()), phi_prime_(sr.phi_prime_.data()),
+    T1_(sr.T1_.data()),
     delayed_fission_source_(sr.delayed_fission_source_.data()),
     precursors_old_(sr.precursors_old_.data()),
     precursors_new_(sr.precursors_new_.data()),
@@ -75,8 +76,11 @@ SourceRegion::SourceRegion(int negroups, int ndgroups, bool is_linear)
     scalar_flux_bd_.resize(negroups);
     scalar_flux_rhs_bd_.resize(negroups);
 
-    // Source Derivative Propogation arrays
-    if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
+    if (RandomRay::time_method_ == RandomRayTimeMethod::ISOTROPIC) {
+      // Time Isotropic arrays
+      phi_prime_.assign(negroups, 0.0);
+    } else if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
+      // Source Derivative Propogation arrays
       source_final_.assign(negroups, 0.0);
 
       T1_.assign(negroups, 0.0);
@@ -160,8 +164,11 @@ void SourceRegionContainer::push_back(const SourceRegion& sr)
       scalar_flux_bd_.push_back(sr.scalar_flux_bd_[g]);
       scalar_flux_rhs_bd_.push_back(sr.scalar_flux_rhs_bd_[g]);
 
-      // Source Derivative Propogation arrays
-      if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
+      if (RandomRay::time_method_ == RandomRayTimeMethod::ISOTROPIC) {
+        // Time Isotropic arrays
+        phi_prime_.push_back(sr.phi_prime_[g]);
+      } else if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
+        // Source Derivative Propogation arrays
         source_final_.push_back(sr.source_final_[g]);
 
         T1_.push_back(sr.T1_[g]);
@@ -237,7 +244,9 @@ void SourceRegionContainer::assign(
     scalar_flux_bd_.clear();
     scalar_flux_rhs_bd_.clear();
 
-    if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
+    if (RandomRay::time_method_ == RandomRayTimeMethod::ISOTROPIC) {
+      phi_prime_.clear();
+    } else if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
       source_final_.clear();
 
       T1_.clear();
@@ -320,7 +329,9 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr)
     handle.scalar_flux_bd_ = &scalar_flux_bd(sr, 0);
     handle.scalar_flux_rhs_bd_ = &scalar_flux_rhs_bd(sr, 0);
 
-    if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
+    if (RandomRay::time_method_ == RandomRayTimeMethod::ISOTROPIC) {
+      handle.phi_prime_ = &phi_prime(sr, 0);
+    } else if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
       handle.source_final_ = &source_final(sr, 0);
 
       handle.T1_ = &T1(sr, 0);
@@ -394,7 +405,9 @@ void SourceRegionContainer::adjoint_reset()
     // BD Vectors
     std::fill(scalar_flux_rhs_bd_.begin(), scalar_flux_rhs_bd_.end(), 0.0);
 
-    if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
+    if (RandomRay::time_method_ == RandomRayTimeMethod::ISOTROPIC) {
+      std::fill(phi_prime_.begin(), phi_prime_.end(), 0.0);
+    } else if (RandomRay::time_method_ == RandomRayTimeMethod::PROPAGATION) {
       std::fill(T1_.begin(), T1_.end(), 0.0);
 
       std::fill(source_rhs_bd_.begin(), source_rhs_bd_.end(), 0.0);

tests/regression_tests/random_ray_k_eff_kinetic/isotropic/results_true_1.dat

@@ -1,12 +1,12 @@
 k-combined:
 1.325787E+00 5.102917E-04
 tally 1:
-2.630190E+03
-3.458951E+04
+2.630191E+03
+3.458954E+04
 1.079890E+02
-5.830986E+01
-2.651431E+02
-3.515147E+02
+5.830990E+01
+2.651432E+02
+3.515149E+02
 tally 2:
 1.077582E+00
 5.805916E-03