Move solver to C++

CMakeLists.txt

@@ -384,6 +384,7 @@ add_library(libopenmc SHARED
   src/dagmc.cpp
   src/cell.cpp
   src/cmfd_execute.cpp
+  src/cmfd_solver.cpp
   src/distribution.cpp
   src/distribution_angle.cpp
   src/distribution_energy.cpp

include/openmc/cmfd_solver.h

@@ -0,0 +1,113 @@
+#ifndef OPENMC_CMFD_SOLVER_H
+#define OPENMC_CMFD_SOLVER_H
+
+#include <cmath>
+
+#include "xtensor/xtensor.hpp"
+
+namespace openmc {
+
+//===============================================================================
+// Global variables
+//===============================================================================
+
+// CSR format index pointer array of loss matrix
+extern std::vector<int> indptr;
+
+// CSR format index array of loss matrix
+extern std::vector<int> indices;
+
+// Dimension n of nxn CMFD loss matrix
+extern int dim;
+
+// Spectral radius of CMFD matrices and tolerances
+extern double spectral;
+
+// Maximum dimension in x, y, and z directions
+extern int nx;
+extern int ny;
+extern int nz;
+
+// Number of energy groups
+extern int ng;
+
+// Indexmap storing all x, y, z positions of accelerated regions
+extern xt::xtensor<int, 2> indexmap;
+
+//===============================================================================
+// Non-member functions
+//===============================================================================
+
+//! returns the index in CSR index array corresponding to the diagonal element
+//! of a specified row
+//! \param[in] row of interest
+//! \return index in CSR index array corresponding to diagonal element
+int get_diagonal_index(int row);
+
+//! sets the elements of indexmap based on input coremap
+//! \param[in] user-defined coremap
+void set_indexmap(int* coremap);
+
+//! solves a one group CMFD linear system
+//! \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
+//! \return number of inner iterations required to reach convergence
+int cmfd_linsolver_1g(double* A_data, double* b, double* x, double tol);
+
+//! solves a two group CMFD linear system
+//! \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
+//! \return number of inner iterations required to reach convergence
+int cmfd_linsolver_2g(double* A_data, double* b, double* x, double tol);
+
+//! solves a general CMFD linear system
+//! \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
+//! \return number of inner iterations required to reach convergence
+int cmfd_linsolver_ng(double* A_data, double* b, double* x, double tol);
+
+//! converts a matrix index to spatial and group indices
+//! \param[in] iteration counter over row
+//! \param[out] iteration counter for groups
+//! \param[out] iteration counter for x
+//! \param[out] iteration counter for y
+//! \param[out] iteration counter for z
+void matrix_to_indices(int irow, int& g, int& i, int& j, int& k);
+
+//===============================================================================
+// External functions
+//===============================================================================
+
+//! sets the fixed variables that are used for the 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
+extern "C" void openmc_initialize_linsolver(int* indptr, int len_indptr,
+                                            int* indices, int n_elements,
+                                            int dim, double spectral,
+                                            int* cmfd_indices, int* map);
+
+//! 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
+//! \return number of inner iterations required to reach convergence
+extern "C" int openmc_run_linsolver(double* A_data, double* b, double* x,
+                                    double tol);
+
+
+} // namespace openmc
+#endif // OPENMC_CMFD_SOLVER_H

include/openmc/constants.h

@@ -434,6 +434,12 @@ constexpr int RUN_MODE_PLOTTING    {3};
 constexpr int RUN_MODE_PARTICLE    {4};
 constexpr int RUN_MODE_VOLUME      {5};
 
+// ============================================================================
+// CMFD CONSTANTS
+
+// For non-accelerated regions on coarse mesh overlay
+constexpr int CMFD_NOACCEL {99999};
+
 } // namespace openmc
 
 #endif // OPENMC_CONSTANTS_H

openmc/cmfd.py

@@ -18,8 +18,10 @@ import numpy as np
 # Line below is added to suppress warnings when using numpy.divide to
 # divide by numpy arrays that contain entries with zero
 np.seterr(divide='ignore', invalid='ignore')
+import numpy.ctypeslib as npct
 from scipy import sparse
 import time
+from ctypes import c_int, c_double
 # See if mpi4py module can be imported, define have_mpi global variable
 try:
     from mpi4py import MPI
@@ -33,6 +35,20 @@ from openmc.checkvalue import (check_type, check_length, check_value,
                                check_greater_than, check_less_than)
 from openmc.exceptions import OpenMCError
 
+
+# Define input type for numpy arrays that will be passed into C++ functions
+# Must be an int or double array, with single dimension that is contiguous
+array_1d_int = npct.ndpointer(dtype=np.int32, ndim=1, flags='CONTIGUOUS')
+array_1d_dble = npct.ndpointer(dtype=np.double, ndim=1, flags='CONTIGUOUS')
+
+# Setup the return types and argument types for C++ functions
+openmc.capi._dll.openmc_initialize_linsolver.restype = None
+openmc.capi._dll.openmc_initialize_linsolver.argtypes = [array_1d_int, c_int,
+    array_1d_int, c_int, c_int, c_double, array_1d_int, array_1d_int]
+openmc.capi._dll.openmc_run_linsolver.restype = c_int
+openmc.capi._dll.openmc_run_linsolver.argtypes = [array_1d_dble, array_1d_dble,
+    array_1d_dble, c_double]
+
 """
 --------------
 CMFD CONSTANTS
@@ -758,6 +774,9 @@ class CMFDRun(object):
     k_cmfd : list of floats
         List of CMFD k-effectives, stored for each generation that CMFD is
         invoked
+    spectral : float
+        Optional spectral radius that can be used to accelerate the convergence
+        of Gauss-Seidel iterations during CMFD power iteration.
     resnb : numpy.ndarray
         Residual from solving neutron balance equations
     time_cmfd : float
@@ -798,9 +817,11 @@ class CMFDRun(object):
         self._cmfd_stol = 1.e-8
         self._cmfd_reset = []
         self._cmfd_write_matrices = False
+        self._cmfd_spectral = 0.0
+        self._gauss_seidel_tolerance = [1.e-10, 1.e-5]
 
         # External variables used during runtime but users cannot control
-        self._indices = np.zeros(4, dtype=int)
+        self._indices = np.zeros(4, dtype=np.int32)
         self._egrid = None
         self._albedo = None
         self._coremap = None
@@ -852,12 +873,12 @@ class CMFDRun(object):
         self._notlast_y_accel = None
         self._notfirst_z_accel = None
         self._notlast_z_accel = None
-        self._adj_reflector_left = None
-        self._adj_reflector_right = None
-        self._adj_reflector_back = None
-        self._adj_reflector_front = None
-        self._adj_reflector_bottom = None
-        self._adj_reflector_top = None
+        self._is_adj_ref_left = None
+        self._is_adj_ref_right = None
+        self._is_adj_ref_back = None
+        self._is_adj_ref_front = None
+        self._is_adj_ref_bottom = None
+        self._is_adj_ref_top = None
         self._accel_idxs = None
         self._accel_neig_left_idxs = None
         self._accel_neig_right_idxs = None
@@ -867,6 +888,8 @@ class CMFDRun(object):
         self._accel_neig_top_idxs = None
         self._loss_row = None
         self._loss_col = None
+        self._prod_row = None
+        self._prod_col = None
         self._intracomm = None
 
 
@@ -922,6 +945,10 @@ class CMFDRun(object):
     def cmfd_stol(self):
         return self._cmfd_stol
 
+    @property
+    def cmfd_spectral(self):
+        return self._cmfd_spectral
+
     @property
     def cmfd_reset(self):
         return self._cmfd_reset
@@ -930,6 +957,10 @@ class CMFDRun(object):
     def cmfd_write_matrices(self):
         return self._cmfd_write_matrices
 
+    @property
+    def gauss_seidel_tolerance(self):
+        return self._gauss_seidel_tolerance
+
     @cmfd_begin.setter
     def cmfd_begin(self, cmfd_begin):
         check_type('CMFD begin batch', cmfd_begin, Integral)
@@ -1039,6 +1070,11 @@ class CMFDRun(object):
         check_type('CMFD fission source tolerance', cmfd_stol, Real)
         self._cmfd_stol = cmfd_stol
 
+    @cmfd_spectral.setter
+    def cmfd_spectral(self, spectral):
+        check_type('CMFD spectral radius', spectral, Real)
+        self._cmfd_spectral = spectral
+
     @cmfd_reset.setter
     def cmfd_reset(self, cmfd_reset):
         check_type('tally reset batches', cmfd_reset, Iterable, Integral)
@@ -1049,7 +1085,14 @@ class CMFDRun(object):
         check_type('CMFD write matrices', cmfd_write_matrices, bool)
         self._cmfd_write_matrices = cmfd_write_matrices
 
-    def run(self, omp_num_threads=None, intracomm=None, vectorized=True):
+    @gauss_seidel_tolerance.setter
+    def gauss_seidel_tolerance(self, gauss_seidel_tolerance):
+        check_type('CMFD Gauss-Seidel tolerance', gauss_seidel_tolerance,
+                   Iterable, Real)
+        check_length('Gauss-Seidel tolerance', gauss_seidel_tolerance, 2)
+        self._gauss_seidel_tolerance = gauss_seidel_tolerance
+
+    def run(self, omp_num_threads=None, intracomm=None, vectorized=True, cpp_solver=False):
         """Public method to run OpenMC with CMFD
 
         This method is called by user to run CMFD once instance variables of
@@ -1071,6 +1114,8 @@ class CMFDRun(object):
         elif intracomm is None and have_mpi:
             self._intracomm = MPI.COMM_WORLD
 
+        self._cpp_solver = cpp_solver
+
         # Check number of OpenMP threads is valid input and initialize C API
         if omp_num_threads is not None:
             check_type('OpenMP num threads', omp_num_threads, Integral)
@@ -1090,6 +1135,10 @@ class CMFDRun(object):
         # Compute and store row and column indices used to build CMFD matrices
         self._precompute_matrix_indices()
 
+        # Initialize all variables used for linear solver in C++
+        if self._cpp_solver:
+            self._initialize_linsolver()
+
         # Initialize simulation
         openmc.capi.simulation_init()
 
@@ -1120,6 +1169,23 @@ class CMFDRun(object):
         # Finalize and free memory
         openmc.capi.finalize()
 
+    def _initialize_linsolver(self):
+        # Determine number of rows in CMFD matrix
+        ng = self._indices[3]
+        n = self._mat_dim*ng
+
+        # Create temp loss matrix to pass row/col indices to C++ linear solver
+        temp_data = np.ones(len(self._loss_row))
+        temp_loss = sparse.csr_matrix((temp_data, (self._loss_row, self._loss_col)),
+                                      shape=(n, n))
+
+        # Pass coremap as 1-d array of 32-bit integers
+        coremap = np.swapaxes(self._coremap, 0, 2).flatten().astype(np.int32)
+
+        return openmc.capi._dll.openmc_initialize_linsolver(temp_loss.indptr,
+            len(temp_loss.indptr), temp_loss.indices, len(temp_loss.indices), n,
+            self._cmfd_spectral, self._indices, coremap)
+
     def _write_cmfd_output(self):
         """Write CMFD output to buffer at the end of each batch"""
         # Display CMFD k-effective
@@ -1935,6 +2001,12 @@ class CMFDRun(object):
         # Get problem size
         n = loss.shape[0]
 
+        # Set up tolerances for C++ solver
+        if self._cpp_solver:
+            atoli = self._gauss_seidel_tolerance[0]
+            rtoli = self._gauss_seidel_tolerance[1]
+            toli = rtoli * 100
+
         # Set up flux vectors, intital guess set to 1
         phi_n = np.ones((n,))
         phi_o = np.ones((n,))
@@ -1942,7 +2014,6 @@ class CMFDRun(object):
         # Set up source vectors
         s_n = np.zeros((n,))
         s_o = np.zeros((n,))
-        serr_v = np.zeros((n,))
 
         # Set initial guess
         k_n = openmc.capi.keff_temp()[0]
@@ -1975,8 +2046,13 @@ class CMFDRun(object):
             # Normalize source vector
             s_o /= k_lo
 
-            # Compute new flux vector with scipy sparse solver
-            phi_n = sparse.linalg.spsolve(loss, s_o)
+            # Compute new flux with either C++ solver or scipy sparse solver
+            if self._cpp_solver:
+                innerits = openmc.capi._dll.openmc_run_linsolver(loss.data,
+                    s_o, phi_n, toli)
+            else:
+                phi_n = sparse.linalg.spsolve(loss, s_o)
+                innerits = 0
 
             # Compute new source vector
             s_n = prod.dot(phi_n)
@@ -1991,7 +2067,8 @@ class CMFDRun(object):
             s_o *= k_lo
 
             # Check convergence
-            iconv, norm_n = self._check_convergence(s_n, s_o, k_n, k_o, i+1)
+            iconv, norm_n = self._check_convergence(s_n, s_o, k_n, k_o, i+1,
+                                                    innerits=innerits)
 
             # If converged, calculate dominance ratio and break from loop
             if iconv:
@@ -2004,7 +2081,11 @@ class CMFDRun(object):
             k_lo = k_ln
             norm_o = norm_n
 
-    def _check_convergence(self, s_n, s_o, k_n, k_o, iter):
+            # Update tolerance for inner iterations
+            if self._cpp_solver:
+                toli = max(atoli, rtoli*norm_n)
+
+    def _check_convergence(self, s_n, s_o, k_n, k_o, iter, innerits=0):
         """Checks the convergence of the CMFD problem
 
         Parameters
@@ -2044,14 +2125,19 @@ class CMFDRun(object):
             str2 = 'k-eff: {:0.8f}'.format(k_n)
             str3 = 'k-error:  {0:.5E}'.format(kerr)
             str4 = 'src-error:  {0:.5E}'.format(serr)
-            print('{0:8s}{1:20s}{2:25s}{3:s}'.format(str1, str2, str3, str4))
+            if innerits:
+                str5 = '  {:d}'.format(innerits)
+                print('{0:8s}{1:20s}{2:25s}{3:s}{4:s}'.format(str1, str2, str3,
+                                                              str4, str5))
+            else:
+                print('{0:8s}{1:20s}{2:25s}{3:s}'.format(str1, str2, str3, str4))
             sys.stdout.flush()
 
         return iconv, serr
 
     def _set_coremap(self):
         """Sets the core mapping information. All regions marked with zero
-        are set to CMFD_NO_ACCEL, while all regions marked with 1 are set to a
+        are set to CMFD_NOACCEL, while all regions marked with 1 are set to a
         unique index that maps each fuel region to a row number when building
         CMFD matrices
 
@@ -2485,8 +2571,8 @@ class CMFDRun(object):
                 mode='constant', constant_values=_CMFD_NOACCEL)[:,:,1:]
 
         # Create empty row and column vectors to store for loss matrix
-        row = np.array([], dtype=int)
-        col = np.array([], dtype=int)
+        row = np.array([])
+        col = np.array([])
 
         # Store all indices used to populate production and loss matrix
         self._accel_idxs = np.where(self._coremap != _CMFD_NOACCEL)
@@ -3552,4 +3638,4 @@ class CMFDRun(object):
         if current[2*l] < 1.0e-10:
             return 1.0
         else:
-          return current[2*l+1]/current[2*l]
+          return current[2*l+1]/current[2*l]

src/cmfd_solver.cpp

@@ -0,0 +1,339 @@
+//TODO remove
+#include <iostream>
+#include <vector>
+
+#include "openmc/cmfd_solver.h"
+#include "openmc/error.h"
+#include "openmc/constants.h"
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+// TODO check which variables actually necessary
+
+std::vector<int> indptr;
+
+std::vector<int> indices;
+
+int dim;
+
+double spectral;
+
+int nx, ny, nz, ng;
+
+xt::xtensor<int, 2> indexmap;
+
+//==============================================================================
+// GET_DIAGONAL_INDEX returns the index in CSR index array corresponding to
+// the diagonal element of a specified row
+//==============================================================================
+
+int get_diagonal_index(int row) {
+  for (int j = indptr[row]; j < indptr[row+1]; j++) {
+    if (indices[j] == row)
+      return j;
+  }
+
+  // Return -1 if not found
+  return -1;
+}
+
+//==============================================================================
+// SET_INDEXMAP sets the elements of indexmap based on input coremap
+//==============================================================================
+
+void set_indexmap(int* coremap) {
+  for (int z = 0; z < nz; z++) {
+    for (int y = 0; y < ny; y++) {
+      for (int x = 0; x < nx; x++) {
+        if (coremap[(z*ny*nx) + (y*nx) + x] != CMFD_NOACCEL) {
+          int counter = coremap[(z*ny*nx) + (y*nx) + x];
+          indexmap(counter, 0) = x;
+          indexmap(counter, 1) = y;
+          indexmap(counter, 2) = z;
+        }
+      }
+    }
+  }
+}
+
+//==============================================================================
+// CMFD_LINSOLVER_1G solves a one group CMFD linear system
+//==============================================================================
+
+int cmfd_linsolver_1g(double* A_data, double* b, double* x, double tol) {
+  // Set overrelaxation parameter
+  double w = 1.0;
+
+  // Perform Gauss-Seidel iterations
+  for (int igs = 1; igs <= 10000; igs++) {
+    double tmpx[dim];
+    double err = 0.0;
+
+    // Copy over x vector
+    std::copy(x, x+dim, tmpx);
+
+    // Perform red/black Gauss-Seidel iterations
+    for (int irb = 0; irb < 2; irb++) {
+
+      // Loop around matrix rows
+      for (int irow = 0; irow < dim; irow++) {
+        int g, i, j, k;
+        matrix_to_indices(irow, g, i, j, k);
+
+        // Filter out black cells
+        if ((i+j+k) % 2 != irb) continue;
+
+        // Get index of diagonal for current row
+        int didx = get_diagonal_index(irow);
+
+        // Perform temporary sums, first do left of diag, then right of diag
+        double tmp1 = 0.0;
+        for (int icol = indptr[irow]; icol < didx; icol++)
+          tmp1 += A_data[icol] * x[indices[icol]];
+        for (int icol = didx + 1; icol < indptr[irow + 1]; icol++)
+          tmp1 += A_data[icol] * x[indices[icol]];
+
+        // Solve for new x
+        double x1 = (b[irow] - tmp1) / A_data[didx];
+
+        // Perform overrelaxation
+        x[irow] = (1.0 - w) * x[irow] + w * x1;
+
+        // Compute residual and update error
+        double res = (tmpx[irow] - x[irow]) / tmpx[irow];
+        err += res * res;
+      }
+    }
+
+    // Check convergence
+    err = std::sqrt(err / dim);
+    if (err < tol)
+      return igs;
+
+    // Calculate new overrelaxation parameter
+    w = 1.0/(1.0 - 0.25 * spectral * w);
+  }
+
+  // Throw error, as max iterations met
+  fatal_error("Maximum Gauss-Seidel iterations encountered.");
+
+  // Return -1 by default, although error thrown before reaching this point
+  return -1;
+}
+
+//==============================================================================
+// CMFD_LINSOLVER_2G solves a two group CMFD linear system
+//==============================================================================
+
+int cmfd_linsolver_2g(double* A_data, double* b, double* x, double tol) {
+  // Set overrelaxation parameter
+  double w = 1.0;
+
+  // Perform Gauss-Seidel iterations
+  for (int igs = 1; igs <= 10000; igs++) {
+    double tmpx[dim];
+    double err = 0.0;
+
+    // Copy over x vector
+    std::copy(x, x+dim, tmpx);
+
+    // Perform red/black Gauss-Seidel iterations
+    for (int irb = 0; irb < 2; irb++) {
+
+      // Loop around matrix rows
+      for (int irow = 0; irow < dim; irow+=2) {
+        int g, i, j, k;
+        matrix_to_indices(irow, g, i, j, k);
+
+        // Filter out black cells
+        if ((i+j+k) % 2 != irb) continue;
+
+        // Get index of diagonals for current row and next row
+        int d1idx = get_diagonal_index(irow);
+        int d2idx = get_diagonal_index(irow+1);
+
+        // Get block diagonal
+        double m11 = A_data[d1idx];     // group 1 diagonal
+        double m12 = A_data[d1idx + 1]; // group 1 right of diagonal (sorted by col)
+        double m21 = A_data[d2idx - 1]; // group 2 left of diagonal (sorted by col)
+        double m22 = A_data[d2idx];     // group 2 diagonal
+
+        // Analytically invert the diagonal
+        double dm = m11*m22 - m12*m21;
+        double d11 = m22/dm;
+        double d12 = -m12/dm;
+        double d21 = -m21/dm;
+        double d22 = m11/dm;
+
+        // Perform temporary sums, first do left of diag, then right of diag
+        double tmp1 = 0.0;
+        double tmp2 = 0.0;
+        for (int icol = indptr[irow]; icol < d1idx; icol++)
+          tmp1 += A_data[icol] * x[indices[icol]];
+        for (int icol = indptr[irow+1]; icol < d2idx-1; icol++)
+          tmp2 += A_data[icol] * x[indices[icol]];
+        for (int icol = d1idx + 2; icol < indptr[irow + 1]; icol++)
+          tmp1 += A_data[icol] * x[indices[icol]];
+        for (int icol = d2idx + 1; icol < indptr[irow + 2]; icol++)
+          tmp1 += A_data[icol] * x[indices[icol]];
+
+        // Adjust with RHS vector
+        tmp1 = b[irow] - tmp1;
+        tmp2 = b[irow + 1] - tmp2;
+
+        // Solve for new x
+        double x1 = d11*tmp1 + d12*tmp2;
+        double x2 = d21*tmp1 + d22*tmp2;
+
+        // Perform overrelaxation
+        x[irow] = (1.0 - w) * x[irow] + w * x1;
+        x[irow + 1] = (1.0 - w) * x[irow + 1] + w * x2;
+
+        // Compute residual and update error
+        double res = (tmpx[irow] - x[irow]) / tmpx[irow];
+        err += res * res;
+      }
+    }
+
+    // Check convergence
+    err = std::sqrt(err / dim);
+    if (err < tol)
+      return igs;
+
+    // Calculate new overrelaxation parameter
+    w = 1.0/(1.0 - 0.25 * spectral * w);
+  }
+
+  // Throw error, as max iterations met
+  fatal_error("Maximum Gauss-Seidel iterations encountered.");
+
+  // Return -1 by default, although error thrown before reaching this point
+  return -1;
+}
+
+//==============================================================================
+// CMFD_LINSOLVER_NG solves a general CMFD linear system
+//==============================================================================
+
+int cmfd_linsolver_ng(double* A_data, double* b, double* x, double tol) {
+  // Set overrelaxation parameter
+  double w = 1.0;
+
+  // Perform Gauss-Seidel iterations
+  for (int igs = 1; igs <= 10000; igs++) {
+    double tmpx[dim];
+    double err = 0.0;
+
+    // Copy over x vector
+    std::copy(x, x+dim, tmpx);
+
+    // Loop around matrix rows
+    for (int irow = 0; irow < dim; irow++) {
+      // Get index of diagonal for current row
+      int didx = get_diagonal_index(irow);
+
+      // Perform temporary sums, first do left of diag, then right of diag
+      double tmp1 = 0.0;
+      for (int icol = indptr[irow]; icol < didx; icol++)
+        tmp1 += A_data[icol] * x[indices[icol]];
+      for (int icol = didx + 1; icol < indptr[irow + 1]; icol++)
+        tmp1 += A_data[icol] * x[indices[icol]];
+
+      // Solve for new x
+      double x1 = (b[irow] - tmp1) / A_data[didx];
+
+      // Perform overrelaxation
+      x[irow] = (1.0 - w) * x[irow] + w * x1;
+
+      // Compute residual and update error
+      double res = (tmpx[irow] - x[irow]) / tmpx[irow];
+      err += res * res;
+    }
+
+    // Check convergence
+    err = std::sqrt(err / dim);
+    if (err < tol)
+      return igs;
+
+    // Calculate new overrelaxation parameter
+    w = 1.0/(1.0 - 0.25 * spectral * w);
+  }
+
+  // Throw error, as max iterations met
+  fatal_error("Maximum Gauss-Seidel iterations encountered.");
+
+  // Return -1 by default, although error thrown before reaching this point
+  return -1;
+}
+
+//==============================================================================
+// MATRIX_TO_INDICES converts a matrix index to spatial and group
+// indices
+//==============================================================================
+
+void matrix_to_indices(int irow, int& g, int& i, int& j, int& k) {
+  g = irow % ng;
+  i = indexmap(irow/ng, 0);
+  j = indexmap(irow/ng, 1);
+  k = indexmap(irow/ng, 2);
+}
+
+//==============================================================================
+// OPENMC_INITIALIZE_LINSOLVER sets the fixed variables that are used for the
+// linear solver
+//==============================================================================
+
+extern "C"
+void openmc_initialize_linsolver(int* indptr, int len_indptr, int* indices,
+                                 int n_elements, int dim, double spectral,
+                                 int* cmfd_indices, int* map) {
+  // Store elements of indptr
+  for (int i = 0; i < len_indptr; i++)
+    openmc::indptr.push_back(indptr[i]);
+
+  // Store elements of indices
+  for (int i = 0; i < n_elements; i++)
+    openmc::indices.push_back(indices[i]);
+
+  // Set dimenion of CMFD problem and specral radius
+  openmc::dim = dim;
+  openmc::spectral = spectral;
+
+  // Set number of groups
+  openmc::ng = cmfd_indices[3];
+
+  // Set problem dimensions and indexmap if 1 or 2 group problem
+  if (openmc::ng == 1 || openmc::ng == 2) {
+    openmc::nx = cmfd_indices[0];
+    openmc::ny = cmfd_indices[1];
+    openmc::nz = cmfd_indices[2];
+
+    // Resize indexmap and set its elements
+    openmc::indexmap.resize({static_cast<size_t>(dim), 3});
+    set_indexmap(map);
+  }
+}
+
+//==============================================================================
+// OPENMC_RUN_LINSOLVER runs a Gauss Seidel linear solver to solve CMFD matrix
+// equations
+//==============================================================================
+
+extern "C"
+int openmc_run_linsolver(double* A_data, double* b, double* x, double tol) {
+  switch (ng) {
+  case 1:
+    return cmfd_linsolver_1g(A_data, b, x, tol);
+  case 2:
+    return cmfd_linsolver_2g(A_data, b, x, tol);
+  default:
+    return cmfd_linsolver_ng(A_data, b, x, tol);
+  }
+}
+
+
+} // namespace openmc