Fixed Source Random Ray (#2988)

Co-authored-by: Gavin Ridley <gavin.keith.ridley@gmail.com> Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
2026-07-28 06:05:58 -04:00 · 2024-06-17 11:02:20 -05:00 · 2024-06-17 11:02:20 -05:00 · 5222b343a4
commit 5222b343a4
parent b1b8a4c328
21 changed files with 2064 additions and 369 deletions
--- a/docs/source/methods/random_ray.rst
+++ b/docs/source/methods/random_ray.rst
@ -94,8 +94,8 @@ Method of Characteristics
 The Boltzmann neutron transport equation is a partial differential equation
 (PDE) that describes the angular flux within a system. It is a balance equation,
 with the streaming and absorption terms typically appearing on the left hand
-side, which are balanced by the scattering source and fission source terms on
-the right hand side.
+side, which are balanced by the scattering source, fission, and fixed source
+terms on the right hand side.

 .. math::
    :label: transport
@ -522,8 +522,8 @@ make their traversals, and summing these contributions up as in Equation
 improve the estimate of the source and scalar flux over many iterations, given
 that our initial starting source will just be a guess?

-The source :math:`Q^{n}` for iteration :math:`n` can be inferred
-from the scalar flux from the previous iteration :math:`n-1` as:
+In an eigenvalue simulation, the source :math:`Q^{n}` for iteration :math:`n`
+can be inferred from the scalar flux from the previous iteration :math:`n-1` as:

 .. math::
    :label: source_update
@ -535,7 +535,7 @@ where :math:`Q^{n}(i, g)` is the total source (fission + scattering) in region
 :math:`g` must be computed by summing over the contributions from all groups
 :math:`g' \in G`.

-In a similar manner, the eigenvalue for iteration :math:`n` can be computed as:
+The eigenvalue for iteration :math:`n` can be computed as:

 .. math::
    :label: eigenvalue_update
@ -576,6 +576,18 @@ and a similar substitution can be made to update Equation
 estimate is used, such that the total fission source from the previous iteration
 (:math:`n-1`) is also recomputed each iteration.

+In a fixed source simulation, the fission source is replaced by a user specified
+fixed source term :math:`Q_\text{fixed}(i,E)`, which is defined for each FSR and
+energy group. This additional source term is applied at this stage for
+generating the next iteration's source estimate as:
+
+.. math::
+    :label: fixed_source_update
+
+    Q^{n}(i, g) = Q_\text{fixed}(i,g) + \sum\limits^{G}_{g'} \Sigma_{s}(i,g,g') \phi^{n-1}(g')
+
+and no eigenvalue is computed.
+
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 Ray Starting Conditions and Inactive Length
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -742,6 +754,32 @@ behavior if a single simulation cell is able to score to multiple filter mesh
 cells. In the future, the capability to fully support mesh tallies may be added
 to OpenMC, but for now this restriction needs to be respected.

+.. _usersguide_fixed_source_methods:
+
+------------
+Fixed Source
+------------
+
+The random ray solver in OpenMC can be used for both eigenvalue and fixed source
+problems. There are a few key differences between fixed source transport with
+random ray and Monte Carlo, however.
+
+- **Source definition:** In Monte Carlo, it is relatively easy to define various
+  source distributions, including point sources, surface sources, volume
+  sources, and even custom user sources -- all with varying angular and spatial
+  statistical distributions. In random ray, the natural way to include a fixed
+  source term is by adding a fixed (flat) contribution to specific flat source
+  regions. Thus, in the OpenMC implementation of random ray, particle sources
+  are restricted to being volumetric and isotropic, although different energy
+  spectrums are supported. Fixed sources can be applied to specific materials,
+  cells, or universes.
+
+- **Inactive batches:** In Monte Carlo, use of a fixed source implies that all
+  batches are active batches, as there is no longer a need to develop a fission
+  source distribution. However, in random ray mode, there is still a need to
+  develop the scattering source by way of inactive batches before beginning
+  active batches.
+
 ---------------------------
 Fundamental Sources of Bias
 ---------------------------
--- a/docs/source/usersguide/random_ray.rst
+++ b/docs/source/usersguide/random_ray.rst
@ -40,13 +40,15 @@ Carlo, **inactive batches are required for both eigenvalue and fixed source
 solves in random ray mode** due to this additional need to converge the
 scattering source.

+.. warning::
+    Unlike Monte Carlo, the random ray solver still requires usage of inactive
+    batches when in fixed source mode so as to develop the scattering source.
+
 The additional burden of converging the scattering source generally results in a
 higher requirement for the number of inactive batches---often by an order of
 magnitude or more. For instance, it may be reasonable to only use 50 inactive
 batches for a light water reactor simulation with Monte Carlo, but random ray
-might require 500 or more inactive batches. Similar to Monte Carlo,
-:ref:`Shannon entropy <usersguide_entropy>` can be used to gauge whether the
-combined scattering and fission source has fully developed.
+might require 500 or more inactive batches.

 Similar to Monte Carlo, active batches are used in the random ray solver mode to
 accumulate and converge statistics on unknown quantities (i.e., the random ray
@ -248,6 +250,8 @@ a larger value until the "low ray density" messages go away.
    ray lengths are sufficiently long to allow for transport to occur between
    source and target regions of interest.

+.. _usersguide_ray_source:
+
 ----------
 Ray Source
 ----------
@ -261,7 +265,7 @@ that the source must not be limited to only fissionable regions. Additionally,
 the source box must cover the entire simulation domain. In the case of a
 simulation domain that is not box shaped, a box source should still be used to
 bound the domain but with the source limited to rejection sampling the actual
-simulation universe (which can be specified via the ``domains`` field of the
+simulation universe (which can be specified via the ``domains`` constraint of the
 :class:`openmc.IndependentSource` Python class). Similar to Monte Carlo sources,
 for two-dimensional problems (e.g., a 2D pincell) it is desirable to make the
 source bounded near the origin of the infinite dimension. An example of an
@ -411,11 +415,78 @@ in the `OpenMC Jupyter notebook collection
    separate materials can be defined each with a separate multigroup dataset
    corresponding to a given temperature.

+---------------------------------
+Fixed Source and Eigenvalue Modes
+---------------------------------
+
+Both fixed source and eigenvalue modes are supported with the random ray solver
+in OpenMC. Modes can be selected as described in the :ref:`run modes section
+<usersguide_run_modes>`. In both modes, a ray source must be provided to let
+OpenMC know where to sample ray starting locations from, as discussed in the
+:ref:`ray source section <usersguide_ray_source>`. In fixed source mode, at
+least one regular source must be provided as well that represents the physical
+particle fixed source. As discussed in the :ref:`fixed source methodology
+section <usersguide_fixed_source_methods>`, the types of fixed sources supported
+in the random ray solver mode are limited compared to what is possible with the
+Monte Carlo solver.
+
+Currently, all of the following conditions must be met for the particle source
+to be valid in random ray mode:
+
+- One or more domain ids must be specified that indicate which cells, universes,
+  or materials the source applies to. This implicitly limits the source type to
+  being volumetric. This is specified via the ``domains`` constraint placed on the
+  :class:`openmc.IndependentSource` Python class.
+- The source must be isotropic (default for a source)
+- The source must use a discrete (i.e., multigroup) energy distribution. The
+  discrete energy distribution is input by defining a
+  :class:`openmc.stats.Discrete` Python class, and passed as the ``energy``
+  field of the :class:`openmc.IndependentSource` Python class.
+
+Any other spatial distribution information contained in a particle source will
+be ignored. Only the specified cell, material, or universe domains will be used
+to define the spatial location of the source, as the source will be applied
+during a pre-processing stage of OpenMC to all source regions that are contained
+within the specified domains for the source.
+
+When defining a :class:`openmc.stats.Discrete` object, note that the ``x`` field
+will correspond to the discrete energy points, and the ``p`` field will
+correspond to the discrete probabilities. It is recommended to select energy
+points that fall within energy groups rather than on boundaries between the
+groups. That is, if the problem contains two energy groups (with bin edges of
+1.0e-5, 1.0e-1, 1.0e7), then a good selection for the ``x`` field might be
+points of 1.0e-2 and 1.0e1.
+
+::
+
+    # Define geometry, etc.
+    ...
+    source_cell = openmc.Cell(fill=source_mat, name='cell where fixed source will be')
+    ...
+    # Define physical neutron fixed source
+    energy_points = [1.0e-2, 1.0e1]
+    strengths = [0.25, 0.75]
+    energy_distribution = openmc.stats.Discrete(x=energy_points, p=strengths)
+    neutron_source = openmc.IndependentSource(
+        energy=energy_distribution,
+        constraints={'domains': [source_cell]}
+    )
+
+    # Add fixed source and ray sampling source to settings file
+    settings.source = [neutron_source]
+
 ---------------------------------------
 Putting it All Together: Example Inputs
 ---------------------------------------

-An example of a settings definition for random ray is given below::
+~~~~~~~~~~~~~~~~~~
+Eigenvalue Example
+~~~~~~~~~~~~~~~~~~
+
+An example of a settings definition for an eigenvalue random ray simulation is
+given below:
+
+::

    # Geometry and MGXS material definition of 2x2 lattice (not shown)
    pitch = 1.26
@ -478,3 +549,84 @@ Monte Carlo run (see the :ref:`geometry <usersguide_geometry>` and

 There is also a complete example of a pincell available in the
 ``openmc/examples/pincell_random_ray`` folder.
+
+~~~~~~~~~~~~~~~~~~~~
+Fixed Source Example
+~~~~~~~~~~~~~~~~~~~~
+
+An example of a settings definition for a fixed source random ray simulation is
+given below:
+
+::
+
+    # Geometry and MGXS material definition of 2x2 lattice (not shown)
+    pitch = 1.26
+    source_cell = openmc.Cell(fill=source_mat, name='cell where fixed source will be')
+    ebins = [1e-5, 1e-1, 20.0e6]
+    ...
+
+    # Instantiate a settings object for a random ray solve
+    settings = openmc.Settings()
+    settings.energy_mode = "multi-group"
+    settings.batches = 1200
+    settings.inactive = 600
+    settings.particles = 2000
+    settings.run_mode = 'fixed source'
+    settings.random_ray['distance_inactive'] = 40.0
+    settings.random_ray['distance_active'] = 400.0
+
+    # Create an initial uniform spatial source distribution for sampling rays
+    lower_left  = (-pitch, -pitch, -pitch)
+    upper_right = ( pitch,  pitch,  pitch)
+    uniform_dist = openmc.stats.Box(lower_left, upper_right)
+    settings.random_ray['ray_source'] = openmc.IndependentSource(space=uniform_dist)
+
+    # Define physical neutron fixed source
+    energy_points = [1.0e-2, 1.0e1]
+    strengths = [0.25, 0.75]
+    energy_distribution = openmc.stats.Discrete(x=energy_points, p=strengths)
+    neutron_source = openmc.IndependentSource(
+        energy=energy_distribution,
+        constraints={'domains': [source_cell]}
+    )
+
+    # Add fixed source and ray sampling source to settings file
+    settings.source = [neutron_source]
+
+    settings.export_to_xml()
+
+    # Define tallies
+
+    # Create a mesh filter
+    mesh = openmc.RegularMesh()
+    mesh.dimension = (2, 2)
+    mesh.lower_left = (-pitch/2, -pitch/2)
+    mesh.upper_right = (pitch/2, pitch/2)
+    mesh_filter = openmc.MeshFilter(mesh)
+
+    # Create a multigroup energy filter
+    energy_filter = openmc.EnergyFilter(ebins)
+
+    # Create tally using our two filters and add scores
+    tally = openmc.Tally()
+    tally.filters = [mesh_filter, energy_filter]
+    tally.scores = ['flux']
+
+    # Instantiate a Tallies collection and export to XML
+    tallies = openmc.Tallies([tally])
+    tallies.export_to_xml()
+
+    # Create voxel plot
+    plot = openmc.Plot()
+    plot.origin = [0, 0, 0]
+    plot.width = [2*pitch, 2*pitch, 1]
+    plot.pixels = [1000, 1000, 1]
+    plot.type = 'voxel'
+
+    # Instantiate a Plots collection and export to XML
+    plots = openmc.Plots([plot])
+    plots.export_to_xml()
+
+All other inputs (e.g., geometry, material) will be unchanged from a typical
+Monte Carlo run (see the :ref:`geometry <usersguide_geometry>` and
+:ref:`multigroup materials <create_mgxs>` user guides for more information).
--- a/include/openmc/mgxs_interface.h
+++ b/include/openmc/mgxs_interface.h
@ -46,6 +46,9 @@ public:
  // Get the kT values which are used in the OpenMC model
  vector<vector<double>> get_mat_kTs();

+  // Get the group index corresponding to a continuous energy
+  int get_group_index(double E);
+
  int num_energy_groups_;
  int num_delayed_groups_;
  vector<std::string> xs_names_;            // available names in HDF5 file
--- a/include/openmc/random_ray/flat_source_domain.h
+++ b/include/openmc/random_ray/flat_source_domain.h
@ -3,9 +3,83 @@

 #include "openmc/openmp_interface.h"
 #include "openmc/position.h"
+#include "openmc/source.h"

 namespace openmc {

+//----------------------------------------------------------------------------
+// Helper Functions
+
+// The hash_combine function is the standard hash combine function from boost
+// that is typically used for combining multiple hash values into a single hash
+// as is needed for larger objects being stored in a hash map. The function is
+// taken from:
+// https://www.boost.org/doc/libs/1_55_0/doc/html/hash/reference.html#boost.hash_combine
+// which carries the following license:
+//
+// Boost Software License - Version 1.0 - August 17th, 2003
+// Permission is hereby granted, free of charge, to any person or organization
+// obtaining a copy of the software and accompanying documentation covered by
+// this license (the "Software") to use, reproduce, display, distribute,
+// execute, and transmit the Software, and to prepare derivative works of the
+// Software, and to permit third-parties to whom the Software is furnished to
+// do so, all subject to the following:
+// The copyright notices in the Software and this entire statement, including
+// the above license grant, this restriction and the following disclaimer,
+// must be included in all copies of the Software, in whole or in part, and
+// all derivative works of the Software, unless such copies or derivative
+// works are solely in the form of machine-executable object code generated by
+// a source language processor.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
+// SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
+// FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
+// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+// DEALINGS IN THE SOFTWARE.
+inline void hash_combine(size_t& seed, const size_t v)
+{
+  seed ^= (v + 0x9e3779b9 + (seed << 6) + (seed >> 2));
+}
+
+//----------------------------------------------------------------------------
+// Helper Structs
+
+// A mapping object that is used to map between a specific random ray
+// source region and an OpenMC native tally bin that it should score to
+// every iteration.
+struct TallyTask {
+  int tally_idx;
+  int filter_idx;
+  int score_idx;
+  int score_type;
+  TallyTask(int tally_idx, int filter_idx, int score_idx, int score_type)
+    : tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx),
+      score_type(score_type)
+  {}
+  TallyTask() = default;
+
+  // Comparison and Hash operators are defined to allow usage of the
+  // TallyTask struct as a key in an unordered_set
+  bool operator==(const TallyTask& other) const
+  {
+    return tally_idx == other.tally_idx && filter_idx == other.filter_idx &&
+           score_idx == other.score_idx && score_type == other.score_type;
+  }
+
+  struct HashFunctor {
+    size_t operator()(const TallyTask& task) const
+    {
+      size_t seed = 0;
+      hash_combine(seed, task.tally_idx);
+      hash_combine(seed, task.filter_idx);
+      hash_combine(seed, task.score_idx);
+      hash_combine(seed, task.score_type);
+      return seed;
+    }
+  };
+};
+
 /*
 * The FlatSourceDomain class encompasses data and methods for storing
 * scalar flux and source region for all flat source regions in a
@ -14,23 +88,6 @@ namespace openmc {

 class FlatSourceDomain {
 public:
-  //----------------------------------------------------------------------------
-  // Helper Structs
-
-  // A mapping object that is used to map between a specific random ray
-  // source region and an OpenMC native tally bin that it should score to
-  // every iteration.
-  struct TallyTask {
-    int tally_idx;
-    int filter_idx;
-    int score_idx;
-    int score_type;
-    TallyTask(int tally_idx, int filter_idx, int score_idx, int score_type)
-      : tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx),
-        score_type(score_type)
-    {}
-  };
-
  //----------------------------------------------------------------------------
  // Constructors
  FlatSourceDomain();
@ -44,10 +101,13 @@ public:
  int64_t add_source_to_scalar_flux();
  void batch_reset();
  void convert_source_regions_to_tallies();
-  void random_ray_tally() const;
+  void reset_tally_volumes();
+  void random_ray_tally();
  void accumulate_iteration_flux();
  void output_to_vtk() const;
  void all_reduce_replicated_source_regions();
+  void convert_external_sources();
+  void count_external_source_regions();

  //----------------------------------------------------------------------------
  // Public Data members
@ -55,6 +115,8 @@ public:
  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
+  int64_t n_external_source_regions_ {0}; // Total number of source regions with
+                                          // non-zero external source terms

  // 1D array representing source region starting offset for each OpenMC Cell
  // in model::cells
@ -72,18 +134,39 @@ public:
  vector<float> scalar_flux_old_;
  vector<float> scalar_flux_new_;
  vector<float> source_;
+  vector<float> external_source_;
+
+private:
+  //----------------------------------------------------------------------------
+  // Methods
+  void apply_external_source_to_source_region(
+    Discrete* discrete, double strength_factor, int64_t source_region);
+  void apply_external_source_to_cell_instances(int32_t i_cell,
+    Discrete* discrete, double strength_factor, int target_material_id,
+    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);

  //----------------------------------------------------------------------------
  // Private data members
-private:
  int negroups_;                  // Number of energy groups in simulation
  int64_t n_source_elements_ {0}; // Total number of source regions in the model
                                  // times the number of energy groups

-  // 2D array representing values for all source regions x energy groups x tally
+  double
+    simulation_volume_; // Total physical volume of the simulation domain, as
+                        // defined by the 3D box of the random ray source
+
+  // 2D array representing values for all source elements x tally
  // tasks
  vector<vector<TallyTask>> tally_task_;

+  // 1D array representing values for all source regions, with each region
+  // containing a set of volume tally tasks. This more complicated data
+  // structure is convenient for ensuring that volumes are only tallied once per
+  // source region, regardless of how many energy groups are used for tallying.
+  vector<std::unordered_set<TallyTask, TallyTask::HashFunctor>> volume_task_;
+
  // 1D arrays representing values for all source regions
  vector<int> material_;
  vector<double> volume_t_;
@ -92,6 +175,14 @@ private:
  // groups
  vector<float> scalar_flux_final_;

+  // Volumes for each tally and bin/score combination. This intermediate data
+  // structure is used when tallying quantities that must be normalized by
+  // volume (i.e., flux). The vector is index by tally index, while the inner 2D
+  // xtensor is indexed by bin index and score index in a similar manner to the
+  // results tensor in the Tally class, though without the third dimension, as
+  // SUM and SUM_SQ do not need to be tracked.
+  vector<xt::xtensor<double, 2>> tally_volumes_;
+
 }; // class FlatSourceDomain

 //============================================================================
--- a/include/openmc/random_ray/random_ray_simulation.h
+++ b/include/openmc/random_ray/random_ray_simulation.h
@ -24,7 +24,8 @@ public:
  void instability_check(
    int64_t n_hits, double k_eff, double& avg_miss_rate) const;
  void print_results_random_ray(uint64_t total_geometric_intersections,
-    double avg_miss_rate, int negroups, int64_t n_source_regions) const;
+    double avg_miss_rate, int negroups, int64_t n_source_regions,
+    int64_t n_external_source_regions) const;

  //----------------------------------------------------------------------------
  // Data members
--- a/include/openmc/source.h
+++ b/include/openmc/source.h
@ -52,6 +52,8 @@ extern vector<unique_ptr<Source>> external_sources;

 class Source {
 public:
+  // Domain types
+  enum class DomainType { UNIVERSE, MATERIAL, CELL };
  // Constructors, destructors
  Source() = default;
  explicit Source(pugi::xml_node node);
@ -76,9 +78,6 @@ public:
  static unique_ptr<Source> create(pugi::xml_node node);

 protected:
-  // Domain types
-  enum class DomainType { UNIVERSE, MATERIAL, CELL };
-
  // Strategy used for rejecting sites when constraints are applied. KILL means
  // that sites are always accepted but if they don't satisfy constraints, they
  // are given weight 0. RESAMPLE means that a new source site will be sampled
@ -134,6 +133,10 @@ public:
  Distribution* energy() const { return energy_.get(); }
  Distribution* time() const { return time_.get(); }

+  // Make domain type and ids available
+  DomainType domain_type() const { return domain_type_; }
+  const std::unordered_set<int32_t>& domain_ids() const { return domain_ids_; }
+
 protected:
  // Indicates whether derived class already handles constraints
  bool constraints_applied() const override { return true; }
--- a/src/mgxs_interface.cpp
+++ b/src/mgxs_interface.cpp
@ -15,6 +15,7 @@
 #include "openmc/material.h"
 #include "openmc/math_functions.h"
 #include "openmc/nuclide.h"
+#include "openmc/search.h"
 #include "openmc/settings.h"

 namespace openmc {
@ -183,6 +184,15 @@ vector<vector<double>> MgxsInterface::get_mat_kTs()

 //==============================================================================

+int MgxsInterface::get_group_index(double E)
+{
+  int g =
+    lower_bound_index(rev_energy_bins_.begin(), rev_energy_bins_.end(), E);
+  return num_energy_groups_ - g - 1.;
+}
+
+//==============================================================================
+
 void MgxsInterface::read_header(const std::string& path_cross_sections)
 {
  // Save name of HDF5 file to be read to struct data
--- a/src/random_ray/flat_source_domain.cpp
+++ b/src/random_ray/flat_source_domain.cpp
@ -2,6 +2,7 @@

 #include "openmc/cell.h"
 #include "openmc/geometry.h"
+#include "openmc/material.h"
 #include "openmc/message_passing.h"
 #include "openmc/mgxs_interface.h"
 #include "openmc/output.h"
@ -48,10 +49,19 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)

  // Initialize element-wise arrays
  scalar_flux_new_.assign(n_source_elements_, 0.0);
-  scalar_flux_old_.assign(n_source_elements_, 1.0);
  scalar_flux_final_.assign(n_source_elements_, 0.0);
  source_.resize(n_source_elements_);
+  external_source_.assign(n_source_elements_, 0.0);
  tally_task_.resize(n_source_elements_);
+  volume_task_.resize(n_source_regions_);
+
+  if (settings::run_mode == RunMode::EIGENVALUE) {
+    // If in eigenvalue mode, set starting flux to guess of unity
+    scalar_flux_old_.assign(n_source_elements_, 1.0);
+  } else {
+    // If in fixed source mode, set starting flux to guess of zero
+    scalar_flux_old_.assign(n_source_elements_, 0.0);
+  }

  // Initialize material array
  int64_t source_region_id = 0;
@ -68,6 +78,25 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
  if (source_region_id != n_source_regions_) {
    fatal_error("Unexpected number of source regions");
  }
+
+  // Initialize tally volumes
+  tally_volumes_.resize(model::tallies.size());
+  for (int i = 0; i < model::tallies.size(); i++) {
+    //  Get the shape of the 3D result tensor
+    auto shape = model::tallies[i]->results().shape();
+
+    // Create a new 2D tensor with the same size as the first
+    // two dimensions of the 3D tensor
+    tally_volumes_[i] =
+      xt::xtensor<double, 2>::from_shape({shape[0], shape[1]});
+  }
+
+  // Compute simulation domain volume based on ray source
+  auto* is = dynamic_cast<IndependentSource*>(RandomRay::ray_source_.get());
+  SpatialDistribution* space_dist = is->space();
+  SpatialBox* sb = dynamic_cast<SpatialBox*>(space_dist);
+  Position dims = sb->upper_right() - sb->lower_left();
+  simulation_volume_ = dims.x * dims.y * dims.z;
 }

 void FlatSourceDomain::batch_reset()
@ -97,11 +126,11 @@ void FlatSourceDomain::update_neutron_source(double k_eff)

  // Temperature and angle indices, if using multiple temperature
  // data sets and/or anisotropic data sets.
-  // TODO: Currently assumes we are only using single temp/single
-  // angle data.
+  // TODO: Currently assumes we are only using single temp/single angle data.
  const int t = 0;
  const int a = 0;

+  // Add scattering source
 #pragma omp parallel for
  for (int sr = 0; sr < n_source_regions_; sr++) {
    int material = material_[sr];
@ -110,23 +139,47 @@ void FlatSourceDomain::update_neutron_source(double k_eff)
      float sigma_t = data::mg.macro_xs_[material].get_xs(
        MgxsType::TOTAL, e_out, nullptr, nullptr, nullptr, t, a);
      float scatter_source = 0.0f;
-      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 sigma_s = data::mg.macro_xs_[material].get_xs(
          MgxsType::NU_SCATTER, e_in, &e_out, nullptr, nullptr, t, a);
-        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);
        scatter_source += sigma_s * scalar_flux;
-        fission_source += nu_sigma_f * scalar_flux * chi;
      }

-      fission_source *= inverse_k_eff;
-      float new_isotropic_source = (scatter_source + fission_source) / sigma_t;
-      source_[sr * negroups_ + e_out] = new_isotropic_source;
+      source_[sr * negroups_ + e_out] = scatter_source / sigma_t;
+    }
+  }
+
+  if (settings::run_mode == RunMode::EIGENVALUE) {
+    // Add fission source if in eigenvalue mode
+#pragma omp parallel for
+    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_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
+#pragma omp parallel for
+    for (int se = 0; se < n_source_elements_; se++) {
+      source_[se] += external_source_[se];
    }
  }

@ -355,10 +408,14 @@ void FlatSourceDomain::convert_source_regions_to_tallies()
          for (auto score_index = 0; score_index < tally.scores_.size();
               score_index++) {
            auto score_bin = tally.scores_[score_index];
-            // If a valid tally, filter, and score cobination has been found,
+            // If a valid tally, filter, and score combination has been found,
            // then add it to the list of tally tasks for this source element.
-            tally_task_[source_element].emplace_back(
-              i_tally, filter_index, score_index, score_bin);
+            TallyTask task(i_tally, filter_index, score_index, score_bin);
+            tally_task_[source_element].push_back(task);
+
+            // Also add this task to the list of volume tasks for this source
+            // region.
+            volume_task_[sr].insert(task);
          }
        }
      }
@ -372,6 +429,16 @@ void FlatSourceDomain::convert_source_regions_to_tallies()
  mapped_all_tallies_ = all_source_regions_mapped;
 }

+// Set the volume accumulators to zero for all tallies
+void FlatSourceDomain::reset_tally_volumes()
+{
+#pragma omp parallel for
+  for (int i = 0; i < tally_volumes_.size(); i++) {
+    auto& tensor = tally_volumes_[i];
+    tensor.fill(0.0); // Set all elements of the tensor to 0.0
+  }
+}
+
 // Tallying in random ray is not done directly during transport, rather,
 // it is done only once after each power iteration. This is made possible
 // by way of a mapping data structure that relates spatial source regions
@ -381,10 +448,13 @@ void FlatSourceDomain::convert_source_regions_to_tallies()
 // tally function simply traverses the mapping data structure and executes
 // the scoring operations to OpenMC's native tally result arrays.

-void FlatSourceDomain::random_ray_tally() const
+void FlatSourceDomain::random_ray_tally()
 {
  openmc::simulation::time_tallies.start();

+  // Reset our tally volumes to zero
+  reset_tally_volumes();
+
  // Temperature and angle indices, if using multiple temperature
  // data sets and/or anisotropic data sets.
  // TODO: Currently assumes we are only using single temp/single
@ -397,32 +467,46 @@ void FlatSourceDomain::random_ray_tally() const
 // them.
 #pragma omp parallel for
  for (int sr = 0; sr < n_source_regions_; sr++) {
-    double volume = volume_[sr];
+    // The fsr.volume_ is the unitless fractional simulation averaged volume
+    // (i.e., it is the FSR's fraction of the overall simulation volume). The
+    // simulation_volume_ is the total 3D physical volume in cm^3 of the entire
+    // global simulation domain (as defined by the ray source box). Thus, the
+    // FSR's true 3D spatial volume in cm^3 is found by multiplying its fraction
+    // of the total volume by the total volume. Not important in eigenvalue
+    // solves, but useful in fixed source solves for returning the flux shape
+    // with a magnitude that makes sense relative to the fixed source strength.
+    double volume = volume_[sr] * simulation_volume_;
+
    double material = material_[sr];
    for (int g = 0; g < negroups_; g++) {
      int idx = sr * negroups_ + g;
-      double flux = scalar_flux_new_[idx] * volume;
+      double flux = scalar_flux_new_[idx];
+
+      // Determine numerical score value
      for (auto& task : tally_task_[idx]) {
        double score;
        switch (task.score_type) {

        case SCORE_FLUX:
-          score = flux;
+          score = flux * volume;
          break;

        case SCORE_TOTAL:
-          score = flux * data::mg.macro_xs_[material].get_xs(
-                           MgxsType::TOTAL, g, NULL, NULL, NULL, t, a);
+          score = flux * volume *
+                  data::mg.macro_xs_[material].get_xs(
+                    MgxsType::TOTAL, g, NULL, NULL, NULL, t, a);
          break;

        case SCORE_FISSION:
-          score = flux * data::mg.macro_xs_[material].get_xs(
-                           MgxsType::FISSION, g, NULL, NULL, NULL, t, a);
+          score = flux * volume *
+                  data::mg.macro_xs_[material].get_xs(
+                    MgxsType::FISSION, g, NULL, NULL, NULL, t, a);
          break;

        case SCORE_NU_FISSION:
-          score = flux * data::mg.macro_xs_[material].get_xs(
-                           MgxsType::NU_FISSION, g, NULL, NULL, NULL, t, a);
+          score = flux * volume *
+                  data::mg.macro_xs_[material].get_xs(
+                    MgxsType::NU_FISSION, g, NULL, NULL, NULL, t, a);
          break;

        case SCORE_EVENTS:
@ -435,13 +519,49 @@ void FlatSourceDomain::random_ray_tally() const
                      "random ray mode.");
          break;
        }
+
+        // Apply score to the appropriate tally bin
        Tally& tally {*model::tallies[task.tally_idx]};
 #pragma omp atomic
        tally.results_(task.filter_idx, task.score_idx, TallyResult::VALUE) +=
          score;
+      } // end tally task loop
+    }   // end energy group loop
+
+    // For flux tallies, the total volume of the spatial region is needed
+    // for normalizing the flux. We store this volume in a separate tensor.
+    // We only contribute to each volume tally bin once per FSR.
+    for (const auto& task : volume_task_[sr]) {
+      if (task.score_type == SCORE_FLUX) {
+#pragma omp atomic
+        tally_volumes_[task.tally_idx](task.filter_idx, task.score_idx) +=
+          volume;
+      }
+    }
+  } // end FSR loop
+
+  // Normalize any flux scores by the total volume of the FSRs scoring to that
+  // bin. To do this, we loop over all tallies, and then all filter bins,
+  // and then scores. For each score, we check the tally data structure to
+  // see what index that score corresponds to. If that score is a flux score,
+  // then we divide it by volume.
+  for (int i = 0; i < model::tallies.size(); i++) {
+    Tally& tally {*model::tallies[i]};
+#pragma omp parallel for
+    for (int bin = 0; bin < tally.n_filter_bins(); bin++) {
+      for (int score_idx = 0; score_idx < tally.n_scores(); score_idx++) {
+        auto score_type = tally.scores_[score_idx];
+        if (score_type == SCORE_FLUX) {
+          double vol = tally_volumes_[i](bin, score_idx);
+          if (vol > 0.0) {
+            tally.results_(bin, score_idx, TallyResult::VALUE) /= vol;
+          }
+        }
      }
    }
  }
+
+  openmc::simulation::time_tallies.stop();
 }

 void FlatSourceDomain::all_reduce_replicated_source_regions()
@ -683,4 +803,130 @@ void FlatSourceDomain::output_to_vtk() const
  }
 }

+void FlatSourceDomain::apply_external_source_to_source_region(
+  Discrete* discrete, double strength_factor, int64_t source_region)
+{
+  const auto& discrete_energies = discrete->x();
+  const auto& discrete_probs = discrete->prob();
+
+  for (int e = 0; e < discrete_energies.size(); e++) {
+    int g = data::mg.get_group_index(discrete_energies[e]);
+    external_source_[source_region * negroups_ + g] +=
+      discrete_probs[e] * strength_factor;
+  }
+}
+
+void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell,
+  Discrete* discrete, double strength_factor, int target_material_id,
+  const vector<int32_t>& instances)
+{
+  Cell& cell = *model::cells[i_cell];
+
+  if (cell.type_ != Fill::MATERIAL)
+    return;
+
+  for (int j : instances) {
+    int cell_material_idx = cell.material(j);
+    int cell_material_id = model::materials[cell_material_idx]->id();
+    if (target_material_id == C_NONE ||
+        cell_material_id == target_material_id) {
+      int64_t source_region = source_region_offsets_[i_cell] + j;
+      apply_external_source_to_source_region(
+        discrete, strength_factor, source_region);
+    }
+  }
+}
+
+void FlatSourceDomain::apply_external_source_to_cell_and_children(
+  int32_t i_cell, Discrete* discrete, double strength_factor,
+  int32_t target_material_id)
+{
+  Cell& cell = *model::cells[i_cell];
+
+  if (cell.type_ == Fill::MATERIAL) {
+    vector<int> instances(cell.n_instances_);
+    std::iota(instances.begin(), instances.end(), 0);
+    apply_external_source_to_cell_instances(
+      i_cell, discrete, strength_factor, target_material_id, instances);
+  } else if (target_material_id == C_NONE) {
+    std::unordered_map<int32_t, vector<int32_t>> cell_instance_list =
+      cell.get_contained_cells(0, nullptr);
+    for (const auto& pair : cell_instance_list) {
+      int32_t i_child_cell = pair.first;
+      apply_external_source_to_cell_instances(i_child_cell, discrete,
+        strength_factor, target_material_id, pair.second);
+    }
+  }
+}
+
+void FlatSourceDomain::count_external_source_regions()
+{
+#pragma omp parallel for reduction(+ : n_external_source_regions_)
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    float total = 0.f;
+    for (int e = 0; e < negroups_; e++) {
+      int64_t se = sr * negroups_ + e;
+      total += external_source_[se];
+    }
+    if (total != 0.f) {
+      n_external_source_regions_++;
+    }
+  }
+}
+
+void FlatSourceDomain::convert_external_sources()
+{
+  // Loop over external sources
+  for (int es = 0; es < model::external_sources.size(); es++) {
+    Source* s = model::external_sources[es].get();
+    IndependentSource* is = dynamic_cast<IndependentSource*>(s);
+    Discrete* energy = dynamic_cast<Discrete*>(is->energy());
+    const std::unordered_set<int32_t>& domain_ids = is->domain_ids();
+
+    double strength_factor = is->strength();
+
+    if (is->domain_type() == Source::DomainType::MATERIAL) {
+      for (int32_t material_id : domain_ids) {
+        for (int i_cell = 0; i_cell < model::cells.size(); i_cell++) {
+          apply_external_source_to_cell_and_children(
+            i_cell, energy, strength_factor, material_id);
+        }
+      }
+    } else if (is->domain_type() == Source::DomainType::CELL) {
+      for (int32_t cell_id : domain_ids) {
+        int32_t i_cell = model::cell_map[cell_id];
+        apply_external_source_to_cell_and_children(
+          i_cell, energy, strength_factor, C_NONE);
+      }
+    } else if (is->domain_type() == Source::DomainType::UNIVERSE) {
+      for (int32_t universe_id : domain_ids) {
+        int32_t i_universe = model::universe_map[universe_id];
+        Universe& universe = *model::universes[i_universe];
+        for (int32_t i_cell : universe.cells_) {
+          apply_external_source_to_cell_and_children(
+            i_cell, energy, strength_factor, C_NONE);
+        }
+      }
+    }
+  } // End loop over external sources
+
+  // Temperature and angle indices, if using multiple temperature
+  // data sets and/or anisotropic data sets.
+  // TODO: Currently assumes we are only using single temp/single angle data.
+  const int t = 0;
+  const int a = 0;
+
+// Divide the fixed source term by sigma t (to save time when applying each
+// iteration)
+#pragma omp parallel for
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    int material = material_[sr];
+    for (int e = 0; e < negroups_; e++) {
+      float sigma_t = data::mg.macro_xs_[material].get_xs(
+        MgxsType::TOTAL, e, nullptr, nullptr, nullptr, t, a);
+      external_source_[sr * negroups_ + e] /= sigma_t;
+    }
+  }
+}
+
 } // namespace openmc
--- a/src/random_ray/random_ray_simulation.cpp
+++ b/src/random_ray/random_ray_simulation.cpp
@ -111,13 +111,6 @@ void validate_random_ray_inputs()
    }
  }

-  // Validate solver mode
-  ///////////////////////////////////////////////////////////////////
-  if (settings::run_mode == RunMode::FIXED_SOURCE) {
-    fatal_error(
-      "Invalid run mode. Fixed source not yet supported in random ray mode.");
-  }
-
  // Validate ray source
  ///////////////////////////////////////////////////////////////////

@ -125,8 +118,8 @@ void validate_random_ray_inputs()
  IndependentSource* is =
    dynamic_cast<IndependentSource*>(RandomRay::ray_source_.get());
  if (!is) {
-    fatal_error(
-      "Invalid ray source definition. Ray source must be IndependentSource.");
+    fatal_error("Invalid ray source definition. Ray source must provided and "
+                "be of type IndependentSource.");
  }

  // Check for box source
@ -134,24 +127,68 @@ void validate_random_ray_inputs()
  SpatialBox* sb = dynamic_cast<SpatialBox*>(space_dist);
  if (!sb) {
    fatal_error(
-      "Invalid source definition -- only box sources are allowed in random "
-      "ray "
-      "mode. If no source is specified, OpenMC default is an isotropic point "
-      "source at the origin, which is invalid in random ray mode.");
+      "Invalid ray source definition -- only box sources are allowed.");
  }

  // Check that box source is not restricted to fissionable areas
  if (sb->only_fissionable()) {
-    fatal_error("Invalid source definition -- fissionable spatial distribution "
-                "not allowed for random ray source.");
+    fatal_error(
+      "Invalid ray source definition -- fissionable spatial distribution "
+      "not allowed.");
  }

  // Check for isotropic source
  UnitSphereDistribution* angle_dist = is->angle();
  Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);
  if (!id) {
-    fatal_error("Invalid source definition -- only isotropic sources are "
-                "allowed for random ray source.");
+    fatal_error("Invalid ray source definition -- only isotropic sources are "
+                "allowed.");
+  }
+
+  // Validate external sources
+  ///////////////////////////////////////////////////////////////////
+  if (settings::run_mode == RunMode::FIXED_SOURCE) {
+    if (model::external_sources.size() < 1) {
+      fatal_error("Must provide a particle source (in addition to ray source) "
+                  "in fixed source random ray mode.");
+    }
+
+    for (int i = 0; i < model::external_sources.size(); i++) {
+      Source* s = model::external_sources[i].get();
+
+      // Check for independent source
+      IndependentSource* is = dynamic_cast<IndependentSource*>(s);
+
+      if (!is) {
+        fatal_error(
+          "Only IndependentSource external source types are allowed in "
+          "random ray mode");
+      }
+
+      // Check for isotropic source
+      UnitSphereDistribution* angle_dist = is->angle();
+      Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);
+      if (!id) {
+        fatal_error(
+          "Invalid source definition -- only isotropic external sources are "
+          "allowed in random ray mode.");
+      }
+
+      // Validate that a domain ID was specified
+      if (is->domain_ids().size() == 0) {
+        fatal_error("Fixed sources must be specified by domain "
+                    "id (cell, material, or universe) in random ray mode.");
+      }
+
+      // Check that a discrete energy distribution was used
+      Distribution* d = is->energy();
+      Discrete* dd = dynamic_cast<Discrete*>(d);
+      if (!dd) {
+        fatal_error(
+          "Only discrete (multigroup) energy distributions are allowed for "
+          "external sources in random ray mode.");
+      }
+    }
  }

  // Validate plotting files
@ -208,6 +245,12 @@ RandomRaySimulation::RandomRaySimulation()

 void RandomRaySimulation::simulate()
 {
+  if (settings::run_mode == RunMode::FIXED_SOURCE) {
+    // Transfer external source user inputs onto random ray source regions
+    domain_.convert_external_sources();
+    domain_.count_external_source_regions();
+  }
+
  // Random ray power iteration loop
  while (simulation::current_batch < settings::n_batches) {

@ -249,11 +292,13 @@ void RandomRaySimulation::simulate()
    // Add source to scalar flux, compute number of FSR hits
    int64_t n_hits = domain_.add_source_to_scalar_flux();

-    // Compute random ray k-eff
-    k_eff_ = domain_.compute_k_eff(k_eff_);
+    if (settings::run_mode == RunMode::EIGENVALUE) {
+      // Compute random ray k-eff
+      k_eff_ = domain_.compute_k_eff(k_eff_);

-    // Store random ray k-eff into OpenMC's native k-eff variable
-    global_tally_tracklength = k_eff_;
+      // Store random ray k-eff into OpenMC's native k-eff variable
+      global_tally_tracklength = k_eff_;
+    }

    // Execute all tallying tasks, if this is an active batch
    if (simulation::current_batch > settings::n_inactive && mpi::master) {
@ -302,7 +347,7 @@ void RandomRaySimulation::output_simulation_results() const
  if (mpi::master) {
    print_results_random_ray(total_geometric_intersections_,
      avg_miss_rate_ / settings::n_batches, negroups_,
-      domain_.n_source_regions_);
+      domain_.n_source_regions_, domain_.n_external_source_regions_);
    if (model::plots.size() > 0) {
      domain_.output_to_vtk();
    }
@ -339,7 +384,7 @@ void RandomRaySimulation::instability_check(
 // Print random ray simulation results
 void RandomRaySimulation::print_results_random_ray(
  uint64_t total_geometric_intersections, double avg_miss_rate, int negroups,
-  int64_t n_source_regions) const
+  int64_t n_source_regions, int64_t n_external_source_regions) const
 {
  using namespace simulation;

@ -355,6 +400,8 @@ void RandomRaySimulation::print_results_random_ray(
    fmt::print(
      " Total Iterations                  = {}\n", settings::n_batches);
    fmt::print(" Flat Source Regions (FSRs)        = {}\n", n_source_regions);
+    fmt::print(
+      " FSRs Containing External Sources  = {}\n", n_external_source_regions);
    fmt::print(" Total Geometric Intersections     = {:.4e}\n",
      static_cast<double>(total_geometric_intersections));
    fmt::print("   Avg per Iteration               = {:.4e}\n",
@ -387,7 +434,7 @@ void RandomRaySimulation::print_results_random_ray(
    show_time("Time per integration", time_per_integration);
  }

-  if (settings::verbosity >= 4) {
+  if (settings::verbosity >= 4 && settings::run_mode == RunMode::EIGENVALUE) {
    header("Results", 4);
    fmt::print(" k-effective                       = {:.5f} +/- {:.5f}\n",
      simulation::keff, simulation::keff_std);
--- a/src/settings.cpp
+++ b/src/settings.cpp
@ -191,7 +191,8 @@ void get_run_parameters(pugi::xml_node node_base)
  }

  // Get number of inactive batches
-  if (run_mode == RunMode::EIGENVALUE) {
+  if (run_mode == RunMode::EIGENVALUE ||
+      solver_type == SolverType::RANDOM_RAY) {
    if (check_for_node(node_base, "inactive")) {
      n_inactive = std::stoi(get_node_value(node_base, "inactive"));
    }
@ -524,8 +525,9 @@ void read_settings_xml(pugi::xml_node root)
  }

  // If no source specified, default to isotropic point source at origin with
-  // Watt spectrum
-  if (model::external_sources.empty()) {
+  // Watt spectrum. No default source is needed in random ray mode.
+  if (model::external_sources.empty() &&
+      settings::solver_type != SolverType::RANDOM_RAY) {
    double T[] {0.0};
    double p[] {1.0};
    model::external_sources.push_back(make_unique<IndependentSource>(
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@ -137,7 +137,11 @@ int openmc_simulation_init()
  // Display header
  if (mpi::master) {
    if (settings::run_mode == RunMode::FIXED_SOURCE) {
-      header("FIXED SOURCE TRANSPORT SIMULATION", 3);
+      if (settings::solver_type == SolverType::MONTE_CARLO) {
+        header("FIXED SOURCE TRANSPORT SIMULATION", 3);
+      } else if (settings::solver_type == SolverType::RANDOM_RAY) {
+        header("FIXED SOURCE TRANSPORT SIMULATION (RANDOM RAY SOLVER)", 3);
+      }
    } else if (settings::run_mode == RunMode::EIGENVALUE) {
      if (settings::solver_type == SolverType::MONTE_CARLO) {
        header("K EIGENVALUE SIMULATION", 3);
@ -343,7 +347,13 @@ void initialize_batch()
  ++simulation::current_batch;

  if (settings::run_mode == RunMode::FIXED_SOURCE) {
-    write_message(6, "Simulating batch {}", simulation::current_batch);
+    if (settings::solver_type == SolverType::RANDOM_RAY &&
+        simulation::current_batch < settings::n_inactive + 1) {
+      write_message(
+        6, "Simulating batch {:<4} (inactive)", simulation::current_batch);
+    } else {
+      write_message(6, "Simulating batch {}", simulation::current_batch);
+    }
  }

  // Reset total starting particle weight used for normalizing tallies
--- a/tests/regression_tests/random_ray_basic/results_true.dat
+++ b/tests/regression_tests/random_ray_basic/results_true.dat
@ -1,171 +1,171 @@
 k-combined:
-8.400322E-01 8.023349E-03
+8.400322E-01 8.023350E-03
 tally 1:
-1.260220E+00
-3.179889E-01
-1.484289E-01
-4.411066E-03
-3.612463E-01
-2.612843E-02
-7.086707E-01
-1.006119E-01
-3.342483E-02
-2.238499E-04
-8.134936E-02
-1.325949E-03
-4.194328E-01
-3.558669E-02
-4.287776E-03
-3.717447E-06
-1.043559E-02
-2.201986E-05
-5.878720E-01
-7.045887E-02
-6.147757E-03
-7.701173E-06
-1.496241E-02
-4.561699E-05
-1.768113E+00
-6.356917E-01
-6.513486E-03
-8.628535E-06
-1.585272E-02
-5.111136E-05
-5.063704E+00
-5.152401E+00
-2.440293E-03
-1.196869E-06
-6.038334E-03
-7.328193E-06
-3.253717E+00
-2.117655E+00
-1.389120E-02
-3.859385E-05
-3.863767E-02
-2.985798E-04
-1.876994E+00
-7.046366E-01
+5.086560E+00
+5.180937E+00
+1.885166E+00
+7.115505E-01
+4.588117E+00
+4.214785E+00
+2.860401E+00
+1.639329E+00
+4.245221E-01
+3.610930E-02
+1.033202E+00
+2.138892E-01
+1.692631E+00
+5.793967E-01
+5.445818E-02
+5.996625E-04
+1.325403E-01
+3.552030E-03
+2.372249E+00
+1.146944E+00
+7.808143E-02
+1.242279E-03
+1.900346E-01
+7.358492E-03
+7.134949E+00
+1.034824E+01
+8.272648E-02
+1.391872E-03
+2.013422E-01
+8.244790E-03
+2.043539E+01
+8.389902E+01
+3.099367E-02
+1.930673E-04
+7.669167E-02
+1.182113E-03
+1.313212E+01
+3.449537E+01
+1.764293E-01
+6.225587E-03
+4.907293E-01
+4.816401E-02
+7.567717E+00
+1.145439E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-8.390875E-01
-1.408791E-01
+3.383194E+00
+2.290469E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-4.513839E-01
-4.139640E-02
+1.819673E+00
+6.726159E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-6.682186E-01
-9.116003E-02
+2.693683E+00
+1.480961E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.849034E+00
-6.944337E-01
+7.453759E+00
+1.128171E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-4.523425E+00
-4.112118E+00
+1.823561E+01
+6.681652E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.821432E+00
-1.592568E+00
+1.137517E+01
+2.588512E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.159618E+00
-2.694138E-01
-1.354028E-01
-3.672094E-03
-3.295432E-01
-2.175122E-02
-6.880334E-01
-9.491215E-02
-3.234611E-02
-2.097782E-04
-7.872396E-02
-1.242596E-03
-4.184841E-01
-3.536436E-02
-4.274305E-03
-3.687209E-06
-1.040280E-02
-2.184074E-05
-5.810180E-01
-6.872944E-02
-6.060273E-03
-7.476500E-06
-1.474949E-02
-4.428617E-05
-1.782580E+00
-6.457892E-01
-6.552384E-03
-8.730345E-06
-1.594739E-02
-5.171444E-05
-5.278155E+00
-5.596601E+00
-2.546878E-03
-1.303010E-06
-6.302072E-03
-7.978072E-06
-3.420419E+00
-2.340454E+00
-1.465798E-02
-4.299061E-05
-4.077042E-02
-3.325951E-04
-1.279417E+00
-3.278133E-01
-1.509073E-01
-4.561836E-03
-3.672782E-01
-2.702150E-02
-7.212777E-01
-1.042487E-01
-3.411552E-02
-2.332877E-04
-8.303035E-02
-1.381852E-03
-4.269473E-01
-3.685202E-02
-4.378540E-03
-3.872997E-06
-1.065649E-02
-2.294124E-05
-5.973530E-01
-7.266946E-02
-6.260881E-03
-7.976490E-06
-1.523773E-02
-4.724780E-05
-1.795373E+00
-6.547440E-01
-6.635941E-03
-8.945067E-06
-1.615075E-02
-5.298634E-05
-5.161876E+00
-5.353441E+00
-2.505311E-03
-1.261399E-06
-6.199218E-03
-7.723296E-06
-3.344042E+00
-2.236603E+00
-1.443089E-02
-4.166228E-05
-4.013879E-02
-3.223186E-04
+4.601917E+00
+4.242626E+00
+1.719723E+00
+5.923467E-01
+4.185463E+00
+3.508696E+00
+2.730305E+00
+1.494324E+00
+4.108214E-01
+3.383938E-02
+9.998572E-01
+2.004436E-01
+1.660852E+00
+5.570433E-01
+5.428709E-02
+5.947848E-04
+1.321239E-01
+3.523137E-03
+2.306069E+00
+1.082856E+00
+7.697032E-02
+1.206037E-03
+1.873304E-01
+7.143816E-03
+7.075194E+00
+1.017519E+01
+8.322052E-02
+1.408295E-03
+2.025446E-01
+8.342072E-03
+2.094832E+01
+8.816716E+01
+3.234739E-02
+2.101889E-04
+8.004135E-02
+1.286945E-03
+1.357413E+01
+3.685984E+01
+1.861680E-01
+6.934828E-03
+5.178169E-01
+5.365103E-02
+5.072150E+00
+5.151431E+00
+1.916644E+00
+7.358713E-01
+4.664727E+00
+4.358847E+00
+2.859465E+00
+1.638250E+00
+4.332944E-01
+3.763171E-02
+1.054552E+00
+2.229070E-01
+1.693008E+00
+5.796672E-01
+5.561096E-02
+6.247543E-04
+1.353459E-01
+3.700658E-03
+2.368860E+00
+1.143296E+00
+7.951820E-02
+1.286690E-03
+1.935314E-01
+7.621558E-03
+7.119587E+00
+1.030023E+01
+8.428176E-02
+1.442932E-03
+2.051275E-01
+8.547244E-03
+2.046758E+01
+8.418768E+01
+3.181946E-02
+2.034766E-04
+7.873502E-02
+1.245847E-03
+1.325834E+01
+3.515919E+01
+1.832838E-01
+6.720556E-03
+5.097947E-01
+5.199331E-02
--- a/tests/regression_tests/random_ray_fixed_source/init.py
+++ b/tests/regression_tests/random_ray_fixed_source/init.py
--- a/tests/regression_tests/random_ray_fixed_source/cell/inputs_true.dat
+++ b/tests/regression_tests/random_ray_fixed_source/cell/inputs_true.dat
@ -0,0 +1,204 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="source">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="source"/>
+    </material>
+    <material id="2" name="void">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="void"/>
+    </material>
+    <material id="3" name="shield">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="shield"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="infinite source region" universe="1"/>
+    <cell id="2" material="2" name="infinite void region" universe="2"/>
+    <cell id="3" material="3" name="infinite shield region" universe="3"/>
+    <cell fill="4" id="4" name="source lattice cell" universe="7"/>
+    <cell fill="5" id="5" name="void lattice cell" universe="8"/>
+    <cell fill="6" id="6" name="shield lattice cell" universe="9"/>
+    <cell fill="10" id="7" name="dogleg lattice cell" universe="11"/>
+    <cell fill="11" id="8" name="full domain" region="1 -2 3 -4 5 -6" universe="12"/>
+    <lattice id="4">
+      <pitch>5.0 5.0 5.0</pitch>
+      <dimension>2 2 2</dimension>
+      <lower_left>-5.0 -5.0 -5.0</lower_left>
+      <universes>
+1 1 
+1 1 
+
+1 1 
+1 1 </universes>
+    </lattice>
+    <lattice id="5">
+      <pitch>5.0 5.0 5.0</pitch>
+      <dimension>2 2 2</dimension>
+      <lower_left>-5.0 -5.0 -5.0</lower_left>
+      <universes>
+2 2 
+2 2 
+
+2 2 
+2 2 </universes>
+    </lattice>
+    <lattice id="6">
+      <pitch>5.0 5.0 5.0</pitch>
+      <dimension>2 2 2</dimension>
+      <lower_left>-5.0 -5.0 -5.0</lower_left>
+      <universes>
+3 3 
+3 3 
+
+3 3 
+3 3 </universes>
+    </lattice>
+    <lattice id="10">
+      <pitch>10.0 10.0 10.0</pitch>
+      <dimension>6 10 6</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <universes>
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+8 8 8 8 9 9 
+8 9 9 9 9 9 
+8 9 9 9 9 9 
+8 9 9 9 9 9 
+8 9 9 9 9 9 
+7 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 8 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 8 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 </universes>
+    </lattice>
+    <surface boundary="reflective" coeffs="0.0" id="1" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="60.0" id="2" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="3" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="100.0" id="4" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="5" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="60.0" id="6" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <energy type="discrete">
+        <parameters>100.0 1.0</parameters>
+      </energy>
+      <constraints>
+        <domain_type>cell</domain_type>
+        <domain_ids>4</domain_ids>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>400.0</distance_active>
+      <distance_inactive>100.0</distance_inactive>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>0.0 0.0 0.0 60.0 100.0 60.0</parameters>
+        </space>
+      </source>
+    </random_ray>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>1 10 1</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <upper_right>10.0 100.0 10.0</upper_right>
+    </mesh>
+    <mesh id="2">
+      <dimension>6 1 1</dimension>
+      <lower_left>0.0 50.0 0.0</lower_left>
+      <upper_right>60.0 60.0 10.0</upper_right>
+    </mesh>
+    <mesh id="3">
+      <dimension>6 1 1</dimension>
+      <lower_left>0.0 90.0 30.0</lower_left>
+      <upper_right>60.0 100.0 40.0</upper_right>
+    </mesh>
+    <filter id="1" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>2</bins>
+    </filter>
+    <filter id="3" type="mesh">
+      <bins>3</bins>
+    </filter>
+    <tally id="1" name="Case 3A">
+      <filters>1</filters>
+      <scores>flux</scores>
+      <estimator>analog</estimator>
+    </tally>
+    <tally id="2" name="Case 3B">
+      <filters>2</filters>
+      <scores>flux</scores>
+      <estimator>analog</estimator>
+    </tally>
+    <tally id="3" name="Case 3C">
+      <filters>3</filters>
+      <scores>flux</scores>
+      <estimator>analog</estimator>
+    </tally>
+  </tallies>
+</model>
--- a/tests/regression_tests/random_ray_fixed_source/cell/results_true.dat
+++ b/tests/regression_tests/random_ray_fixed_source/cell/results_true.dat
@ -0,0 +1,47 @@
+tally 1:
+3.751048E+01
+2.841797E+02
+9.930788E+00
+1.988180E+01
+3.781121E+00
+3.005165E+00
+2.383139E+00
+1.232560E+00
+1.561884E+00
+6.440577E-01
+1.089787E+00
+3.724896E-01
+6.608456E-01
+1.285592E-01
+2.372611E-01
+1.601299E-02
+7.814803E-02
+1.765829E-03
+2.862108E-02
+2.460129E-04
+tally 2:
+1.089787E+00
+3.724896E-01
+3.767926E-01
+3.724399E-02
+8.614121E-02
+1.526889E-03
+3.610725E-02
+2.629885E-04
+1.466261E-02
+4.536997E-05
+4.653106E-03
+4.381672E-06
+tally 3:
+1.617918E-03
+6.317049E-07
+1.161473E-03
+2.789553E-07
+1.198879E-03
+3.189531E-07
+1.031737E-03
+2.207381E-07
+5.466329E-04
+6.166808E-08
+2.146062E-04
+9.937520E-09
--- a/tests/regression_tests/random_ray_fixed_source/material/inputs_true.dat
+++ b/tests/regression_tests/random_ray_fixed_source/material/inputs_true.dat
@ -0,0 +1,204 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="source">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="source"/>
+    </material>
+    <material id="2" name="void">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="void"/>
+    </material>
+    <material id="3" name="shield">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="shield"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="infinite source region" universe="1"/>
+    <cell id="2" material="2" name="infinite void region" universe="2"/>
+    <cell id="3" material="3" name="infinite shield region" universe="3"/>
+    <cell fill="4" id="4" name="source lattice cell" universe="7"/>
+    <cell fill="5" id="5" name="void lattice cell" universe="8"/>
+    <cell fill="6" id="6" name="shield lattice cell" universe="9"/>
+    <cell fill="10" id="7" name="dogleg lattice cell" universe="11"/>
+    <cell fill="11" id="8" name="full domain" region="1 -2 3 -4 5 -6" universe="12"/>
+    <lattice id="4">
+      <pitch>5.0 5.0 5.0</pitch>
+      <dimension>2 2 2</dimension>
+      <lower_left>-5.0 -5.0 -5.0</lower_left>
+      <universes>
+1 1 
+1 1 
+
+1 1 
+1 1 </universes>
+    </lattice>
+    <lattice id="5">
+      <pitch>5.0 5.0 5.0</pitch>
+      <dimension>2 2 2</dimension>
+      <lower_left>-5.0 -5.0 -5.0</lower_left>
+      <universes>
+2 2 
+2 2 
+
+2 2 
+2 2 </universes>
+    </lattice>
+    <lattice id="6">
+      <pitch>5.0 5.0 5.0</pitch>
+      <dimension>2 2 2</dimension>
+      <lower_left>-5.0 -5.0 -5.0</lower_left>
+      <universes>
+3 3 
+3 3 
+
+3 3 
+3 3 </universes>
+    </lattice>
+    <lattice id="10">
+      <pitch>10.0 10.0 10.0</pitch>
+      <dimension>6 10 6</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <universes>
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+8 8 8 8 9 9 
+8 9 9 9 9 9 
+8 9 9 9 9 9 
+8 9 9 9 9 9 
+8 9 9 9 9 9 
+7 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 8 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 8 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 </universes>
+    </lattice>
+    <surface boundary="reflective" coeffs="0.0" id="1" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="60.0" id="2" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="3" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="100.0" id="4" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="5" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="60.0" id="6" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <energy type="discrete">
+        <parameters>100.0 1.0</parameters>
+      </energy>
+      <constraints>
+        <domain_type>material</domain_type>
+        <domain_ids>1</domain_ids>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>400.0</distance_active>
+      <distance_inactive>100.0</distance_inactive>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>0.0 0.0 0.0 60.0 100.0 60.0</parameters>
+        </space>
+      </source>
+    </random_ray>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>1 10 1</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <upper_right>10.0 100.0 10.0</upper_right>
+    </mesh>
+    <mesh id="2">
+      <dimension>6 1 1</dimension>
+      <lower_left>0.0 50.0 0.0</lower_left>
+      <upper_right>60.0 60.0 10.0</upper_right>
+    </mesh>
+    <mesh id="3">
+      <dimension>6 1 1</dimension>
+      <lower_left>0.0 90.0 30.0</lower_left>
+      <upper_right>60.0 100.0 40.0</upper_right>
+    </mesh>
+    <filter id="1" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>2</bins>
+    </filter>
+    <filter id="3" type="mesh">
+      <bins>3</bins>
+    </filter>
+    <tally id="1" name="Case 3A">
+      <filters>1</filters>
+      <scores>flux</scores>
+      <estimator>analog</estimator>
+    </tally>
+    <tally id="2" name="Case 3B">
+      <filters>2</filters>
+      <scores>flux</scores>
+      <estimator>analog</estimator>
+    </tally>
+    <tally id="3" name="Case 3C">
+      <filters>3</filters>
+      <scores>flux</scores>
+      <estimator>analog</estimator>
+    </tally>
+  </tallies>
+</model>
--- a/tests/regression_tests/random_ray_fixed_source/material/results_true.dat
+++ b/tests/regression_tests/random_ray_fixed_source/material/results_true.dat
@ -0,0 +1,47 @@
+tally 1:
+3.751047E+01
+2.841797E+02
+9.930788E+00
+1.988181E+01
+3.781121E+00
+3.005165E+00
+2.383139E+00
+1.232560E+00
+1.561884E+00
+6.440577E-01
+1.089787E+00
+3.724896E-01
+6.608456E-01
+1.285592E-01
+2.372611E-01
+1.601299E-02
+7.814803E-02
+1.765829E-03
+2.862108E-02
+2.460129E-04
+tally 2:
+1.089787E+00
+3.724896E-01
+3.767925E-01
+3.724398E-02
+8.614120E-02
+1.526889E-03
+3.610725E-02
+2.629885E-04
+1.466261E-02
+4.536997E-05
+4.653106E-03
+4.381672E-06
+tally 3:
+1.617918E-03
+6.317049E-07
+1.161473E-03
+2.789553E-07
+1.198879E-03
+3.189531E-07
+1.031737E-03
+2.207381E-07
+5.466329E-04
+6.166809E-08
+2.146062E-04
+9.937520E-09
--- a/tests/regression_tests/random_ray_fixed_source/test.py
+++ b/tests/regression_tests/random_ray_fixed_source/test.py
@ -0,0 +1,339 @@
+import os
+
+import numpy as np
+import openmc
+from openmc.utility_funcs import change_directory
+import pytest
+
+from tests.testing_harness import TolerantPyAPITestHarness
+
+def fill_3d_list(n, val):
+    """
+    Generates a 3D list of dimensions nxnxn filled with copies of val.
+
+    Parameters:
+    n (int): The dimension of the 3D list.
+    val (any): The value to fill the 3D list with.
+
+    Returns:
+    list: A 3D list of dimensions nxnxn filled with val.
+    """
+    return [[[val for _ in range(n)] for _ in range(n)] for _ in range(n)]
+
+
+class MGXSTestHarness(TolerantPyAPITestHarness):
+    def _cleanup(self):
+        super()._cleanup()
+        f = 'mgxs.h5'
+        if os.path.exists(f):
+            os.remove(f)
+
+def create_random_ray_model(domain_type):
+    openmc.reset_auto_ids()
+    ###############################################################################
+    # Create multigroup data
+
+    # Instantiate the energy group data
+    ebins = [1e-5, 20.0e6]
+    groups = openmc.mgxs.EnergyGroups(group_edges=ebins)
+
+    # High scattering ratio means system is all scattering
+    # Low means fully absorbing
+    scattering_ratio = 0.5
+
+    source_total_xs = 0.1
+    source_mat_data = openmc.XSdata('source', groups)
+    source_mat_data.order = 0
+    source_mat_data.set_total([source_total_xs])
+    source_mat_data.set_absorption([source_total_xs * (1.0 - scattering_ratio)])
+    source_mat_data.set_scatter_matrix(np.rollaxis(np.array([[[source_total_xs * scattering_ratio]]]),0,3))
+
+    void_total_xs = 1.0e-4
+    void_mat_data = openmc.XSdata('void', groups)
+    void_mat_data.order = 0
+    void_mat_data.set_total([void_total_xs])
+    void_mat_data.set_absorption([void_total_xs * (1.0 - scattering_ratio)])
+    void_mat_data.set_scatter_matrix(np.rollaxis(np.array([[[void_total_xs * scattering_ratio]]]),0,3))
+
+    shield_total_xs = 0.1
+    shield_mat_data = openmc.XSdata('shield', groups)
+    shield_mat_data.order = 0
+    shield_mat_data.set_total([shield_total_xs])
+    shield_mat_data.set_absorption([shield_total_xs * (1.0 - scattering_ratio)])
+    shield_mat_data.set_scatter_matrix(np.rollaxis(np.array([[[shield_total_xs * scattering_ratio]]]),0,3))
+
+    mg_cross_sections_file = openmc.MGXSLibrary(groups)
+    mg_cross_sections_file.add_xsdatas([source_mat_data, void_mat_data, shield_mat_data])
+    mg_cross_sections_file.export_to_hdf5()
+
+    ###############################################################################
+    # Create materials for the problem
+
+    # Instantiate some Macroscopic Data
+    source_data = openmc.Macroscopic('source')
+    void_data   = openmc.Macroscopic('void')
+    shield_data = openmc.Macroscopic('shield')
+
+    # Instantiate some Materials and register the appropriate Macroscopic objects
+    source_mat = openmc.Material(name='source')
+    source_mat.set_density('macro', 1.0)
+    source_mat.add_macroscopic(source_data)
+
+    void_mat = openmc.Material(name='void')
+    void_mat.set_density('macro', 1.0)
+    void_mat.add_macroscopic(void_data)
+
+    shield_mat = openmc.Material(name='shield')
+    shield_mat.set_density('macro', 1.0)
+    shield_mat.add_macroscopic(shield_data)
+
+    # Instantiate a Materials collection and export to XML
+    materials_file = openmc.Materials([source_mat, void_mat, shield_mat])
+    materials_file.cross_sections = "mgxs.h5"
+
+    ###############################################################################
+    # Define problem geometry
+
+    source_cell = openmc.Cell(fill=source_mat, name='infinite source region')
+    void_cell = openmc.Cell(fill=void_mat, name='infinite void region')
+    shield_cell = openmc.Cell(fill=shield_mat, name='infinite shield region')
+
+    sub = openmc.Universe()
+    sub.add_cells([source_cell])
+
+    vub = openmc.Universe()
+    vub.add_cells([void_cell])
+
+    aub = openmc.Universe()
+    aub.add_cells([shield_cell])
+
+    # n controls the dimension of subdivision within each outer lattice element
+    # E.g., n = 10 results in 1cm cubic FSRs
+    n = 2
+    delta = 10.0 / n
+    ll = [-5.0, -5.0, -5.0]
+    pitch = [delta, delta, delta]
+
+    source_lattice = openmc.RectLattice()
+    source_lattice.lower_left = ll
+    source_lattice.pitch = pitch
+    source_lattice.universes = fill_3d_list(n, sub)
+
+    void_lattice = openmc.RectLattice()
+    void_lattice.lower_left = ll
+    void_lattice.pitch = pitch
+    void_lattice.universes = fill_3d_list(n, vub)
+
+    shield_lattice = openmc.RectLattice()
+    shield_lattice.lower_left = ll
+    shield_lattice.pitch = pitch
+    shield_lattice.universes = fill_3d_list(n, aub)
+
+    source_lattice_cell = openmc.Cell(fill=source_lattice, name='source lattice cell')
+    su = openmc.Universe()
+    su.add_cells([source_lattice_cell])
+
+    void_lattice_cell = openmc.Cell(fill=void_lattice, name='void lattice cell')
+    vu = openmc.Universe()
+    vu.add_cells([void_lattice_cell])
+
+    shield_lattice_cell = openmc.Cell(fill=shield_lattice, name='shield lattice cell')
+    au = openmc.Universe()
+    au.add_cells([shield_lattice_cell])
+
+    z_base = [
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [vu, vu, vu, vu, au, au],
+            [vu, au, au, au, au, au],
+            [vu, au, au, au, au, au],
+            [vu, au, au, au, au, au],
+            [vu, au, au, au, au, au],
+            [su, au, au, au, au, au]
+                         ]
+
+    z_col = [
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, vu, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au]
+                         ]
+
+    z_high = [
+            [au, au, au, vu, au, au],
+            [au, au, au, vu, au, au],
+            [au, au, au, vu, au, au],
+            [au, au, au, vu, au, au],
+            [au, au, au, vu, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au]
+                         ]
+
+    z_cap = [
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au],
+            [au, au, au, au, au, au]
+                         ]
+
+    dogleg_pattern = [
+            z_base,
+            z_col,
+            z_col,
+            z_high,
+            z_cap,
+            z_cap
+            ]
+
+    x = 60.0
+    x_dim = 6
+
+    y = 100.0
+    y_dim = 10
+
+    z = 60.0
+    z_dim = 6
+
+    lattice = openmc.RectLattice()
+    lattice.lower_left = [0.0, 0.0, 0.0]
+    lattice.pitch = [x/x_dim, y/y_dim, z/z_dim]
+    lattice.universes = dogleg_pattern
+
+    lattice_cell = openmc.Cell(fill=lattice, name='dogleg lattice cell')
+
+    lattice_uni = openmc.Universe()
+    lattice_uni.add_cells([lattice_cell])
+
+    x_low  = openmc.XPlane(x0=0.0,boundary_type='reflective')
+    x_high = openmc.XPlane(x0=x,boundary_type='vacuum')
+
+    y_low  = openmc.YPlane(y0=0.0,boundary_type='reflective')
+    y_high = openmc.YPlane(y0=y,boundary_type='vacuum')
+
+    z_low  = openmc.ZPlane(z0=0.0,boundary_type='reflective')
+    z_high = openmc.ZPlane(z0=z,boundary_type='vacuum')
+
+    full_domain = openmc.Cell(fill=lattice_uni, region=+x_low & -x_high & +y_low & -y_high & +z_low & -z_high, name='full domain')
+
+    root = openmc.Universe(name='root universe')
+    root.add_cell(full_domain)
+
+    # Create a geometry with the two cells and export to XML
+    geometry = openmc.Geometry(root)
+
+    ###############################################################################
+    # Define problem settings
+
+    # Instantiate a Settings object, set all runtime parameters, and export to XML
+    settings = openmc.Settings()
+    settings.energy_mode = "multi-group"
+    settings.batches = 10
+    settings.inactive = 5
+    settings.particles = 1000
+    settings.run_mode = 'fixed source'
+
+    settings.random_ray['distance_active'] = 400.0
+    settings.random_ray['distance_inactive'] = 100.0
+
+    # Create an initial uniform spatial source for ray integration
+    lower_left = (0.0, 0.0, 0.0)
+    upper_right = (x, y, z)
+    uniform_dist = openmc.stats.Box(lower_left, upper_right)
+    settings.random_ray['ray_source']= openmc.IndependentSource(space=uniform_dist)
+
+    # Create the neutron source in the bottom right of the moderator
+    strengths = [1.0]
+    midpoints = [100.0]
+    energy_distribution = openmc.stats.Discrete(x=midpoints,p=strengths)
+    if domain_type == 'cell':
+        domain = source_lattice_cell
+    elif domain_type == 'material':
+        domain = source_mat
+    elif domain_type == 'universe':
+        domain = sub
+    source = openmc.IndependentSource(
+        energy=energy_distribution,
+        constraints={'domains': [domain]}
+    )
+    settings.source = [source]
+
+    ###############################################################################
+    # Define tallies
+
+    estimator = 'analog'
+
+    # Case 3A
+    mesh_3A = openmc.RegularMesh()
+    mesh_3A.dimension = (1, y_dim, 1)
+    mesh_3A.lower_left = (0.0, 0.0, 0.0)
+    mesh_3A.upper_right = (10.0, y, 10.0)
+    mesh_filter_3A = openmc.MeshFilter(mesh_3A)
+
+    tally_3A = openmc.Tally(name="Case 3A")
+    tally_3A.filters = [mesh_filter_3A]
+    tally_3A.scores = ['flux']
+    tally_3A.estimator = estimator
+
+    # Case 3B
+    mesh_3B = openmc.RegularMesh()
+    mesh_3B.dimension = (x_dim, 1, 1)
+    mesh_3B.lower_left = (0.0, 50.0, 0.0)
+    mesh_3B.upper_right = (x, 60.0, 10.0)
+    mesh_filter_3B = openmc.MeshFilter(mesh_3B)
+
+    tally_3B = openmc.Tally(name="Case 3B")
+    tally_3B.filters = [mesh_filter_3B]
+    tally_3B.scores = ['flux']
+    tally_3B.estimator = estimator
+
+    # Case 3C
+    mesh_3C = openmc.RegularMesh()
+    mesh_3C.dimension = (x_dim, 1, 1)
+    mesh_3C.lower_left = (0.0, 90.0, 30.0)
+    mesh_3C.upper_right = (x, 100.0, 40.0)
+    mesh_filter_3C = openmc.MeshFilter(mesh_3C)
+
+    tally_3C = openmc.Tally(name="Case 3C")
+    tally_3C.filters = [mesh_filter_3C]
+    tally_3C.scores = ['flux']
+    tally_3C.estimator = estimator
+
+    # Instantiate a Tallies collection and export to XML
+    tallies = openmc.Tallies([tally_3A, tally_3B, tally_3C])
+
+    ###############################################################################
+    # Assmble Model
+
+    model = openmc.model.Model()
+    model.geometry = geometry
+    model.materials = materials_file
+    model.settings = settings
+    model.xs_data = mg_cross_sections_file
+    model.tallies = tallies
+
+    return model
+
+@pytest.mark.parametrize("domain_type", ["cell", "material", "universe"])
+def test_random_ray_fixed_source(domain_type):
+    with change_directory(domain_type):
+        model = create_random_ray_model(domain_type)
+
+        harness = MGXSTestHarness('statepoint.10.h5', model)
+        harness.main()
--- a/tests/regression_tests/random_ray_fixed_source/universe/inputs_true.dat
+++ b/tests/regression_tests/random_ray_fixed_source/universe/inputs_true.dat
@ -0,0 +1,204 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="source">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="source"/>
+    </material>
+    <material id="2" name="void">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="void"/>
+    </material>
+    <material id="3" name="shield">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="shield"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="infinite source region" universe="1"/>
+    <cell id="2" material="2" name="infinite void region" universe="2"/>
+    <cell id="3" material="3" name="infinite shield region" universe="3"/>
+    <cell fill="4" id="4" name="source lattice cell" universe="7"/>
+    <cell fill="5" id="5" name="void lattice cell" universe="8"/>
+    <cell fill="6" id="6" name="shield lattice cell" universe="9"/>
+    <cell fill="10" id="7" name="dogleg lattice cell" universe="11"/>
+    <cell fill="11" id="8" name="full domain" region="1 -2 3 -4 5 -6" universe="12"/>
+    <lattice id="4">
+      <pitch>5.0 5.0 5.0</pitch>
+      <dimension>2 2 2</dimension>
+      <lower_left>-5.0 -5.0 -5.0</lower_left>
+      <universes>
+1 1 
+1 1 
+
+1 1 
+1 1 </universes>
+    </lattice>
+    <lattice id="5">
+      <pitch>5.0 5.0 5.0</pitch>
+      <dimension>2 2 2</dimension>
+      <lower_left>-5.0 -5.0 -5.0</lower_left>
+      <universes>
+2 2 
+2 2 
+
+2 2 
+2 2 </universes>
+    </lattice>
+    <lattice id="6">
+      <pitch>5.0 5.0 5.0</pitch>
+      <dimension>2 2 2</dimension>
+      <lower_left>-5.0 -5.0 -5.0</lower_left>
+      <universes>
+3 3 
+3 3 
+
+3 3 
+3 3 </universes>
+    </lattice>
+    <lattice id="10">
+      <pitch>10.0 10.0 10.0</pitch>
+      <dimension>6 10 6</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <universes>
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+8 8 8 8 9 9 
+8 9 9 9 9 9 
+8 9 9 9 9 9 
+8 9 9 9 9 9 
+8 9 9 9 9 9 
+7 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 8 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 8 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 8 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 
+9 9 9 9 9 9 </universes>
+    </lattice>
+    <surface boundary="reflective" coeffs="0.0" id="1" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="60.0" id="2" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="3" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="100.0" id="4" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="5" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="60.0" id="6" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <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>400.0</distance_active>
+      <distance_inactive>100.0</distance_inactive>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>0.0 0.0 0.0 60.0 100.0 60.0</parameters>
+        </space>
+      </source>
+    </random_ray>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>1 10 1</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <upper_right>10.0 100.0 10.0</upper_right>
+    </mesh>
+    <mesh id="2">
+      <dimension>6 1 1</dimension>
+      <lower_left>0.0 50.0 0.0</lower_left>
+      <upper_right>60.0 60.0 10.0</upper_right>
+    </mesh>
+    <mesh id="3">
+      <dimension>6 1 1</dimension>
+      <lower_left>0.0 90.0 30.0</lower_left>
+      <upper_right>60.0 100.0 40.0</upper_right>
+    </mesh>
+    <filter id="1" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>2</bins>
+    </filter>
+    <filter id="3" type="mesh">
+      <bins>3</bins>
+    </filter>
+    <tally id="1" name="Case 3A">
+      <filters>1</filters>
+      <scores>flux</scores>
+      <estimator>analog</estimator>
+    </tally>
+    <tally id="2" name="Case 3B">
+      <filters>2</filters>
+      <scores>flux</scores>
+      <estimator>analog</estimator>
+    </tally>
+    <tally id="3" name="Case 3C">
+      <filters>3</filters>
+      <scores>flux</scores>
+      <estimator>analog</estimator>
+    </tally>
+  </tallies>
+</model>
--- a/tests/regression_tests/random_ray_fixed_source/universe/results_true.dat
+++ b/tests/regression_tests/random_ray_fixed_source/universe/results_true.dat
@ -0,0 +1,47 @@
+tally 1:
+3.751047E+01
+2.841797E+02
+9.930788E+00
+1.988180E+01
+3.781121E+00
+3.005165E+00
+2.383139E+00
+1.232560E+00
+1.561884E+00
+6.440577E-01
+1.089787E+00
+3.724896E-01
+6.608456E-01
+1.285592E-01
+2.372611E-01
+1.601299E-02
+7.814803E-02
+1.765829E-03
+2.862107E-02
+2.460129E-04
+tally 2:
+1.089787E+00
+3.724896E-01
+3.767926E-01
+3.724399E-02
+8.614120E-02
+1.526889E-03
+3.610725E-02
+2.629885E-04
+1.466261E-02
+4.536997E-05
+4.653106E-03
+4.381672E-06
+tally 3:
+1.617918E-03
+6.317049E-07
+1.161473E-03
+2.789553E-07
+1.198879E-03
+3.189531E-07
+1.031737E-03
+2.207381E-07
+5.466329E-04
+6.166808E-08
+2.146062E-04
+9.937520E-09
--- a/tests/regression_tests/random_ray_vacuum/results_true.dat
+++ b/tests/regression_tests/random_ray_vacuum/results_true.dat
@ -1,171 +1,171 @@
 k-combined:
 1.010455E-01 1.585558E-02
 tally 1:
-1.849176E-01
-7.634332E-03
-2.181815E-02
-1.062861E-04
-5.310100E-02
-6.295730E-04
-4.048251E-02
-3.851890E-04
-1.893676E-03
-8.448769E-07
-4.608828E-03
-5.004529E-06
-4.063643E-03
-4.022442E-06
-4.112970E-05
-4.186661E-10
-1.001015E-04
-2.479919E-09
-7.467029E-03
-1.178864E-05
-7.688748E-05
-1.266903E-09
-1.871288E-04
-7.504350E-09
-3.870644E-02
-3.010745E-04
-1.375240E-04
-3.807356E-09
-3.347099E-04
-2.255298E-08
-4.524967E-01
-4.098857E-02
-2.437418E-04
-1.190325E-08
-6.031220E-04
-7.288126E-08
-4.989226E-01
-4.993728E-02
-2.374296E-03
-1.135824E-06
-6.603983E-03
-8.787258E-06
-3.899991E-01
-3.308783E-02
+7.466984E-01
+1.245920E-01
+2.771079E-01
+1.714504E-02
+6.744252E-01
+1.015566E-01
+1.634870E-01
+6.288701E-03
+2.405120E-02
+1.362874E-04
+5.853580E-02
+8.072823E-04
+1.641162E-02
+6.568765E-05
+5.223801E-04
+6.753516E-08
+1.271369E-03
+4.000365E-07
+3.014736E-02
+1.922973E-04
+9.765325E-04
+2.043645E-07
+2.376685E-03
+1.210529E-06
+1.562379E-01
+4.906526E-03
+1.746664E-03
+6.141658E-07
+4.251084E-03
+3.638028E-06
+1.826385E+00
+6.678141E-01
+3.095716E-03
+1.920117E-06
+7.660132E-03
+1.175650E-05
+2.013809E+00
+8.136726E-01
+3.015545E-02
+1.832201E-04
+8.387586E-02
+1.417475E-03
+1.573069E+00
+5.388374E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-7.108982E-02
-1.144390E-03
+2.867849E-01
+1.864798E-02
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-5.295259E-03
-6.352159E-06
+2.136372E-02
+1.035499E-04
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-9.852001E-03
-1.984406E-05
+3.973323E-02
+3.229766E-04
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-4.414391E-02
-3.905201E-04
+1.779974E-01
+6.350613E-03
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.571668E-01
-1.323140E-02
+1.036886E+00
+2.151147E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.752932E-01
-1.517930E-02
+1.109991E+00
+2.468005E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.465446E-01
-4.901884E-03
-1.700791E-02
-6.587674E-05
-4.139385E-02
-3.902131E-04
-3.424032E-02
-2.841338E-04
-1.598985E-03
-6.216312E-07
-3.891610E-03
-3.682159E-06
-4.067582E-03
-4.209468E-06
-4.152829E-05
-4.494649E-10
-1.010715E-04
-2.662352E-09
-7.526712E-03
-1.225032E-05
-7.769969E-05
-1.328443E-09
-1.891055E-04
-7.868877E-09
-4.008649E-02
-3.246821E-04
-1.417944E-04
-4.070719E-09
-3.451035E-04
-2.411301E-08
-4.859902E-01
-4.747592E-02
-2.606214E-04
-1.369749E-08
-6.448895E-04
-8.386705E-08
-5.475198E-01
-6.061269E-02
-2.625477E-03
-1.405458E-06
-7.302631E-03
-1.087327E-05
-1.909660E-01
-8.147906E-03
-2.269063E-02
-1.149570E-04
-5.522446E-02
-6.809342E-04
-4.196583E-02
-4.141620E-04
-1.980406E-03
-9.227119E-07
-4.819913E-03
-5.465576E-06
-4.247004E-03
-4.420116E-06
-4.341806E-05
-4.691518E-10
-1.056709E-04
-2.778965E-09
-7.742814E-03
-1.272112E-05
-8.039606E-05
-1.389209E-09
-1.956679E-04
-8.228817E-09
-3.982370E-02
-3.190931E-04
-1.427171E-04
-4.103942E-09
-3.473492E-04
-2.430981E-08
-4.849535E-01
-4.707014E-02
-2.678327E-04
-1.438540E-08
-6.627333E-04
-8.807897E-08
-5.493457E-01
-6.069440E-02
-2.717400E-03
-1.501450E-06
-7.558312E-03
-1.161591E-05
+5.813217E-01
+7.704274E-02
+2.160141E-01
+1.062660E-02
+5.257350E-01
+6.294540E-02
+1.358032E-01
+4.462239E-03
+2.030839E-02
+1.002755E-04
+4.942655E-02
+5.939703E-04
+1.612888E-02
+6.604109E-05
+5.274425E-04
+7.250333E-08
+1.283689E-03
+4.294649E-07
+2.985723E-02
+1.925418E-04
+9.868482E-04
+2.142916E-07
+2.401791E-03
+1.269331E-06
+1.590678E-01
+5.110940E-03
+1.800903E-03
+6.566490E-07
+4.383091E-03
+3.889678E-06
+1.928611E+00
+7.475889E-01
+3.310101E-03
+2.209547E-06
+8.190613E-03
+1.352862E-05
+2.172777E+00
+9.544513E-01
+3.334566E-02
+2.267149E-04
+9.274926E-02
+1.753971E-03
+7.566911E-01
+1.278105E-01
+2.881892E-01
+1.854375E-02
+7.013949E-01
+1.098417E-01
+1.662732E-01
+6.495387E-03
+2.515274E-02
+1.488430E-04
+6.121675E-02
+8.816538E-04
+1.682747E-02
+6.933017E-05
+5.514441E-04
+7.567901E-08
+1.342104E-03
+4.482757E-07
+3.068773E-02
+1.997062E-04
+1.021094E-03
+2.240938E-07
+2.485139E-03
+1.327393E-06
+1.578722E-01
+5.013656E-03
+1.812622E-03
+6.620083E-07
+4.411613E-03
+3.921424E-06
+1.922798E+00
+7.400474E-01
+3.401690E-03
+2.320513E-06
+8.417243E-03
+1.420805E-05
+2.178318E+00
+9.546106E-01
+3.451316E-02
+2.421994E-04
+9.599661E-02
+1.873766E-03