diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index 6d82929d1..8094a2e2f 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -131,27 +131,38 @@ extern "C" {
   int openmc_zernike_filter_set_params(int32_t index, const double* x,
                                        const double* y, const double* r);
 
+  //! Sets the mesh and energy grid for CMFD reweight
+  //! \param[in] meshtyally_id id of CMFD Mesh Tally
+  //! \param[in] cmfd_indices indices storing spatial and energy dimensions of CMFD problem
+  //! \param[in] norm CMFD normalization factor
+  extern "C" void openmc_initialize_mesh_egrid(const int meshtally_id, const int* cmfd_indices,
+                                               const double norm);
+
+  //! Sets the mesh and energy grid for CMFD reweight
+  //! \param[in] feedback whether or not to run CMFD feedback
+  //! \param[in] cmfd_src computed CMFD source
+  extern "C" void openmc_cmfd_reweight(const bool feedback, const double* cmfd_src);
+
   //! Sets the fixed variables that are used for CMFD linear solver
-  //! \param[in] CSR format index pointer array of loss matrix
-  //! \param[in] length of indptr
-  //! \param[in] CSR format index array of loss matrix
-  //! \param[in] number of non-zero elements in CMFD loss matrix
-  //! \param[in] dimension n of nxn CMFD loss matrix
-  //! \param[in] spectral radius of CMFD matrices and tolerances
-  //! \param[in] indices storing spatial and energy dimensions of CMFD problem
-  //! \param[in] coremap for problem, storing accelerated regions
+  //! \param[in] indptr CSR format index pointer array of loss matrix
+  //! \param[in] len_indptr length of indptr
+  //! \param[in] indices CSR format index array of loss matrix
+  //! \param[in] n_elements number of non-zero elements in CMFD loss matrix
+  //! \param[in] dim dimension n of nxn CMFD loss matrix
+  //! \param[in] spectral spectral radius of CMFD matrices and tolerances
+  //! \param[in] map coremap for problem, storing accelerated regions
+  //! \param[in] use_all_threads whether to use all threads when running CMFD solver
   extern "C" void openmc_initialize_linsolver(const int* indptr, int len_indptr,
                                               const int* indices, int n_elements,
                                               int dim, double spectral,
-                                              const int* cmfd_indices,
                                               const int* map, bool use_all_threads);
 
   //! Runs a Gauss Seidel linear solver to solve CMFD matrix equations
   //! linear solver
-  //! \param[in] CSR format data array of coefficient matrix
-  //! \param[in] right hand side vector
-  //! \param[out] unknown vector
-  //! \param[in] tolerance on final error
+  //! \param[in] A_data CSR format data array of coefficient matrix
+  //! \param[in] b right hand side vector
+  //! \param[out] x unknown vector
+  //! \param[in] tol tolerance on final error
   //! \return number of inner iterations required to reach convergence
   extern "C" int openmc_run_linsolver(const double* A_data, const double* b,
                                       double* x, double tol);
diff --git a/openmc/cmfd.py b/openmc/cmfd.py
index 78901b03b..a99b60619 100644
--- a/openmc/cmfd.py
+++ b/openmc/cmfd.py
@@ -366,8 +366,6 @@ class CMFDRun:
         self._current = None
         self._cmfd_src = None
         self._openmc_src = None
-        self._sourcecounts = None
-        self._weightfactors = None
         self._entropy = []
         self._balance = []
         self._src_cmp = []
@@ -914,7 +912,7 @@ class CMFDRun:
 
         args = temp_loss.indptr, len(temp_loss.indptr), \
             temp_loss.indices, len(temp_loss.indices), n, \
-            self._spectral, self._indices, coremap, self._use_all_threads
+            self._spectral, coremap, self._use_all_threads
         return openmc.lib._dll.openmc_initialize_linsolver(*args)
 
     def _write_cmfd_output(self):
@@ -1171,8 +1169,12 @@ class CMFDRun:
                 # Calculate fission source
                 self._calc_fission_source()
 
-            # Calculate weight factors
-            self._cmfd_reweight()
+            # Calculate weight factors through C++ and manipulate CMFD
+            # source into a 1-D vector that matches C++ array ordering
+            src_flipped = np.flip(self._cmfd_src, axis=3)
+            src_swapped = np.swapaxes(src_flipped, 0, 2)
+            args = self._feedback, src_swapped.flatten()
+            openmc.lib._dll.openmc_cmfd_reweight(*args)
 
         # Stop CMFD timer
         if openmc.lib.master():
@@ -1390,151 +1392,6 @@ class CMFDRun:
         self._src_cmp.append(np.sqrt(1.0 / self._norm
                              * np.sum((self._cmfd_src - self._openmc_src)**2)))
 
-    def _cmfd_reweight(self):
-        """Performs weighting of particles in source bank"""
-        # Get spatial dimensions and energy groups
-        nx, ny, nz, ng = self._indices
-
-        # Count bank site in mesh and reverse due to egrid structured
-        outside = self._count_bank_sites()
-
-        # Check and raise error if source sites exist outside of CMFD mesh
-        if openmc.lib.master() and outside:
-            raise OpenMCError('Source sites outside of the CMFD mesh')
-
-        # Have master compute weight factors, ignore any zeros in
-        # sourcecounts or cmfd_src
-        if openmc.lib.master():
-            # Compute normalization factor
-            norm = np.sum(self._sourcecounts) / np.sum(self._cmfd_src)
-
-            # Define target reshape dimensions for sourcecounts. This
-            # defines how self._sourcecounts is ordered by dimension
-            target_shape = [nz, ny, nx, ng]
-
-            # Reshape sourcecounts to target shape. Swap x and z axes so
-            # that the shape is now [nx, ny, nz, ng]
-            sourcecounts = np.swapaxes(
-                    self._sourcecounts.reshape(target_shape), 0, 2)
-
-            # Flip index of energy dimension
-            sourcecounts = np.flip(sourcecounts, axis=3)
-
-            # Compute weight factors
-            div_condition = np.logical_and(sourcecounts > 0,
-                                           self._cmfd_src > 0)
-            self._weightfactors = (np.divide(self._cmfd_src * norm,
-                                   sourcecounts, where=div_condition,
-                                   out=np.ones_like(self._cmfd_src),
-                                   dtype=np.float32))
-
-        if not self._feedback:
-            return
-
-        # Broadcast weight factors to all procs
-        if have_mpi:
-            self._weightfactors = self._intracomm.bcast(
-                                  self._weightfactors)
-
-        m = openmc.lib.meshes[self._mesh_id]
-        energy = self._egrid
-        ng = self._indices[3]
-
-        # Get locations and energies of all particles in source bank
-        source_xyz = openmc.lib.source_bank()['r']
-        source_energies = openmc.lib.source_bank()['E']
-
-        # Convert xyz location to the CMFD mesh index
-        mesh_ijk = np.floor((source_xyz - m.lower_left)/m.width).astype(int)
-
-        # Determine which energy bin each particle's energy belongs to
-        # Separate into cases bases on where source energies lies on egrid
-        energy_bins = np.zeros(len(source_energies), dtype=int)
-        idx = np.where(source_energies < energy[0])
-        energy_bins[idx] = ng - 1
-        idx = np.where(source_energies > energy[-1])
-        energy_bins[idx] = 0
-        idx = np.where((source_energies >= energy[0]) &
-                       (source_energies <= energy[-1]))
-        energy_bins[idx] = ng - np.digitize(source_energies[idx], energy)
-
-        # Determine weight factor of each particle based on its mesh index
-        # and energy bin and updates its weight
-        openmc.lib.source_bank()['wgt'] *= self._weightfactors[
-                mesh_ijk[:,0], mesh_ijk[:,1], mesh_ijk[:,2], energy_bins]
-
-        if openmc.lib.master() and np.any(source_energies < energy[0]):
-            print(' WARNING: Source point below energy grid')
-            sys.stdout.flush()
-        if openmc.lib.master() and np.any(source_energies > energy[-1]):
-            print(' WARNING: Source point above energy grid')
-            sys.stdout.flush()
-
-    def _count_bank_sites(self):
-        """Determines the number of fission bank sites in each cell of a given
-        mesh and energy group structure.
-        Returns
-        -------
-        bool
-            Wheter any source sites outside of CMFD mesh were found
-
-        """
-        # Initialize variables
-        m = openmc.lib.meshes[self._mesh_id]
-        bank = openmc.lib.source_bank()
-        energy = self._egrid
-        sites_outside = np.zeros(1, dtype=bool)
-        nxnynz = np.prod(self._indices[0:3])
-        ng = self._indices[3]
-
-        outside = np.zeros(1, dtype=bool)
-        self._sourcecounts = np.zeros((nxnynz, ng))
-        count = np.zeros(self._sourcecounts.shape)
-
-        # Get location and energy of each particle in source bank
-        source_xyz = openmc.lib.source_bank()['r']
-        source_energies = openmc.lib.source_bank()['E']
-
-        # Convert xyz location to mesh index and ravel index to scalar
-        mesh_locations = np.floor((source_xyz - m.lower_left) / m.width)
-        mesh_bins = mesh_locations[:,2] * m.dimension[1] * m.dimension[0] + \
-            mesh_locations[:,1] * m.dimension[0] + mesh_locations[:,0]
-
-        # Check if any source locations lie outside of defined CMFD mesh
-        if np.any(mesh_bins < 0) or np.any(mesh_bins >= np.prod(m.dimension)):
-            outside[0] = True
-
-        # Determine which energy bin each particle's energy belongs to
-        # Separate into cases bases on where source energies lies on egrid
-        energy_bins = np.zeros(len(source_energies), dtype=int)
-        idx = np.where(source_energies < energy[0])
-        energy_bins[idx] = 0
-        idx = np.where(source_energies > energy[-1])
-        energy_bins[idx] = ng - 1
-        idx = np.where((source_energies >= energy[0]) &
-                       (source_energies <= energy[-1]))
-        energy_bins[idx] = np.digitize(source_energies[idx], energy) - 1
-
-        # Determine all unique combinations of mesh bin and energy bin, and
-        # count number of particles that belong to these combinations
-        idx, counts = np.unique(np.array([mesh_bins, energy_bins]), axis=1,
-                                return_counts=True)
-
-        # Store counts to appropriate mesh-energy combination
-        count[idx[0].astype(int), idx[1].astype(int)] = counts
-
-        if have_mpi:
-            # Collect values of count from all processors
-            self._intracomm.Reduce(count, self._sourcecounts, MPI.SUM)
-            # Check if there were sites outside the mesh for any processor
-            self._intracomm.Reduce(outside, sites_outside, MPI.LOR)
-        # Deal with case if MPI not defined (only one proc)
-        else:
-            sites_outside = outside
-            self._sourcecounts = count
-
-        return sites_outside[0]
-
     def _build_loss_matrix(self, adjoint):
         # Extract spatial and energy indices and define matrix dimension
         ng = self._indices[3]
@@ -3045,3 +2902,7 @@ class CMFDRun:
 
             # Set all tallies to be active from beginning
             cmfd_tally.active = True
+
+        # Initialize CMFD mesh and energy grid in C++ for CMFD reweight
+        args = self._tally_ids[0], self._indices, self._norm
+        openmc.lib._dll.openmc_initialize_mesh_egrid(*args)
diff --git a/openmc/lib/core.py b/openmc/lib/core.py
index 3fe75dccf..7be94e03d 100644
--- a/openmc/lib/core.py
+++ b/openmc/lib/core.py
@@ -31,6 +31,8 @@ _array_1d_dble = np.ctypeslib.ndpointer(dtype=np.double, ndim=1,
 
 _dll.openmc_calculate_volumes.restype = c_int
 _dll.openmc_calculate_volumes.errcheck = _error_handler
+_dll.openmc_cmfd_reweight.argtypes = c_bool, _array_1d_dble
+_dll.openmc_cmfd_reweight.restype = None
 _dll.openmc_finalize.restype = c_int
 _dll.openmc_finalize.errcheck = _error_handler
 _dll.openmc_find_cell.argtypes = [POINTER(c_double*3), POINTER(c_int32),
@@ -45,8 +47,12 @@ _dll.openmc_init.errcheck = _error_handler
 _dll.openmc_get_keff.argtypes = [POINTER(c_double*2)]
 _dll.openmc_get_keff.restype = c_int
 _dll.openmc_get_keff.errcheck = _error_handler
+_dll.openmc_initialize_mesh_egrid.argtypes = [
+    c_int, _array_1d_int, c_double
+]
+_dll.openmc_initialize_mesh_egrid.restype = None
 _init_linsolver_argtypes = [_array_1d_int, c_int, _array_1d_int, c_int, c_int,
-                            c_double, _array_1d_int, _array_1d_int, c_bool]
+                            c_double, _array_1d_int, c_bool]
 _dll.openmc_initialize_linsolver.argtypes = _init_linsolver_argtypes
 _dll.openmc_initialize_linsolver.restype = None
 _dll.openmc_is_statepoint_batch.restype = c_bool
diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp
index 8406c9899..98ba0f6b4 100644
--- a/src/cmfd_solver.cpp
+++ b/src/cmfd_solver.cpp
@@ -8,9 +8,15 @@
 #endif
 #include "xtensor/xtensor.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/error.h"
 #include "openmc/constants.h"
 #include "openmc/capi.h"
+#include "openmc/mesh.h"
+#include "openmc/message_passing.h"
+#include "openmc/tallies/filter_energy.h"
+#include "openmc/tallies/filter_mesh.h"
+#include "openmc/tallies/tally.h"
 
 namespace openmc {
 
@@ -34,8 +40,200 @@ xt::xtensor<int, 2> indexmap;
 
 int use_all_threads;
 
+RegularMesh* mesh;
+
+std::vector<double> egrid;
+
+double norm;
+
 } // namespace cmfd
 
+//==============================================================================
+// GET_CMFD_ENERGY_BIN returns the energy bin for a source site energy
+//==============================================================================
+
+int get_cmfd_energy_bin(const double E)
+{
+  // Check if energy is out of grid bounds
+  if (E < cmfd::egrid[0]) {
+    // throw warning message
+    warning("Detected source point below energy grid");
+    return 0;
+  } else if (E >= cmfd::egrid[cmfd::ng]) {
+    // throw warning message
+    warning("Detected source point above energy grid");
+    return cmfd::ng - 1;
+  } else {
+    // Iterate through energy grid to find matching bin
+    for (int g = 0; g < cmfd::ng; g++) {
+      if (E >= cmfd::egrid[g] && E < cmfd::egrid[g+1]) {
+        return g;
+      }
+    }
+  }
+  // Return -1 by default
+  return -1;
+}
+
+//==============================================================================
+// COUNT_BANK_SITES bins fission sites according to CMFD mesh and energy
+//==============================================================================
+
+xt::xtensor<double, 1> count_bank_sites(xt::xtensor<int, 1>& bins, bool* outside)
+{
+  // Determine shape of array for counts
+  std::size_t cnt_size = cmfd::nx * cmfd::ny * cmfd::nz * cmfd::ng;
+  std::vector<std::size_t> cnt_shape = {cnt_size};
+
+  // Create array of zeros
+  xt::xarray<double> cnt {cnt_shape, 0.0};
+  bool outside_ = false;
+
+  auto bank_size = simulation::source_bank.size();
+  for (int i = 0; i < bank_size; i++) {
+    const auto& site = simulation::source_bank[i];
+
+    // determine scoring bin for CMFD mesh
+    int mesh_bin = cmfd::mesh->get_bin(site.r);
+
+    // if outside mesh, skip particle
+    if (mesh_bin < 0) {
+      outside_ = true;
+      continue;
+    }
+
+    // determine scoring bin for CMFD energy
+    int energy_bin = get_cmfd_energy_bin(site.E);
+
+    // add to appropriate bin
+    cnt(mesh_bin*cmfd::ng+energy_bin) += site.wgt;
+
+    // store bin index which is used again when updating weights
+    bins[i] = mesh_bin*cmfd::ng+energy_bin;
+  }
+
+  // Create copy of count data. Since ownership will be acquired by xtensor,
+  // std::allocator must be used to avoid Valgrind mismatched free() / delete
+  // warnings.
+  int total = cnt.size();
+  double* cnt_reduced = std::allocator<double>{}.allocate(total);
+
+#ifdef OPENMC_MPI
+  // collect values from all processors
+  MPI_Reduce(cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0,
+    mpi::intracomm);
+
+  // Check if there were sites outside the mesh for any processor
+  MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm);
+
+#else
+  std::copy(cnt.data(), cnt.data() + total, cnt_reduced);
+  *outside = outside_;
+#endif
+
+  // Adapt reduced values in array back into an xarray
+  auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), cnt_shape);
+  xt::xarray<double> counts = arr;
+
+  return counts;
+}
+
+//==============================================================================
+// OPENMC_CMFD_REWEIGHT performs reweighting of particles in source bank
+//==============================================================================
+
+extern "C"
+void openmc_cmfd_reweight(const bool feedback, const double* cmfd_src)
+{
+  // Get size of source bank and cmfd_src
+  auto bank_size = simulation::source_bank.size();
+  std::size_t src_size = cmfd::nx * cmfd::ny * cmfd::nz * cmfd::ng;
+
+  // count bank sites for CMFD mesh, store bins in bank_bins for reweighting
+  xt::xtensor<int, 1> bank_bins({bank_size}, 0);
+  bool sites_outside;
+  xt::xtensor<double, 1> sourcecounts = count_bank_sites(bank_bins,
+                                                         &sites_outside);
+
+  // Compute CMFD weightfactors
+  xt::xtensor<double, 1> weightfactors = xt::xtensor<double, 1>({src_size}, 1.);
+  if (mpi::master) {
+    if (sites_outside) {
+      fatal_error("Source sites outside of the CMFD mesh");
+    }
+
+    double norm = xt::sum(sourcecounts)()/cmfd::norm;
+    for (int i = 0; i < src_size; i++) {
+      if (sourcecounts[i] > 0 && cmfd_src[i] > 0) {
+        weightfactors[i] = cmfd_src[i] * norm / sourcecounts[i];
+      }
+    }
+  }
+
+  if (!feedback) return;
+
+#ifdef OPENMC_MPI
+  // Send weightfactors to all processors
+  MPI_Bcast(weightfactors.data(), src_size, MPI_DOUBLE, 0, mpi::intracomm);
+#endif
+
+  // Iterate through fission bank and update particle weights
+  for (int64_t i = 0; i < bank_size; i++) {
+    auto& site = simulation::source_bank[i];
+    site.wgt *= weightfactors(bank_bins(i));
+  }
+}
+
+//==============================================================================
+// OPENMC_INITIALIZE_MESH_EGRID sets the mesh and energy grid for CMFD reweight
+//==============================================================================
+
+extern "C"
+void openmc_initialize_mesh_egrid(const int meshtally_id, const int* cmfd_indices,
+                                  const double norm)
+{
+  // Make sure all CMFD memory is freed
+  free_memory_cmfd();
+
+  // Set CMFD indices
+  cmfd::nx = cmfd_indices[0];
+  cmfd::ny = cmfd_indices[1];
+  cmfd::nz = cmfd_indices[2];
+  cmfd::ng = cmfd_indices[3];
+
+  // Set CMFD reweight properties
+  cmfd::norm = norm;
+
+  // Find index corresponding to tally id
+  int32_t tally_index;
+  openmc_get_tally_index(meshtally_id, &tally_index);
+
+  // Get filters assocaited with tally
+  const auto& tally_filters = model::tallies[tally_index]->filters();
+
+  // Get mesh filter index
+  auto meshfilter_index = tally_filters[0];
+
+  // Store energy filter index if defined, otherwise set to -1
+  auto energy_index = (tally_filters.size() == 2) ? tally_filters[1] : -1;
+
+  // Get mesh index from mesh filter index
+  int32_t mesh_index;
+  openmc_mesh_filter_get_mesh(meshfilter_index, &mesh_index);
+
+  // Get mesh from mesh index
+  cmfd::mesh = get_regular_mesh(mesh_index);
+
+  // Get energy bins from energy index, otherwise use default
+  if (energy_index != -1) {
+    auto efilt_base = model::tally_filters[energy_index].get();
+    auto* efilt = dynamic_cast<EnergyFilter*>(efilt_base);
+    cmfd::egrid = efilt->bins();
+  } else {
+    cmfd::egrid = {0.0, INFTY};
+  }
+}
+
 //==============================================================================
 // MATRIX_TO_INDICES converts a matrix index to spatial and group
 // indices
@@ -309,12 +507,9 @@ int cmfd_linsolver_ng(const double* A_data, const double* b, double* x,
 extern "C"
 void openmc_initialize_linsolver(const int* indptr, int len_indptr,
                                  const int* indices, int n_elements, int dim,
-                                 double spectral, const int* cmfd_indices,
-                                 const int* map, bool use_all_threads)
+                                 double spectral, const int* map,
+                                 bool use_all_threads)
 {
-  // Make sure vectors are empty
-  free_memory_cmfd();
-
   // Store elements of indptr
   for (int i = 0; i < len_indptr; i++)
     cmfd::indptr.push_back(indptr[i]);
@@ -327,15 +522,8 @@ void openmc_initialize_linsolver(const int* indptr, int len_indptr,
   cmfd::dim = dim;
   cmfd::spectral = spectral;
 
-  // Set number of groups
-  cmfd::ng = cmfd_indices[3];
-
-  // Set problem dimensions and indexmap if 1 or 2 group problem
+  // Set indexmap if 1 or 2 group problem
   if (cmfd::ng == 1 || cmfd::ng == 2) {
-    cmfd::nx = cmfd_indices[0];
-    cmfd::ny = cmfd_indices[1];
-    cmfd::nz = cmfd_indices[2];
-
     // Resize indexmap and set its elements
     cmfd::indexmap.resize({static_cast<size_t>(dim), 3});
     set_indexmap(map);
@@ -366,10 +554,16 @@ int openmc_run_linsolver(const double* A_data, const double* b, double* x,
 
 void free_memory_cmfd()
 {
+  // Clear std::vectors
   cmfd::indptr.clear();
   cmfd::indices.clear();
-  // Resize indexmap to be an empty array
+  cmfd::egrid.clear();
+
+  // Resize xtensors to be empty
   cmfd::indexmap.resize({0});
+
+  // Set pointers to null
+  cmfd::mesh = nullptr;
 }
 
 } // namespace openmc
diff --git a/tests/regression_tests/cmfd_feed/results_true.dat b/tests/regression_tests/cmfd_feed/results_true.dat
index ea25c8856..85498af19 100644
--- a/tests/regression_tests/cmfd_feed/results_true.dat
+++ b/tests/regression_tests/cmfd_feed/results_true.dat
@@ -16,7 +16,7 @@ tally 1:
 3.347757E+01
 1.124264E+02
 2.931336E+01
-8.607238E+01
+8.607239E+01
 2.182947E+01
 4.789565E+01
 1.147668E+01
@@ -378,7 +378,7 @@ tally 5:
 5.253064E+00
 1.396224E+00
 2.996497E+01
-4.508839E+01
+4.508840E+01
 3.818076E+00
 7.509442E-01
 1.574994E+01
@@ -461,17 +461,17 @@ cmfd dominance ratio
 cmfd openmc source comparison
 6.959834E-03
 5.655657E-03
-3.886186E-03
+3.886185E-03
 4.035116E-03
 3.043277E-03
-5.455474E-03
-4.515310E-03
+5.455475E-03
+4.515311E-03
 2.439840E-03
 2.114032E-03
-2.673131E-03
+2.673132E-03
 2.431749E-03
 4.330928E-03
-3.404646E-03
+3.404647E-03
 3.680298E-03
 3.309620E-03
 3.705541E-03
diff --git a/tests/regression_tests/cmfd_feed_2g/results_true.dat b/tests/regression_tests/cmfd_feed_2g/results_true.dat
index d928b33e8..4b3dc3c96 100644
--- a/tests/regression_tests/cmfd_feed_2g/results_true.dat
+++ b/tests/regression_tests/cmfd_feed_2g/results_true.dat
@@ -80,7 +80,7 @@ tally 3:
 0.000000E+00
 0.000000E+00
 1.796382E-02
-3.190855E-05
+3.190854E-05
 4.343238E+00
 9.514040E-01
 3.504683E+00
@@ -90,7 +90,7 @@ tally 3:
 9.818871E+01
 4.822799E+02
 8.401346E-01
-3.664038E-02
+3.664037E-02
 6.130607E+01
 1.882837E+02
 0.000000E+00
@@ -212,7 +212,7 @@ tally 4:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-9.042232E+00
+9.042233E+00
 4.097018E+00
 3.711597E+01
 6.891834E+01
@@ -256,7 +256,7 @@ tally 4:
 4.161211E+00
 3.715168E+01
 6.904544E+01
-9.042232E+00
+9.042233E+00
 4.097018E+00
 3.711597E+01
 6.891834E+01
@@ -314,7 +314,7 @@ tally 5:
 1.007512E+02
 5.078162E+02
 1.403236E+01
-9.884346E+00
+9.884347E+00
 4.598486E+01
 1.058070E+02
 6.194677E+01
@@ -377,13 +377,13 @@ cmfd balance
 4.16856E-04
 6.38469E-04
 3.92822E-04
-3.78984E-04
+3.78983E-04
 2.68486E-04
 4.84991E-04
 1.08402E-03
-1.09177E-03
+1.09178E-03
 5.45977E-04
-4.45554E-04
+4.45555E-04
 4.01147E-04
 3.71025E-04
 3.57715E-04
@@ -408,21 +408,21 @@ cmfd dominance ratio
 5.996E-03
 cmfd openmc source comparison
 1.931386E-05
-2.839161E-05
-1.407963E-05
-6.718156E-06
-9.164199E-06
-1.382540E-05
+2.839162E-05
+1.407962E-05
+6.718148E-06
+9.164193E-06
+1.382539E-05
 2.871606E-05
 4.192300E-05
-5.327515E-05
+5.327516E-05
 6.566500E-05
-6.655540E-05
+6.655541E-05
 6.034977E-05
 6.004677E-05
-5.711622E-05
-6.271263E-05
-6.425362E-05
+5.711623E-05
+6.271264E-05
+6.425363E-05
 cmfd source
 2.510278E-01
 2.480592E-01
diff --git a/tests/regression_tests/cmfd_feed_expanding_window/results_true.dat b/tests/regression_tests/cmfd_feed_expanding_window/results_true.dat
index 897f9e9a8..8253d82c3 100644
--- a/tests/regression_tests/cmfd_feed_expanding_window/results_true.dat
+++ b/tests/regression_tests/cmfd_feed_expanding_window/results_true.dat
@@ -37,7 +37,7 @@ tally 2:
 3.708857E+00
 1.375562E+01
 2.698010E+00
-7.279259E+00
+7.279260E+00
 4.233192E+00
 1.791991E+01
 3.023990E+00
@@ -51,11 +51,11 @@ tally 2:
 2.497751E+00
 6.238760E+00
 2.974675E+00
-8.848691E+00
+8.848690E+00
 2.121163E+00
 4.499332E+00
 2.329248E+00
-5.425398E+00
+5.425397E+00
 1.658233E+00
 2.749736E+00
 1.158552E+00
@@ -109,7 +109,7 @@ tally 4:
 0.000000E+00
 0.000000E+00
 1.420707E-01
-2.018410E-02
+2.018409E-02
 2.633150E-01
 6.933479E-02
 0.000000E+00
@@ -131,7 +131,7 @@ tally 4:
 2.633150E-01
 6.933479E-02
 1.420707E-01
-2.018410E-02
+2.018409E-02
 2.628613E-01
 6.909608E-02
 3.590382E-01
@@ -229,7 +229,7 @@ tally 4:
 4.909661E-01
 2.410477E-01
 5.152690E-01
-2.655022E-01
+2.655021E-01
 4.788649E-01
 2.293116E-01
 0.000000E+00
@@ -251,7 +251,7 @@ tally 4:
 4.788649E-01
 2.293116E-01
 5.152690E-01
-2.655022E-01
+2.655021E-01
 4.266876E-01
 1.820623E-01
 3.401342E-01
@@ -277,7 +277,7 @@ tally 4:
 4.266876E-01
 1.820623E-01
 3.724186E-01
-1.386957E-01
+1.386956E-01
 2.560013E-01
 6.553669E-02
 0.000000E+00
@@ -299,7 +299,7 @@ tally 4:
 2.560013E-01
 6.553669E-02
 3.724186E-01
-1.386957E-01
+1.386956E-01
 2.892134E-01
 8.364439E-02
 1.556176E-01
@@ -441,13 +441,13 @@ cmfd dominance ratio
 cmfd openmc source comparison
 7.459013E-03
 5.012869E-03
-1.770224E-03
-5.242540E-03
-3.888027E-03
+1.770225E-03
+5.242539E-03
+3.888028E-03
 6.653433E-03
 8.839928E-03
-5.456904E-03
-5.668412E-03
+5.456903E-03
+5.668411E-03
 4.016377E-03
 4.179381E-03
 cmfd source
diff --git a/tests/regression_tests/cmfd_feed_ng/results_true.dat b/tests/regression_tests/cmfd_feed_ng/results_true.dat
index d9660f7bb..707411459 100644
--- a/tests/regression_tests/cmfd_feed_ng/results_true.dat
+++ b/tests/regression_tests/cmfd_feed_ng/results_true.dat
@@ -572,7 +572,7 @@ cmfd entropy
 1.999693E+00
 cmfd balance
 8.98944E-04
-4.45874E-04
+4.45875E-04
 2.18562E-04
 2.83412E-04
 3.60924E-04
@@ -599,13 +599,13 @@ cmfd openmc source comparison
 2.585546E-05
 2.217204E-05
 2.121913E-05
-9.253928E-06
-2.940947E-05
-3.083596E-05
-2.547491E-05
-2.602691E-05
-2.561169E-05
-2.816825E-05
+9.253920E-06
+2.940945E-05
+3.083595E-05
+2.547490E-05
+2.602689E-05
+2.561168E-05
+2.816824E-05
 cmfd source
 2.417661E-01
 2.547768E-01
diff --git a/tests/regression_tests/cmfd_feed_ref_d/results_true.dat b/tests/regression_tests/cmfd_feed_ref_d/results_true.dat
index b1f14ac58..2f0d4713c 100644
--- a/tests/regression_tests/cmfd_feed_ref_d/results_true.dat
+++ b/tests/regression_tests/cmfd_feed_ref_d/results_true.dat
@@ -1,7 +1,7 @@
 k-combined:
 1.165408E+00 1.127320E-02
 tally 1:
-1.141008E+01
+1.141009E+01
 1.305436E+01
 2.074878E+01
 4.311543E+01
@@ -29,7 +29,7 @@ tally 2:
 1.991651E+00
 3.966673E+00
 1.411802E+00
-1.993185E+00
+1.993186E+00
 2.978294E+00
 8.870233E+00
 2.130084E+00
@@ -63,7 +63,7 @@ tally 2:
 8.282626E-01
 6.860190E-01
 tally 3:
-7.459524E-01
+7.459525E-01
 5.564451E-01
 4.877161E-02
 2.378670E-03
@@ -345,10 +345,10 @@ tally 4:
 0.000000E+00
 0.000000E+00
 tally 5:
-7.459524E-01
+7.459525E-01
 5.564451E-01
 1.149859E-01
-1.322176E-02
+1.322177E-02
 1.356729E+00
 1.840713E+00
 1.801152E-01
@@ -362,7 +362,7 @@ tally 5:
 3.360312E-01
 1.129170E-01
 2.909427E+00
-8.464764E+00
+8.464765E+00
 4.095398E-01
 1.677229E-01
 2.805095E+00
@@ -442,13 +442,13 @@ cmfd openmc source comparison
 7.433111E-03
 5.006211E-03
 1.766072E-03
-5.184425E-03
-3.864321E-03
-6.617152E-03
+5.184426E-03
+3.864323E-03
+6.617153E-03
 8.841210E-03
-5.454747E-03
-5.652690E-03
-4.006767E-03
+5.454745E-03
+5.652688E-03
+4.006766E-03
 4.167617E-03
 cmfd source
 4.116746E-02
diff --git a/tests/regression_tests/cmfd_feed_rolling_window/results_true.dat b/tests/regression_tests/cmfd_feed_rolling_window/results_true.dat
index 0974fcb8b..5d3ea4ec9 100644
--- a/tests/regression_tests/cmfd_feed_rolling_window/results_true.dat
+++ b/tests/regression_tests/cmfd_feed_rolling_window/results_true.dat
@@ -31,7 +31,7 @@ tally 2:
 1.324803E+00
 1.755104E+00
 2.585263E+00
-6.683587E+00
+6.683586E+00
 1.840969E+00
 3.389167E+00
 3.737263E+00
@@ -41,7 +41,7 @@ tally 2:
 3.989056E+00
 1.591257E+01
 2.843382E+00
-8.084820E+00
+8.084819E+00
 3.951710E+00
 1.561601E+01
 2.838216E+00
@@ -54,7 +54,7 @@ tally 2:
 1.050543E+01
 2.315684E+00
 5.362392E+00
-2.459493E+00
+2.459494E+00
 6.049108E+00
 1.734344E+00
 3.007949E+00
@@ -72,7 +72,7 @@ tally 3:
 7.284456E-02
 5.306329E-03
 1.779821E+00
-3.167762E+00
+3.167761E+00
 1.110012E-01
 1.232127E-02
 2.586020E+00
@@ -92,7 +92,7 @@ tally 3:
 1.503142E-01
 2.259435E-02
 2.233122E+00
-4.986832E+00
+4.986833E+00
 1.456891E-01
 2.122532E-02
 1.682929E+00
@@ -301,7 +301,7 @@ tally 4:
 4.201532E-01
 1.765287E-01
 2.898817E-01
-8.403142E-02
+8.403143E-02
 1.509185E-01
 2.277639E-02
 0.000000E+00
@@ -323,7 +323,7 @@ tally 4:
 1.509185E-01
 2.277639E-02
 2.898817E-01
-8.403142E-02
+8.403143E-02
 1.482175E-01
 2.196844E-02
 0.000000E+00
@@ -354,7 +354,7 @@ tally 5:
 1.458271E-01
 2.126554E-02
 1.779821E+00
-3.167762E+00
+3.167761E+00
 2.705917E-01
 7.321987E-02
 2.584069E+00
@@ -374,9 +374,9 @@ tally 5:
 3.296244E-01
 1.086523E-01
 2.233122E+00
-4.986832E+00
+4.986833E+00
 3.106254E-01
-9.648813E-02
+9.648814E-02
 1.682929E+00
 2.832250E+00
 2.709632E-01
@@ -440,16 +440,16 @@ cmfd dominance ratio
 5.321E-01
 cmfd openmc source comparison
 7.693485E-03
-4.158805E-03
+4.158806E-03
 2.962505E-03
-4.415044E-03
-2.304495E-03
-7.921580E-03
-8.609203E-03
-8.945198E-03
-8.054204E-03
+4.415043E-03
+2.304496E-03
+7.921579E-03
+8.609204E-03
+8.945200E-03
+8.054206E-03
 1.164189E-02
-5.945645E-03
+5.945644E-03
 cmfd source
 4.077779E-02
 6.659143E-02
diff --git a/tests/regression_tests/cmfd_restart/results_true.dat b/tests/regression_tests/cmfd_restart/results_true.dat
index ea25c8856..85498af19 100644
--- a/tests/regression_tests/cmfd_restart/results_true.dat
+++ b/tests/regression_tests/cmfd_restart/results_true.dat
@@ -16,7 +16,7 @@ tally 1:
 3.347757E+01
 1.124264E+02
 2.931336E+01
-8.607238E+01
+8.607239E+01
 2.182947E+01
 4.789565E+01
 1.147668E+01
@@ -378,7 +378,7 @@ tally 5:
 5.253064E+00
 1.396224E+00
 2.996497E+01
-4.508839E+01
+4.508840E+01
 3.818076E+00
 7.509442E-01
 1.574994E+01
@@ -461,17 +461,17 @@ cmfd dominance ratio
 cmfd openmc source comparison
 6.959834E-03
 5.655657E-03
-3.886186E-03
+3.886185E-03
 4.035116E-03
 3.043277E-03
-5.455474E-03
-4.515310E-03
+5.455475E-03
+4.515311E-03
 2.439840E-03
 2.114032E-03
-2.673131E-03
+2.673132E-03
 2.431749E-03
 4.330928E-03
-3.404646E-03
+3.404647E-03
 3.680298E-03
 3.309620E-03
 3.705541E-03