From b22656e57f9469d25f14a203831bcdb390cd84db Mon Sep 17 00:00:00 2001
From: John Tramm <john.tramm@gmail.com>
Date: Thu, 15 Aug 2024 20:47:05 -0500
Subject: [PATCH] Fix random ray solver to correctly simulate fixed source
 problems with fissionable materials (#3106)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
---
 openmc/examples.py                            |   2 +-
 src/random_ray/flat_source_domain.cpp         |  40 ++---
 .../__init__.py                               |   0
 .../flat/inputs_true.dat                      | 140 +++++++++++++++
 .../flat/results_true.dat                     | 169 ++++++++++++++++++
 .../linear_xy/inputs_true.dat                 | 140 +++++++++++++++
 .../linear_xy/results_true.dat                | 169 ++++++++++++++++++
 .../test.py                                   | 133 ++++++++++++++
 8 files changed, 772 insertions(+), 21 deletions(-)
 create mode 100644 tests/regression_tests/random_ray_fixed_source_subcritical/__init__.py
 create mode 100644 tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat
 create mode 100644 tests/regression_tests/random_ray_fixed_source_subcritical/flat/results_true.dat
 create mode 100644 tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat
 create mode 100644 tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/results_true.dat
 create mode 100644 tests/regression_tests/random_ray_fixed_source_subcritical/test.py

diff --git a/openmc/examples.py b/openmc/examples.py
index 038c75ad2..538b6ea4a 100644
--- a/openmc/examples.py
+++ b/openmc/examples.py
@@ -804,7 +804,7 @@ def random_ray_lattice() -> openmc.Model:
         azimuthal_cells.append(azimuthal_cell)
 
     # Create a geometry with the azimuthal universes
-    pincell = openmc.Universe(cells=azimuthal_cells)
+    pincell = openmc.Universe(cells=azimuthal_cells, name='pincell')
 
     ########################################
     # Define a moderator lattice universe
diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp
index 8b6cda93f..584b3a7ed 100644
--- a/src/random_ray/flat_source_domain.cpp
+++ b/src/random_ray/flat_source_domain.cpp
@@ -158,31 +158,31 @@ void FlatSourceDomain::update_neutron_source(double k_eff)
     }
   }
 
-  if (settings::run_mode == RunMode::EIGENVALUE) {
-    // Add fission source if in eigenvalue mode
+  // Add fission source
 #pragma omp parallel for
-    for (int sr = 0; sr < n_source_regions_; sr++) {
-      int material = material_[sr];
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    int material = material_[sr];
 
-      for (int e_out = 0; e_out < negroups_; e_out++) {
-        float sigma_t = data::mg.macro_xs_[material].get_xs(
-          MgxsType::TOTAL, e_out, nullptr, nullptr, nullptr, t, a);
-        float fission_source = 0.0f;
+    for (int e_out = 0; e_out < negroups_; e_out++) {
+      float sigma_t = data::mg.macro_xs_[material].get_xs(
+        MgxsType::TOTAL, e_out, nullptr, nullptr, nullptr, t, a);
+      float fission_source = 0.0f;
 
-        for (int e_in = 0; e_in < negroups_; e_in++) {
-          float scalar_flux = scalar_flux_old_[sr * negroups_ + e_in];
-          float nu_sigma_f = data::mg.macro_xs_[material].get_xs(
-            MgxsType::NU_FISSION, e_in, nullptr, nullptr, nullptr, t, a);
-          float chi = data::mg.macro_xs_[material].get_xs(
-            MgxsType::CHI_PROMPT, e_in, &e_out, nullptr, nullptr, t, a);
-          fission_source += nu_sigma_f * scalar_flux * chi;
-        }
-        source_[sr * negroups_ + e_out] +=
-          fission_source * inverse_k_eff / sigma_t;
+      for (int e_in = 0; e_in < negroups_; e_in++) {
+        float scalar_flux = scalar_flux_old_[sr * negroups_ + e_in];
+        float nu_sigma_f = data::mg.macro_xs_[material].get_xs(
+          MgxsType::NU_FISSION, e_in, nullptr, nullptr, nullptr, t, a);
+        float chi = data::mg.macro_xs_[material].get_xs(
+          MgxsType::CHI_PROMPT, e_in, &e_out, nullptr, nullptr, t, a);
+        fission_source += nu_sigma_f * scalar_flux * chi;
       }
+      source_[sr * negroups_ + e_out] +=
+        fission_source * inverse_k_eff / sigma_t;
     }
-  } else {
-// Add external source if in fixed source mode
+  }
+
+  // Add external source if in fixed source mode
+  if (settings::run_mode == RunMode::FIXED_SOURCE) {
 #pragma omp parallel for
     for (int se = 0; se < n_source_elements_; se++) {
       source_[se] += external_source_[se];
diff --git a/tests/regression_tests/random_ray_fixed_source_subcritical/__init__.py b/tests/regression_tests/random_ray_fixed_source_subcritical/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat
new file mode 100644
index 000000000..0bf4cfdc9
--- /dev/null
+++ b/tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat
@@ -0,0 +1,140 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="UO2 fuel">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="UO2"/>
+    </material>
+    <material id="2" name="Water">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="LWTR"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="fuel inner a" region="-2" universe="1"/>
+    <cell id="2" material="1" name="fuel inner b" region="2 -3" universe="1"/>
+    <cell id="3" material="1" name="fuel inner c" region="3 -1" universe="1"/>
+    <cell id="4" material="2" name="moderator inner a" region="1 -4" universe="1"/>
+    <cell id="5" material="2" name="moderator outer b" region="4 -5" universe="1"/>
+    <cell id="6" material="2" name="moderator outer c" region="5" universe="1"/>
+    <cell fill="1" id="7" name="azimuthal_cell_0" region="6 -7" universe="2"/>
+    <cell fill="1" id="8" name="azimuthal_cell_1" region="7 -8" universe="2"/>
+    <cell fill="1" id="9" name="azimuthal_cell_2" region="8 -9" universe="2"/>
+    <cell fill="1" id="10" name="azimuthal_cell_3" region="9 -10" universe="2"/>
+    <cell fill="1" id="11" name="azimuthal_cell_4" region="10 -11" universe="2"/>
+    <cell fill="1" id="12" name="azimuthal_cell_5" region="11 -12" universe="2"/>
+    <cell fill="1" id="13" name="azimuthal_cell_6" region="12 -13" universe="2"/>
+    <cell fill="1" id="14" name="azimuthal_cell_7" region="13 -6" universe="2"/>
+    <cell id="18" material="2" name="moderator infinite" universe="8"/>
+    <cell id="19" material="2" name="moderator infinite 2" universe="9"/>
+    <cell fill="10" id="20" universe="12"/>
+    <cell fill="11" id="21" universe="13"/>
+    <cell fill="14" id="22" name="assembly" region="18 -19 20 -21" universe="15"/>
+    <lattice id="10">
+      <pitch>0.126 0.126</pitch>
+      <dimension>10 10</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <universes>
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 </universes>
+    </lattice>
+    <lattice id="11">
+      <pitch>0.126 0.126</pitch>
+      <dimension>10 10</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <universes>
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 9 </universes>
+    </lattice>
+    <lattice id="14">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+2 12 
+12 13 </universes>
+    </lattice>
+    <surface coeffs="0.0 0.0 0.54" id="1" name="Fuel OR" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.33" id="2" name="inner ring a" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.45" id="3" name="inner ring b" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.6" id="4" name="outer ring a" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.69" id="5" name="outer ring b" type="z-cylinder"/>
+    <surface coeffs="-0.0 1.0 0 0" id="6" type="plane"/>
+    <surface coeffs="-0.7071067811865475 0.7071067811865476 0 0" id="7" type="plane"/>
+    <surface coeffs="-1.0 6.123233995736766e-17 0 0" id="8" type="plane"/>
+    <surface coeffs="-0.7071067811865476 -0.7071067811865475 0 0" id="9" type="plane"/>
+    <surface coeffs="-1.2246467991473532e-16 -1.0 0 0" id="10" type="plane"/>
+    <surface coeffs="0.7071067811865475 -0.7071067811865477 0 0" id="11" type="plane"/>
+    <surface coeffs="1.0 -1.8369701987210297e-16 0 0" id="12" type="plane"/>
+    <surface coeffs="0.7071067811865477 0.7071067811865474 0 0" id="13" type="plane"/>
+    <surface boundary="reflective" coeffs="-1.26" id="18" name="minimum x" type="x-plane"/>
+    <surface boundary="reflective" coeffs="1.26" id="19" name="maximum x" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-1.26" id="20" name="minimum y" type="y-plane"/>
+    <surface boundary="reflective" coeffs="1.26" id="21" name="maximum y" type="y-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>30</particles>
+    <batches>125</batches>
+    <inactive>100</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="box">
+        <parameters>1.134 -1.26 -1.0 1.26 -1.134 1.0</parameters>
+      </space>
+      <energy type="discrete">
+        <parameters>2e-05 0.0735 20.0 200.0 2000.0 750000.0 2000000.0 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285</parameters>
+      </energy>
+      <constraints>
+        <domain_type>universe</domain_type>
+        <domain_ids>9</domain_ids>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>40.0</distance_active>
+      <distance_inactive>40.0</distance_inactive>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
+        </space>
+      </source>
+      <volume_normalized_flux_tallies>False</volume_normalized_flux_tallies>
+      <source_shape>flat</source_shape>
+    </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_fixed_source_subcritical/flat/results_true.dat b/tests/regression_tests/random_ray_fixed_source_subcritical/flat/results_true.dat
new file mode 100644
index 000000000..e9d2a733d
--- /dev/null
+++ b/tests/regression_tests/random_ray_fixed_source_subcritical/flat/results_true.dat
@@ -0,0 +1,169 @@
+tally 1:
+1.582116E+02
+1.005480E+03
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.897563E+01
+1.396228E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.038901E+01
+1.667203E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.465551E+01
+2.438101E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.260230E+01
+1.110141E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+9.438906E+01
+3.579643E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.698287E+01
+1.314208E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.002152E+02
+1.608182E+03
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+7.374897E+01
+2.181176E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.636086E+01
+2.787512E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+3.000395E+01
+3.611990E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.824715E+01
+1.360847E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+9.621963E+01
+3.714397E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.954480E+01
+1.436984E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+1.067903E+02
+4.585529E+02
+1.257560E+01
+6.364831E+00
+3.060650E+01
+3.770132E+01
+4.983480E+01
+9.972910E+01
+2.349541E+00
+2.217300E-01
+5.718312E+00
+1.313392E+00
+1.965462E+01
+1.548888E+01
+2.011494E-01
+1.622334E-03
+4.895573E-01
+9.609705E-03
+2.318671E+01
+2.155841E+01
+2.445032E-01
+2.397352E-03
+5.950720E-01
+1.420043E-02
+5.322867E+01
+1.136138E+02
+1.968847E-01
+1.554358E-03
+4.791838E-01
+9.207284E-03
+1.191175E+02
+5.692803E+02
+5.806982E-02
+1.354810E-04
+1.436897E-01
+8.295235E-04
+8.022963E+01
+2.589171E+02
+3.465139E-01
+4.882327E-03
+9.638109E-01
+3.777193E-02
+1.555305E+02
+9.711601E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.836421E+01
+1.367538E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.042089E+01
+1.673401E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.477708E+01
+2.463124E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.317851E+01
+1.134205E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+9.634473E+01
+3.724422E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.879813E+01
+1.396454E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
diff --git a/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat
new file mode 100644
index 000000000..f7572a072
--- /dev/null
+++ b/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat
@@ -0,0 +1,140 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="UO2 fuel">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="UO2"/>
+    </material>
+    <material id="2" name="Water">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="LWTR"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="fuel inner a" region="-2" universe="1"/>
+    <cell id="2" material="1" name="fuel inner b" region="2 -3" universe="1"/>
+    <cell id="3" material="1" name="fuel inner c" region="3 -1" universe="1"/>
+    <cell id="4" material="2" name="moderator inner a" region="1 -4" universe="1"/>
+    <cell id="5" material="2" name="moderator outer b" region="4 -5" universe="1"/>
+    <cell id="6" material="2" name="moderator outer c" region="5" universe="1"/>
+    <cell fill="1" id="7" name="azimuthal_cell_0" region="6 -7" universe="2"/>
+    <cell fill="1" id="8" name="azimuthal_cell_1" region="7 -8" universe="2"/>
+    <cell fill="1" id="9" name="azimuthal_cell_2" region="8 -9" universe="2"/>
+    <cell fill="1" id="10" name="azimuthal_cell_3" region="9 -10" universe="2"/>
+    <cell fill="1" id="11" name="azimuthal_cell_4" region="10 -11" universe="2"/>
+    <cell fill="1" id="12" name="azimuthal_cell_5" region="11 -12" universe="2"/>
+    <cell fill="1" id="13" name="azimuthal_cell_6" region="12 -13" universe="2"/>
+    <cell fill="1" id="14" name="azimuthal_cell_7" region="13 -6" universe="2"/>
+    <cell id="18" material="2" name="moderator infinite" universe="8"/>
+    <cell id="19" material="2" name="moderator infinite 2" universe="9"/>
+    <cell fill="10" id="20" universe="12"/>
+    <cell fill="11" id="21" universe="13"/>
+    <cell fill="14" id="22" name="assembly" region="18 -19 20 -21" universe="15"/>
+    <lattice id="10">
+      <pitch>0.126 0.126</pitch>
+      <dimension>10 10</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <universes>
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 </universes>
+    </lattice>
+    <lattice id="11">
+      <pitch>0.126 0.126</pitch>
+      <dimension>10 10</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <universes>
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 8 
+8 8 8 8 8 8 8 8 8 9 </universes>
+    </lattice>
+    <lattice id="14">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+2 12 
+12 13 </universes>
+    </lattice>
+    <surface coeffs="0.0 0.0 0.54" id="1" name="Fuel OR" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.33" id="2" name="inner ring a" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.45" id="3" name="inner ring b" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.6" id="4" name="outer ring a" type="z-cylinder"/>
+    <surface coeffs="0.0 0.0 0.69" id="5" name="outer ring b" type="z-cylinder"/>
+    <surface coeffs="-0.0 1.0 0 0" id="6" type="plane"/>
+    <surface coeffs="-0.7071067811865475 0.7071067811865476 0 0" id="7" type="plane"/>
+    <surface coeffs="-1.0 6.123233995736766e-17 0 0" id="8" type="plane"/>
+    <surface coeffs="-0.7071067811865476 -0.7071067811865475 0 0" id="9" type="plane"/>
+    <surface coeffs="-1.2246467991473532e-16 -1.0 0 0" id="10" type="plane"/>
+    <surface coeffs="0.7071067811865475 -0.7071067811865477 0 0" id="11" type="plane"/>
+    <surface coeffs="1.0 -1.8369701987210297e-16 0 0" id="12" type="plane"/>
+    <surface coeffs="0.7071067811865477 0.7071067811865474 0 0" id="13" type="plane"/>
+    <surface boundary="reflective" coeffs="-1.26" id="18" name="minimum x" type="x-plane"/>
+    <surface boundary="reflective" coeffs="1.26" id="19" name="maximum x" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-1.26" id="20" name="minimum y" type="y-plane"/>
+    <surface boundary="reflective" coeffs="1.26" id="21" name="maximum y" type="y-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>30</particles>
+    <batches>125</batches>
+    <inactive>100</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="box">
+        <parameters>1.134 -1.26 -1.0 1.26 -1.134 1.0</parameters>
+      </space>
+      <energy type="discrete">
+        <parameters>2e-05 0.0735 20.0 200.0 2000.0 750000.0 2000000.0 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285</parameters>
+      </energy>
+      <constraints>
+        <domain_type>universe</domain_type>
+        <domain_ids>9</domain_ids>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>40.0</distance_active>
+      <distance_inactive>40.0</distance_inactive>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
+        </space>
+      </source>
+      <volume_normalized_flux_tallies>False</volume_normalized_flux_tallies>
+      <source_shape>linear_xy</source_shape>
+    </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_fixed_source_subcritical/linear_xy/results_true.dat b/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/results_true.dat
new file mode 100644
index 000000000..90be5c569
--- /dev/null
+++ b/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/results_true.dat
@@ -0,0 +1,169 @@
+tally 1:
+1.573422E+02
+9.945961E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.868035E+01
+1.382578E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.028997E+01
+1.651549E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.453137E+01
+2.414255E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.232506E+01
+1.098683E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+9.385720E+01
+3.539613E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.664229E+01
+1.298333E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.001324E+02
+1.606858E+03
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+7.358736E+01
+2.171786E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.625416E+01
+2.765491E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.987675E+01
+3.582038E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.797311E+01
+1.348247E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+9.573509E+01
+3.677339E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.926615E+01
+1.423707E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+1.060513E+02
+4.521166E+02
+1.247077E+01
+6.257701E+00
+3.035137E+01
+3.706675E+01
+4.962504E+01
+9.887387E+01
+2.338861E+00
+2.196709E-01
+5.692320E+00
+1.301195E+00
+1.957564E+01
+1.536441E+01
+2.003053E-01
+1.608650E-03
+4.875030E-01
+9.528647E-03
+2.309275E+01
+2.138375E+01
+2.434816E-01
+2.377213E-03
+5.925856E-01
+1.408114E-02
+5.300741E+01
+1.126617E+02
+1.960322E-01
+1.540737E-03
+4.771089E-01
+9.126597E-03
+1.184809E+02
+5.630856E+02
+5.774249E-02
+1.339082E-04
+1.428798E-01
+8.198938E-04
+7.971031E+01
+2.554943E+02
+3.441954E-01
+4.814135E-03
+9.573621E-01
+3.724437E-02
+1.549617E+02
+9.639834E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.808547E+01
+1.354629E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.029953E+01
+1.654123E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.462331E+01
+2.433425E+01
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.283851E+01
+1.120036E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+9.572799E+01
+3.677288E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+5.842566E+01
+1.378719E+02
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
diff --git a/tests/regression_tests/random_ray_fixed_source_subcritical/test.py b/tests/regression_tests/random_ray_fixed_source_subcritical/test.py
new file mode 100644
index 000000000..6c2c569c6
--- /dev/null
+++ b/tests/regression_tests/random_ray_fixed_source_subcritical/test.py
@@ -0,0 +1,133 @@
+import os
+
+import openmc
+from openmc.examples import random_ray_lattice
+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("shape", ["flat", "linear_xy"])
+def test_random_ray_source(shape):
+    with change_directory(shape):
+        openmc.reset_auto_ids()
+
+        # The general strategy is to reuse the random_ray_lattice model,
+        # but redfine some of the geometry to make it a good
+        # subcritical multiplication problem. We then also add in
+        # a fixed source term.
+
+        model = random_ray_lattice()
+
+        # Begin by updating the random ray settings for fixed source
+        settings = model.settings
+        settings.random_ray['source_shape'] = shape
+        settings.run_mode = 'fixed source'
+        settings.particles = 30
+        settings.random_ray['distance_active'] = 40.0
+        settings.random_ray['distance_inactive'] = 40.0
+        settings.random_ray['volume_normalized_flux_tallies'] = False
+
+        # This problem needs about 2k iterations to converge,
+        # but for regression testing we only need a few hundred
+        # to ensure things are working as expected. With
+        # only 100 inactive batches, tallies will still be off
+        # by 3x or more. For validation against MGMC, be sure
+        # to increase the batch counts.
+        settings.batches = 125
+        settings.inactive = 100
+
+        ########################################
+        # Define the alternative geometry
+
+        pitch = 1.26
+
+        for material in model.materials:
+            if material.name == 'Water':
+                water = material
+
+        # The new geometry replaces two of the fuel pins with
+        # moderator, reducing k-eff to around 0.84. We also
+        # add a special universe in the corner of one of the moderator
+        # regions to use as a domain constraint for the source
+        moderator_infinite = openmc.Cell(fill=water, name='moderator infinite')
+        mu = openmc.Universe(cells=[moderator_infinite])
+
+        moderator_infinite2 = openmc.Cell(fill=water, name='moderator infinite 2')
+        mu2 = openmc.Universe(cells=[moderator_infinite2])
+
+        n_sub = 10
+
+        lattice = openmc.RectLattice()
+        lattice.lower_left = [-pitch/2.0, -pitch/2.0]
+        lattice.pitch = [pitch/n_sub, pitch/n_sub]
+        lattice.universes = [[mu] * n_sub for _ in range(n_sub)]
+
+        lattice2 = openmc.RectLattice()
+        lattice2.lower_left = [-pitch/2.0, -pitch/2.0]
+        lattice2.pitch = [pitch/n_sub, pitch/n_sub]
+        lattice2.universes = [[mu] * n_sub for _ in range(n_sub)]
+        lattice2.universes[n_sub-1][n_sub-1] = mu2
+
+        mod_lattice_cell = openmc.Cell(fill=lattice)
+        mod_lattice_uni = openmc.Universe(cells=[mod_lattice_cell])
+
+        mod_lattice_cell2 = openmc.Cell(fill=lattice2)
+        mod_lattice_uni2 = openmc.Universe(cells=[mod_lattice_cell2])
+
+        lattice2x2 = openmc.RectLattice()
+        lattice2x2.lower_left = [-pitch, -pitch]
+        lattice2x2.pitch = [pitch, pitch]
+
+        universes = model.geometry.get_all_universes()
+        for universe in universes.values():
+            if universe.name == 'pincell':
+                pincell = universe
+
+        lattice2x2.universes = [
+            [pincell, mod_lattice_uni],
+            [mod_lattice_uni, mod_lattice_uni2]
+        ]
+
+        box = openmc.model.RectangularPrism(
+            pitch*2, pitch*2, boundary_type='reflective')
+
+        assembly = openmc.Cell(fill=lattice2x2, region=-box, name='assembly')
+
+        root = openmc.Universe(name='root universe', cells=[assembly])
+        model.geometry = openmc.Geometry(root)
+
+        ########################################
+        # Define the fixed source term
+
+        s = 1.0 / 7.0
+        strengths = [s, s, s, s, s, s, s]
+        midpoints = [2.0e-5, 0.0735, 20.0, 2.0e2, 2.0e3, 0.75e6, 2.0e6]
+        energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths)
+
+        lower_left_src = [pitch - pitch/10.0, -pitch, -1.0]
+        upper_right_src = [pitch, -pitch + pitch/10.0, 1.0]
+        spatial_distribution = openmc.stats.Box(
+            lower_left_src, upper_right_src, only_fissionable=False)
+
+        settings.source = openmc.IndependentSource(
+            space=spatial_distribution,
+            energy=energy_distribution,
+            constraints={'domains': [mu2]},
+            strength=1.0
+        )
+
+        ########################################
+        # Run test
+
+        harness = MGXSTestHarness('statepoint.125.h5', model)
+        harness.main()