From e2c2a0322095df20355e27c99cd3db46ec6e80e2 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Sat, 10 Nov 2018 12:11:42 -0500 Subject: [PATCH] Move solver to C++ --- CMakeLists.txt | 1 + include/openmc/cmfd_solver.h | 113 ++++++++++++ include/openmc/constants.h | 6 + openmc/cmfd.py | 122 +++++++++++-- src/cmfd_solver.cpp | 339 +++++++++++++++++++++++++++++++++++ 5 files changed, 563 insertions(+), 18 deletions(-) create mode 100644 include/openmc/cmfd_solver.h create mode 100644 src/cmfd_solver.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index fded5876c..b5d1946d9 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/openmc/cmfd_solver.h b/include/openmc/cmfd_solver.h new file mode 100644 index 000000000..47400b2a8 --- /dev/null +++ b/include/openmc/cmfd_solver.h @@ -0,0 +1,113 @@ +#ifndef OPENMC_CMFD_SOLVER_H +#define OPENMC_CMFD_SOLVER_H + +#include + +#include "xtensor/xtensor.hpp" + +namespace openmc { + +//=============================================================================== +// Global variables +//=============================================================================== + +// CSR format index pointer array of loss matrix +extern std::vector indptr; + +// CSR format index array of loss matrix +extern std::vector 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 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 \ No newline at end of file diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 75da02c0a..888978c59 100644 --- a/include/openmc/constants.h +++ b/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 diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 62f28ff6c..e3bdd6d4b 100644 --- a/openmc/cmfd.py +++ b/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] \ No newline at end of file + return current[2*l+1]/current[2*l] diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp new file mode 100644 index 000000000..07fdd3c3a --- /dev/null +++ b/src/cmfd_solver.cpp @@ -0,0 +1,339 @@ +//TODO remove +#include +#include + +#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 indptr; + +std::vector indices; + +int dim; + +double spectral; + +int nx, ny, nz, ng; + +xt::xtensor 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(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 \ No newline at end of file