From 5bce5adabc42c9ebcfd3d5b9ef474b5eb9fc29aa Mon Sep 17 00:00:00 2001
From: Kevin Sawatzky <66632997+nuclearkevin@users.noreply.github.com>
Date: Wed, 14 Jan 2026 12:02:17 -0600
Subject: [PATCH] Support cell densities in the random ray solver (#3720)

---
 include/openmc/material.h                     |   8 +-
 include/openmc/random_ray/source_region.h     |  12 +-
 src/random_ray/flat_source_domain.cpp         |  53 ++--
 src/random_ray/linear_source_domain.cpp       |  11 +-
 src/random_ray/random_ray.cpp                 |   6 +-
 src/random_ray/source_region.cpp              |  12 +-
 .../random_ray_cell_density/__init__.py       |   0
 .../eigen/inputs_true.dat                     | 109 ++++++++
 .../eigen/results_true.dat                    | 171 ++++++++++++
 .../fs/inputs_true.dat                        | 244 ++++++++++++++++++
 .../fs/results_true.dat                       |   9 +
 .../random_ray_cell_density/test.py           |  43 +++
 12 files changed, 648 insertions(+), 30 deletions(-)
 create mode 100644 tests/regression_tests/random_ray_cell_density/__init__.py
 create mode 100644 tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat
 create mode 100644 tests/regression_tests/random_ray_cell_density/eigen/results_true.dat
 create mode 100644 tests/regression_tests/random_ray_cell_density/fs/inputs_true.dat
 create mode 100644 tests/regression_tests/random_ray_cell_density/fs/results_true.dat
 create mode 100644 tests/regression_tests/random_ray_cell_density/test.py

diff --git a/include/openmc/material.h b/include/openmc/material.h
index e36946c71..c10f25551 100644
--- a/include/openmc/material.h
+++ b/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 @@ public:
   //! \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]
diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h
index 0f5a747ff..c20d46abe 100644
--- a/include/openmc/random_ray/source_region.h
+++ b/include/openmc/random_ray/source_region.h
@@ -146,6 +146,7 @@ public:
 
   // Scalar fields
   int* material_;
+  double* density_mult_;
   int* is_small_;
   int* n_hits_;
   int* birthday_;
@@ -195,6 +196,9 @@ public:
   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 @@ public:
   //---------------------------------------
   // 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 @@ public:
   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 @@ private:
 
   // 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_;
diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp
index ec14795dd..2f5007fc9 100644
--- a/src/random_ray/flat_source_domain.cpp
+++ b/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,6 +1615,7 @@ 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;
     }
@@ -1605,9 +1623,10 @@ void FlatSourceDomain::apply_transport_stabilization()
       // 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);
 
diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp
index 47ffbb727..02f4c9e23 100644
--- a/src/random_ray/linear_source_domain.cpp
+++ b/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
diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp
index c19d136a4..1b61d8c20 100644
--- a/src/random_ray/random_ray.cpp
+++ b/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.
diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp
index 3b06f0ed0..15c65221a 100644
--- a/src/random_ray/source_region.cpp
+++ b/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();
diff --git a/tests/regression_tests/random_ray_cell_density/__init__.py b/tests/regression_tests/random_ray_cell_density/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat b/tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat
new file mode 100644
index 000000000..0dd354cf0
--- /dev/null
+++ b/tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat
@@ -0,0 +1,109 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="UO2 fuel">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="UO2"/>
+    </material>
+    <material id="2" name="Water">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="LWTR"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" name="fuel inner a" material="1" region="-2" density="2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0" universe="1"/>
+    <cell id="2" name="fuel inner b" material="1" region="2 -3" density="2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0" universe="1"/>
+    <cell id="3" name="fuel inner c" material="1" region="3 -1" density="2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0" universe="1"/>
+    <cell id="4" name="moderator inner a" material="2" region="1 -4" universe="1"/>
+    <cell id="5" name="moderator outer b" material="2" region="4 -5" universe="1"/>
+    <cell id="6" name="moderator outer c" material="2" region="5" universe="1"/>
+    <cell id="7" name="azimuthal_cell_0" fill="1" region="6 -7" universe="2"/>
+    <cell id="8" name="azimuthal_cell_1" fill="1" region="7 -8" universe="2"/>
+    <cell id="9" name="azimuthal_cell_2" fill="1" region="8 -9" universe="2"/>
+    <cell id="10" name="azimuthal_cell_3" fill="1" region="9 -10" universe="2"/>
+    <cell id="11" name="azimuthal_cell_4" fill="1" region="10 -11" universe="2"/>
+    <cell id="12" name="azimuthal_cell_5" fill="1" region="11 -12" universe="2"/>
+    <cell id="13" name="azimuthal_cell_6" fill="1" region="12 -13" universe="2"/>
+    <cell id="14" name="azimuthal_cell_7" fill="1" region="13 -6" universe="2"/>
+    <cell id="15" name="moderator infinite" material="2" universe="3"/>
+    <cell id="16" fill="4" universe="5"/>
+    <cell id="17" name="assembly" fill="6" region="14 -15 16 -17" universe="7"/>
+    <lattice id="4">
+      <pitch>0.126 0.126</pitch>
+      <dimension>10 10</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <universes>
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 </universes>
+    </lattice>
+    <lattice id="6">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+2 2 
+2 5 </universes>
+    </lattice>
+    <surface id="1" name="Fuel OR" type="z-cylinder" coeffs="0.0 0.0 0.54"/>
+    <surface id="2" name="inner ring a" type="z-cylinder" coeffs="0.0 0.0 0.33"/>
+    <surface id="3" name="inner ring b" type="z-cylinder" coeffs="0.0 0.0 0.45"/>
+    <surface id="4" name="outer ring a" type="z-cylinder" coeffs="0.0 0.0 0.6"/>
+    <surface id="5" name="outer ring b" type="z-cylinder" coeffs="0.0 0.0 0.69"/>
+    <surface id="6" type="plane" coeffs="-0.0 1.0 0 0"/>
+    <surface id="7" type="plane" coeffs="-0.7071067811865475 0.7071067811865476 0 0"/>
+    <surface id="8" type="plane" coeffs="-1.0 6.123233995736766e-17 0 0"/>
+    <surface id="9" type="plane" coeffs="-0.7071067811865476 -0.7071067811865475 0 0"/>
+    <surface id="10" type="plane" coeffs="-1.2246467991473532e-16 -1.0 0 0"/>
+    <surface id="11" type="plane" coeffs="0.7071067811865475 -0.7071067811865477 0 0"/>
+    <surface id="12" type="plane" coeffs="1.0 -1.8369701987210297e-16 0 0"/>
+    <surface id="13" type="plane" coeffs="0.7071067811865477 0.7071067811865474 0 0"/>
+    <surface id="14" name="minimum x" type="x-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="15" name="maximum x" type="x-plane" boundary="reflective" coeffs="1.26"/>
+    <surface id="16" name="minimum y" type="y-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="17" name="maximum y" type="y-plane" boundary="reflective" coeffs="1.26"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>100.0</distance_active>
+      <distance_inactive>20.0</distance_inactive>
+      <source type="independent" strength="1.0" particle="neutron">
+        <space type="box">
+          <parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
+        </space>
+      </source>
+      <volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
+    </random_ray>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="energy">
+      <bins>1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0</bins>
+    </filter>
+    <tally id="1" name="Mesh tally">
+      <filters>1 2</filters>
+      <scores>flux fission nu-fission</scores>
+      <estimator>analog</estimator>
+    </tally>
+  </tallies>
+</model>
diff --git a/tests/regression_tests/random_ray_cell_density/eigen/results_true.dat b/tests/regression_tests/random_ray_cell_density/eigen/results_true.dat
new file mode 100644
index 000000000..3e6ce4039
--- /dev/null
+++ b/tests/regression_tests/random_ray_cell_density/eigen/results_true.dat
@@ -0,0 +1,171 @@
+k-combined:
+7.606488E-01 9.025978E-03
+tally 1:
+6.209503E-01
+7.732732E-02
+4.335729E-01
+3.769026E-02
+1.055230E+00
+2.232538E-01
+4.077492E-01
+3.335434E-02
+1.188353E-01
+2.833644E-03
+2.892214E-01
+1.678476E-02
+2.715245E-01
+1.488194E-02
+1.736084E-02
+6.079207E-05
+4.225281E-02
+3.600946E-04
+3.700491E-01
+2.796457E-02
+2.369715E-02
+1.145314E-04
+5.767413E-02
+6.784132E-04
+1.223083E+00
+3.058118E-01
+2.790791E-02
+1.593458E-04
+6.792309E-02
+9.438893E-04
+4.095779E+00
+3.377358E+00
+1.251201E-02
+3.153331E-05
+3.096010E-02
+1.930723E-04
+2.608638E+00
+1.361111E+00
+7.059742E-02
+9.968342E-04
+1.963632E-01
+7.711968E-03
+1.164050E+00
+2.710289E-01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.325486E-01
+5.675375E-02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+3.058117E-01
+1.905099E-02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+4.645650E-01
+4.428410E-02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+1.350036E+00
+3.718611E-01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+3.538916E+00
+2.521670E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.201812E+00
+9.698678E-01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+6.776093E-01
+9.211106E-02
+2.527022E-01
+1.280619E-02
+6.150266E-01
+7.585593E-02
+4.194143E-01
+3.528585E-02
+6.303447E-02
+7.972233E-04
+1.534133E-01
+4.722259E-03
+2.743408E-01
+1.520611E-02
+8.961769E-03
+1.621446E-05
+2.181115E-02
+9.604444E-05
+3.871710E-01
+3.063162E-02
+1.288801E-02
+3.392550E-05
+3.136683E-02
+2.009537E-04
+1.259769E+00
+3.241545E-01
+1.478360E-02
+4.466102E-05
+3.598077E-02
+2.645508E-04
+4.032169E+00
+3.273866E+00
+6.215187E-03
+7.776862E-06
+1.537904E-02
+4.761620E-05
+2.583647E+00
+1.335275E+00
+3.533526E-02
+2.498005E-04
+9.828323E-02
+1.932572E-03
+7.851145E-01
+1.234783E-01
+2.966638E-01
+1.763751E-02
+7.220204E-01
+1.044737E-01
+4.505673E-01
+4.067697E-02
+6.823342E-02
+9.332472E-04
+1.660665E-01
+5.527980E-03
+2.832534E-01
+1.626181E-02
+9.297014E-03
+1.749663E-05
+2.262707E-02
+1.036392E-04
+4.056699E-01
+3.370803E-02
+1.358439E-02
+3.774180E-05
+3.306168E-02
+2.235591E-04
+1.279744E+00
+3.345476E-01
+1.512005E-02
+4.668112E-05
+3.679962E-02
+2.765169E-04
+3.917888E+00
+3.090981E+00
+6.082371E-03
+7.449625E-06
+1.505040E-02
+4.561259E-05
+2.520906E+00
+1.271106E+00
+3.477796E-02
+2.419730E-04
+9.673314E-02
+1.872015E-03
diff --git a/tests/regression_tests/random_ray_cell_density/fs/inputs_true.dat b/tests/regression_tests/random_ray_cell_density/fs/inputs_true.dat
new file mode 100644
index 000000000..e90f25973
--- /dev/null
+++ b/tests/regression_tests/random_ray_cell_density/fs/inputs_true.dat
@@ -0,0 +1,244 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="source">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="source"/>
+    </material>
+    <material id="2" name="void">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="void"/>
+    </material>
+    <material id="3" name="absorber">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="absorber"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" name="infinite source region" material="1" density="1000.0" universe="1"/>
+    <cell id="2" name="infinite void region" material="2" universe="2"/>
+    <cell id="3" name="infinite absorber region" material="3" universe="3"/>
+    <cell id="4" fill="4" universe="5"/>
+    <cell id="5" name="full domain" fill="5" region="1 -2 3 -4 5 -6" universe="6"/>
+    <lattice id="4">
+      <pitch>2.5 2.5 2.5</pitch>
+      <dimension>12 12 12</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <universes>
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 </universes>
+    </lattice>
+    <surface id="1" type="x-plane" boundary="reflective" coeffs="0.0"/>
+    <surface id="2" type="x-plane" boundary="vacuum" coeffs="30.0"/>
+    <surface id="3" type="y-plane" boundary="reflective" coeffs="0.0"/>
+    <surface id="4" type="y-plane" boundary="vacuum" coeffs="30.0"/>
+    <surface id="5" type="z-plane" boundary="reflective" coeffs="0.0"/>
+    <surface id="6" type="z-plane" boundary="vacuum" coeffs="30.0"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>90</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source type="independent" strength="3.14" particle="neutron">
+      <energy type="discrete">
+        <parameters>100.0 1.0</parameters>
+      </energy>
+      <constraints>
+        <domain_type>universe</domain_type>
+        <domain_ids>1</domain_ids>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>500.0</distance_active>
+      <distance_inactive>100.0</distance_inactive>
+      <source type="independent" strength="1.0" particle="neutron">
+        <space type="box">
+          <parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
+        </space>
+      </source>
+      <volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
+    </random_ray>
+  </settings>
+  <tallies>
+    <filter id="3" type="material">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="material">
+      <bins>2</bins>
+    </filter>
+    <filter id="1" type="material">
+      <bins>3</bins>
+    </filter>
+    <tally id="3" name="Source Tally">
+      <filters>3</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="2" name="Void Tally">
+      <filters>2</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="1" name="Absorber Tally">
+      <filters>1</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+  </tallies>
+</model>
diff --git a/tests/regression_tests/random_ray_cell_density/fs/results_true.dat b/tests/regression_tests/random_ray_cell_density/fs/results_true.dat
new file mode 100644
index 000000000..e12f476e0
--- /dev/null
+++ b/tests/regression_tests/random_ray_cell_density/fs/results_true.dat
@@ -0,0 +1,9 @@
+tally 1:
+6.203077E-01
+7.706659E-02
+tally 2:
+3.203448E-02
+2.058679E-04
+tally 3:
+2.091970E-03
+8.765168E-07
diff --git a/tests/regression_tests/random_ray_cell_density/test.py b/tests/regression_tests/random_ray_cell_density/test.py
new file mode 100644
index 000000000..cb6062cdc
--- /dev/null
+++ b/tests/regression_tests/random_ray_cell_density/test.py
@@ -0,0 +1,43 @@
+import os
+
+import openmc
+from openmc.examples import random_ray_lattice, random_ray_three_region_cube
+from openmc.utility_funcs import change_directory
+import pytest
+
+from tests.testing_harness import TolerantPyAPITestHarness
+
+
+class MGXSTestHarness(TolerantPyAPITestHarness):
+    def _cleanup(self):
+        super()._cleanup()
+        f = 'mgxs.h5'
+        if os.path.exists(f):
+            os.remove(f)
+
+
+@pytest.mark.parametrize("run_mode", ["eigen", "fs"])
+def test_random_ray_basic(run_mode):
+    with change_directory(run_mode):
+        if run_mode == "eigen":
+            openmc.reset_auto_ids()
+            model = random_ray_lattice()
+            # Double the densities of the lower-left fuel pin -> cell instances [0, 9).
+            for id, cell in model.geometry.get_all_cells().items():
+                if cell.fill.name == "UO2 fuel":
+                    cell.density = [((i < 8) + 1.0) for i in range(24)]
+
+            # Gold file was generated with manually scaled fuel cross sections.
+            harness = MGXSTestHarness('statepoint.10.h5', model)
+            harness.main()
+        else:
+            openmc.reset_auto_ids()
+            model = random_ray_three_region_cube()
+            # Increase the density in the source region.
+            for id, cell in model.geometry.get_all_cells().items():
+                if cell.fill.name == "source":
+                    cell.density = 1e3
+
+            # Gold file was generated with manually scaled source cross sections.
+            harness = MGXSTestHarness('statepoint.10.h5', model)
+            harness.main()