Support cell densities in the random ray solver (#3720)

Author: nuclearkevin

include/openmc/material.h

@@ -14,6 +14,7 @@
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/ncrystal_interface.h"
 #include "openmc/particle.h"
+#include "openmc/settings.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -110,9 +111,12 @@ class Material {
   //! \return Density in [atom/b-cm]
   double density() const { return density_; }
 
-  //! Get density in [g/cm^3]
+  //! Get density in [g/cm^3].
   //! \return Density in [g/cm^3]
-  double density_gpcc() const { return density_gpcc_; }
+  double density_gpcc() const
+  {
+    return settings::run_CE ? density_gpcc_ : density();
+  }
 
   //! Get charge density in [e/b-cm]
   //! \return Charge density in [e/b-cm]

include/openmc/random_ray/source_region.h

@@ -146,6 +146,7 @@ class SourceRegionHandle {
 
   // Scalar fields
   int* material_;
+  double* density_mult_;
   int* is_small_;
   int* n_hits_;
   int* birthday_;
@@ -195,6 +196,9 @@ class SourceRegionHandle {
   int& material() { return *material_; }
   const int material() const { return *material_; }
 
+  double& density_mult() { return *density_mult_; }
+  const double density_mult() const { return *density_mult_; }
+
   int& is_small() { return *is_small_; }
   const int is_small() const { return *is_small_; }
 
@@ -316,7 +320,9 @@ class SourceRegion {
   //---------------------------------------
   // Scalar fields
 
-  int material_ {0}; //!< Index in openmc::model::materials array
+  int material_ {0};          //!< Index in openmc::model::materials array
+  double density_mult_ {1.0}; //!< A density multiplier queried from the cell
+                              //!< corresponding to the source region.
   OpenMPMutex lock_;
   double volume_ {
     0.0}; //!< Volume (computed from the sum of ray crossing lengths)
@@ -394,6 +400,9 @@ class SourceRegionContainer {
   int& material(int64_t sr) { return material_[sr]; }
   const int material(int64_t sr) const { return material_[sr]; }
 
+  double& density_mult(int64_t sr) { return density_mult_[sr]; }
+  const double density_mult(int64_t sr) const { return density_mult_[sr]; }
+
   int& is_small(int64_t sr) { return is_small_[sr]; }
   const int is_small(int64_t sr) const { return is_small_[sr]; }
 
@@ -625,6 +634,7 @@ class SourceRegionContainer {
 
   // SoA storage for scalar fields (one item per source region)
   vector<int> material_;
+  vector<double> density_mult_;
   vector<int> is_small_;
   vector<int> n_hits_;
   vector<int> mesh_;

src/random_ray/flat_source_domain.cpp

@@ -109,18 +109,21 @@ void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
 
   // Add scattering + fission source
   int material = srh.material();
+  double density_mult = srh.density_mult();
   if (material != MATERIAL_VOID) {
     double inverse_k_eff = 1.0 / k_eff_;
     for (int g_out = 0; g_out < negroups_; g_out++) {
-      double sigma_t = sigma_t_[material * negroups_ + g_out];
+      double sigma_t = sigma_t_[material * negroups_ + g_out] * density_mult;
       double scatter_source = 0.0;
       double fission_source = 0.0;
 
       for (int g_in = 0; g_in < negroups_; g_in++) {
         double scalar_flux = srh.scalar_flux_old(g_in);
-        double sigma_s =
-          sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];
-        double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];
+        double sigma_s = sigma_s_[material * negroups_ * negroups_ +
+                                  g_out * negroups_ + g_in] *
+                         density_mult;
+        double nu_sigma_f =
+          nu_sigma_f_[material * negroups_ + g_in] * density_mult;
         double chi = chi_[material * negroups_ + g_out];
 
         scatter_source += sigma_s * scalar_flux;
@@ -198,7 +201,8 @@ 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_[source_regions_.material(sr) * negroups_ + g] *
+                     source_regions_.density_mult(sr);
     source_regions_.scalar_flux_new(sr, g) /= (sigma_t * volume);
     source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g);
   }
@@ -332,7 +336,8 @@ void FlatSourceDomain::compute_k_eff()
     double sr_fission_source_new = 0;
 
     for (int g = 0; g < negroups_; g++) {
-      double nu_sigma_f = nu_sigma_f_[material * negroups_ + g];
+      double nu_sigma_f = nu_sigma_f_[material * negroups_ + g] *
+                          source_regions_.density_mult(sr);
       sr_fission_source_old +=
         nu_sigma_f * source_regions_.scalar_flux_old(sr, g);
       sr_fission_source_new +=
@@ -562,7 +567,8 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const
       // to get the total source strength in the expected units.
       double sigma_t = 1.0;
       if (material != MATERIAL_VOID) {
-        sigma_t = sigma_t_[material * negroups_ + g];
+        sigma_t =
+          sigma_t_[material * negroups_ + g] * source_regions_.density_mult(sr);
       }
       simulation_external_source_strength +=
         source_regions_.external_source(sr, g) * sigma_t * volume;
@@ -618,7 +624,9 @@ void FlatSourceDomain::random_ray_tally()
     // source strength.
     double volume = source_regions_.volume(sr) * simulation_volume_;
 
-    double material = source_regions_.material(sr);
+    int material = source_regions_.material(sr);
+    double density_mult = source_regions_.density_mult(sr);
+
     for (int g = 0; g < negroups_; g++) {
       double flux =
         source_regions_.scalar_flux_new(sr, g) * source_normalization_factor;
@@ -634,19 +642,22 @@ void FlatSourceDomain::random_ray_tally()
 
         case SCORE_TOTAL:
           if (material != MATERIAL_VOID) {
-            score = flux * volume * sigma_t_[material * negroups_ + g];
+            score =
+              flux * volume * sigma_t_[material * negroups_ + g] * density_mult;
           }
           break;
 
         case SCORE_FISSION:
           if (material != MATERIAL_VOID) {
-            score = flux * volume * sigma_f_[material * negroups_ + g];
+            score =
+              flux * volume * sigma_f_[material * negroups_ + g] * density_mult;
           }
           break;
 
         case SCORE_NU_FISSION:
           if (material != MATERIAL_VOID) {
-            score = flux * volume * nu_sigma_f_[material * negroups_ + g];
+            score = flux * volume * nu_sigma_f_[material * negroups_ + g] *
+                    density_mult;
           }
           break;
 
@@ -913,7 +924,8 @@ void FlatSourceDomain::output_to_vtk() const
             for (int g = 0; g < negroups_; g++) {
               int64_t source_element = fsr * negroups_ + g;
               float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g);
-              double sigma_f = sigma_f_[mat * negroups_ + g];
+              double sigma_f = sigma_f_[mat * negroups_ + g] *
+                               source_regions_.density_mult(fsr);
               total_fission += sigma_f * flux;
             }
           }
@@ -934,7 +946,8 @@ void FlatSourceDomain::output_to_vtk() const
             // multiply it back to get the true external source.
             double sigma_t = 1.0;
             if (mat != MATERIAL_VOID) {
-              sigma_t = sigma_t_[mat * negroups_ + g];
+              sigma_t = sigma_t_[mat * negroups_ + g] *
+                        source_regions_.density_mult(fsr);
             }
             total_external += source_regions_.external_source(fsr, g) * sigma_t;
           }
@@ -1244,7 +1257,8 @@ void FlatSourceDomain::set_adjoint_sources()
       continue;
     }
     for (int g = 0; g < negroups_; g++) {
-      double sigma_t = sigma_t_[material * negroups_ + g];
+      double sigma_t =
+        sigma_t_[material * negroups_ + g] * source_regions_.density_mult(sr);
       source_regions_.external_source(sr, g) /= sigma_t;
     }
   }
@@ -1495,6 +1509,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
 
   handle.material() = material;
 
+  handle.density_mult() = cell.density_mult(gs.cell_instance());
+
   // Store the mesh index (if any) assigned to this source region
   handle.mesh() = mesh_idx;
 
@@ -1523,7 +1539,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
     // Divide external source term by sigma_t
     if (material != C_NONE) {
       for (int g = 0; g < negroups_; g++) {
-        double sigma_t = sigma_t_[material * negroups_ + g];
+        double sigma_t =
+          sigma_t_[material * negroups_ + g] * handle.density_mult();
         handle.external_source(g) /= sigma_t;
       }
     }
@@ -1598,16 +1615,18 @@ void FlatSourceDomain::apply_transport_stabilization()
 #pragma omp parallel for
   for (int64_t sr = 0; sr < n_source_regions(); sr++) {
     int material = source_regions_.material(sr);
+    double density_mult = source_regions_.density_mult(sr);
     if (material == MATERIAL_VOID) {
       continue;
     }
     for (int g = 0; g < negroups_; g++) {
       // Only apply stabilization if the diagonal (in-group) scattering XS is
       // negative
       double sigma_s =
-        sigma_s_[material * negroups_ * negroups_ + g * negroups_ + g];
+        sigma_s_[material * negroups_ * negroups_ + g * negroups_ + g] *
+        density_mult;
       if (sigma_s < 0.0) {
-        double sigma_t = sigma_t_[material * negroups_ + g];
+        double sigma_t = sigma_t_[material * negroups_ + g] * density_mult;
         double phi_new = source_regions_.scalar_flux_new(sr, g);
         double phi_old = source_regions_.scalar_flux_old(sr, g);
 

src/random_ray/linear_source_domain.cpp

@@ -43,12 +43,13 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
 
   // Add scattering + fission source
   int material = srh.material();
+  double density_mult = srh.density_mult();
   if (material != MATERIAL_VOID) {
     double inverse_k_eff = 1.0 / k_eff_;
     MomentMatrix invM = srh.mom_matrix().inverse();
 
     for (int g_out = 0; g_out < negroups_; g_out++) {
-      double sigma_t = sigma_t_[material * negroups_ + g_out];
+      double sigma_t = sigma_t_[material * negroups_ + g_out] * density_mult;
 
       double scatter_flat = 0.0f;
       double fission_flat = 0.0f;
@@ -61,9 +62,11 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
         MomentArray flux_linear = srh.flux_moments_old(g_in);
 
         // Handles for cross sections
-        double sigma_s =
-          sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];
-        double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];
+        double sigma_s = sigma_s_[material * negroups_ * negroups_ +
+                                  g_out * negroups_ + g_in] *
+                         density_mult;
+        double nu_sigma_f =
+          nu_sigma_f_[material * negroups_ + g_in] * density_mult;
         double chi = chi_[material * negroups_ + g_out];
 
         // Compute source terms for flat and linear components of the flux

src/random_ray/random_ray.cpp

@@ -435,7 +435,8 @@ void RandomRay::attenuate_flux_flat_source(
 
   // MOC incoming flux attenuation + source contribution/attenuation equation
   for (int g = 0; g < negroups_; g++) {
-    float sigma_t = domain_->sigma_t_[material * negroups_ + g];
+    float sigma_t =
+      domain_->sigma_t_[material * negroups_ + g] * srh.density_mult();
     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;
@@ -558,7 +559,8 @@ void RandomRay::attenuate_flux_linear_source(
   for (int g = 0; g < negroups_; g++) {
 
     // Compute tau, the optical thickness of the ray segment
-    float sigma_t = domain_->sigma_t_[material * negroups_ + g];
+    float sigma_t =
+      domain_->sigma_t_[material * negroups_ + g] * srh.density_mult();
     float tau = sigma_t * distance;
 
     // If tau is very small, set it to zero to avoid numerical issues.

src/random_ray/source_region.cpp

@@ -11,10 +11,11 @@ namespace openmc {
 //==============================================================================
 SourceRegionHandle::SourceRegionHandle(SourceRegion& sr)
   : negroups_(sr.scalar_flux_old_.size()), material_(&sr.material_),
-    is_small_(&sr.is_small_), n_hits_(&sr.n_hits_),
-    is_linear_(sr.source_gradients_.size() > 0), lock_(&sr.lock_),
-    volume_(&sr.volume_), volume_t_(&sr.volume_t_), volume_sq_(&sr.volume_sq_),
-    volume_sq_t_(&sr.volume_sq_t_), volume_naive_(&sr.volume_naive_),
+    density_mult_(&sr.density_mult_), is_small_(&sr.is_small_),
+    n_hits_(&sr.n_hits_), is_linear_(sr.source_gradients_.size() > 0),
+    lock_(&sr.lock_), volume_(&sr.volume_), volume_t_(&sr.volume_t_),
+    volume_sq_(&sr.volume_sq_), volume_sq_t_(&sr.volume_sq_t_),
+    volume_naive_(&sr.volume_naive_),
     position_recorded_(&sr.position_recorded_),
     external_source_present_(&sr.external_source_present_),
     position_(&sr.position_), centroid_(&sr.centroid_),
@@ -70,6 +71,7 @@ void SourceRegionContainer::push_back(const SourceRegion& sr)
 
   // Scalar fields
   material_.push_back(sr.material_);
+  density_mult_.push_back(sr.density_mult_);
   is_small_.push_back(sr.is_small_);
   n_hits_.push_back(sr.n_hits_);
   lock_.push_back(sr.lock_);
@@ -123,6 +125,7 @@ void SourceRegionContainer::assign(
   // Clear existing data
   n_source_regions_ = 0;
   material_.clear();
+  density_mult_.clear();
   is_small_.clear();
   n_hits_.clear();
   lock_.clear();
@@ -180,6 +183,7 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr)
   SourceRegionHandle handle;
   handle.negroups_ = negroups();
   handle.material_ = &material(sr);
+  handle.density_mult_ = &density_mult(sr);
   handle.is_small_ = &is_small(sr);
   handle.n_hits_ = &n_hits(sr);
   handle.is_linear_ = is_linear();