refactor is compiling cleanly

paulromano · Feb 11, 2025 · Jan 29, 2025 · Jan 29, 2025 · Jan 30, 2025 · Jan 30, 2025
commit 412928193e1f67326246985f1a7c4ee86c12b110
diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h
@@ -4,8 +4,8 @@
 #include "openmc/constants.h"
 #include "openmc/openmp_interface.h"
 #include "openmc/position.h"
-#include "openmc/source.h"
 #include "openmc/random_ray/source_region.h"
+#include "openmc/source.h"
 #include <unordered_set>
 
 namespace openmc {
@@ -46,6 +46,7 @@ class FlatSourceDomain {
   double compute_fixed_source_normalization_factor() const;
   void flatten_xs();
   void transpose_scattering_matrix();
+  void serialize_final_fluxes(vector<double>& flux);
 
   //----------------------------------------------------------------------------
   // Static Data members
@@ -58,7 +59,6 @@ class FlatSourceDomain {
 
   //----------------------------------------------------------------------------
   // Public Data members
-
   bool mapped_all_tallies_ {false}; // If all source regions have been visited
 
   int64_t n_source_regions_ {0}; // Total number of source regions in the model
@@ -100,8 +100,7 @@ class FlatSourceDomain {
     const vector<int32_t>& instances);
   void apply_external_source_to_cell_and_children(int32_t i_cell,
     Discrete* discrete, double strength_factor, int32_t target_material_id);
-  virtual void set_flux_to_flux_plus_source(
-    int sr, double volume, int g);
+  virtual void set_flux_to_flux_plus_source(int sr, double volume, int g);
   void set_flux_to_source(int sr, int g);
   virtual void set_flux_to_old_flux(int sr, int g);
 

diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h
@@ -4,7 +4,7 @@
 #include "openmc/openmp_interface.h"
 #include "openmc/position.h"
 #include "openmc/random_ray/moment_matrix.h"
-#include "openmc/random_ray/random_ray.h"
+#include "openmc/settings.h"
 
 namespace openmc {
 
@@ -85,7 +85,7 @@ class SourceRegion {
 public:
   //----------------------------------------------------------------------------
   // Constructors
-  SourceRegion(int negroups)
+  SourceRegion(int negroups, bool is_linear)
   {
     if (settings::run_mode == RunMode::EIGENVALUE) {
       // If in eigenvalue mode, set starting flux to guess of 1
@@ -103,7 +103,7 @@ class SourceRegion {
 
     tally_task_.resize(negroups);
 
-    if (RandomRay::source_shape_ != RandomRaySourceShape::FLAT) {
+    if (is_linear) {
       source_gradients_.resize(negroups);
       flux_moments_old_.resize(negroups);
       flux_moments_new_.resize(negroups);

diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp
@@ -47,7 +47,8 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
 
   // Initialize cell-wise arrays
   volume_task_.resize(n_source_regions_);
-  source_regions_.resize(n_source_regions_, SourceRegion(negroups_));
+  bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT;
+  source_regions_.resize(n_source_regions_, SourceRegion(negroups_, is_linear));
 
   // Initialize material array
   int64_t source_region_id = 0;
@@ -1123,4 +1124,17 @@ void FlatSourceDomain::transpose_scattering_matrix()
   }
 }
 
+void FlatSourceDomain::serialize_final_fluxes(vector<double>& flux)
+{
+  // Ensure array is correct size
+  flux.resize(n_source_regions_ * negroups_);
+  // Serialize the final fluxes for output
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    SourceRegion& region = source_regions_[sr];
+    for (int g = 0; g < negroups_; g++) {
+      flux[sr * negroups_ + g] = region.scalar_flux_final_[g];
+    }
+  }
+}
+
 } // namespace openmc
diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp
@@ -309,20 +309,18 @@ void RandomRay::attenuate_flux_flat_source(double distance, bool is_active)
   int i_cell = lowest_coord().cell;
 
   // The source region is the spatial region index
-  int64_t source_region =
-    domain_->source_region_offsets_[i_cell] + cell_instance();
+  int sr = domain_->source_region_offsets_[i_cell] + cell_instance();
+  SourceRegion& region = domain_->source_regions_[sr];
 
   // The source element is the energy-specific region index
-  int64_t source_element = source_region * negroups_;
   int material = this->material();
 
   // 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 tau = sigma_t * distance;
     float exponential = cjosey_exponential(tau); // exponential = 1 - exp(-tau)
-    float new_delta_psi =
-      (angular_flux_[g] - domain_->source_[source_element + g]) * exponential;
+    float new_delta_psi = (angular_flux_[g] - region.source_[g]) * exponential;
     delta_psi_[g] = new_delta_psi;
     angular_flux_[g] -= new_delta_psi;
   }
@@ -332,28 +330,28 @@ void RandomRay::attenuate_flux_flat_source(double distance, bool is_active)
   if (is_active) {
 
     // Aquire lock for source region
-    domain_->lock_[source_region].lock();
+    region.lock_.lock();
 
     // Accumulate delta psi into new estimate of source region flux for
     // this iteration
     for (int g = 0; g < negroups_; g++) {
-      domain_->scalar_flux_new_[source_element + g] += delta_psi_[g];
+      region.scalar_flux_new_[g] += delta_psi_[g];
     }
 
     // Accomulate volume (ray distance) into this iteration's estimate
     // of the source region's volume
-    domain_->volume_[source_region] += distance;
+    region.volume_ += distance;
 
     // Tally valid position inside the source region (e.g., midpoint of
     // the ray) if not done already
-    if (!domain_->position_recorded_[source_region]) {
+    if (!region.position_recorded_) {
       Position midpoint = r() + u() * (distance / 2.0);
-      domain_->position_[source_region] = midpoint;
-      domain_->position_recorded_[source_region] = 1;
+      region.position_ = midpoint;
+      region.position_recorded_ = 1;
     }
 
     // Release lock
-    domain_->lock_[source_region].unlock();
+    region.lock_.unlock();
   }
 }
 
@@ -373,14 +371,13 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
   int i_cell = lowest_coord().cell;
 
   // The source region is the spatial region index
-  int64_t source_region =
-    domain_->source_region_offsets_[i_cell] + cell_instance();
+  int sr = domain_->source_region_offsets_[i_cell] + cell_instance();
+  SourceRegion& region = domain_->source_regions_[sr];
 
   // The source element is the energy-specific region index
-  int64_t source_element = source_region * negroups_;
   int material = this->material();
 
-  Position& centroid = domain->centroid_[source_region];
+  Position& centroid = region.centroid_;
   Position midpoint = r() + u() * (distance / 2.0);
 
   // Determine the local position of the midpoint and the ray origin
@@ -393,7 +390,7 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
   // be no estimate of its centroid. We detect this by checking if it has
   // any accumulated volume. If its volume is zero, just use the midpoint
   // of the ray as the region's centroid.
-  if (domain->volume_t_[source_region]) {
+  if (region.volume_t_) {
     rm_local = midpoint - centroid;
     r0_local = r() - centroid;
   } else {
@@ -420,9 +417,8 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
     // calculated from the source gradients dot product with local centroid
     // and direction, respectively.
     float spatial_source =
-      domain_->source_[source_element + g] +
-      rm_local.dot(domain->source_gradients_[source_element + g]);
-    float dir_source = u().dot(domain->source_gradients_[source_element + g]);
+      region.source_[g] + rm_local.dot(region.source_gradients_[g]);
+    float dir_source = u().dot(region.source_gradients_[g]);
 
     float gn = exponentialG(tau);
     float f1 = 1.0f - tau * gn;
@@ -465,33 +461,33 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
       rm_local, u(), distance);
 
     // Aquire lock for source region
-    domain_->lock_[source_region].lock();
+    region.lock_.lock();
 
     // Accumulate deltas into the new estimate of source region flux for this
     // iteration
     for (int g = 0; g < negroups_; g++) {
-      domain_->scalar_flux_new_[source_element + g] += delta_psi_[g];
-      domain->flux_moments_new_[source_element + g] += delta_moments_[g];
+      region.scalar_flux_new_[g] += delta_psi_[g];
+      region.flux_moments_new_[g] += delta_moments_[g];
     }
 
     // Accumulate the volume (ray segment distance), centroid, and spatial
     // momement estimates into the running totals for the iteration for this
     // source region. The centroid and spatial momements estimates are scaled by
     // the ray segment length as part of length averaging of the estimates.
-    domain_->volume_[source_region] += distance;
-    domain->centroid_iteration_[source_region] += midpoint * distance;
+    region.volume_ += distance;
+    region.centroid_iteration_ += midpoint * distance;
     moment_matrix_estimate *= distance;
-    domain->mom_matrix_[source_region] += moment_matrix_estimate;
+    region.mom_matrix_ += moment_matrix_estimate;
 
     // Tally valid position inside the source region (e.g., midpoint of
     // the ray) if not done already
-    if (!domain_->position_recorded_[source_region]) {
-      domain_->position_[source_region] = midpoint;
-      domain_->position_recorded_[source_region] = 1;
+    if (!region.position_recorded_) {
+      region.position_ = midpoint;
+      region.position_recorded_ = 1;
     }
 
     // Release lock
-    domain_->lock_[source_region].unlock();
+    region.lock_.unlock();
   }
 }
 
@@ -537,11 +533,11 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
   // Initialize ray's starting angular flux to starting location's isotropic
   // source
   int i_cell = lowest_coord().cell;
-  int64_t source_region_idx =
-    domain_->source_region_offsets_[i_cell] + cell_instance();
+  int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
+  SourceRegion& region = domain_->source_regions_[sr];
 
   for (int g = 0; g < negroups_; g++) {
-    angular_flux_[g] = domain_->source_[source_region_idx * negroups_ + g];
+    angular_flux_[g] = region.source_[g];
   }
 }
 

diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp
@@ -66,7 +66,7 @@ void openmc_run_random_ray()
     simulation::time_total.stop();
 
     // Normalize and save the final forward flux
-    forward_flux = sim.domain()->scalar_flux_final_;
+    sim.domain()->serialize_final_fluxes(forward_flux);
 
     double source_normalization_factor =
       sim.domain()->compute_fixed_source_normalization_factor() /