mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 21:55:41 -04:00
Address @paulromano and @smharper comments
This commit is contained in:
Shikhar Kumar
parent
b4a2687af9
commit
2874b4a2f9
10 changed files with 782 additions and 959 deletions
|
|
@ -1,11 +1,16 @@
|
|||
#include <vector>
|
||||
#include <cmath>
|
||||
|
||||
#include "xtensor/xtensor.hpp"
|
||||
|
||||
#include "openmc/cmfd_solver.h"
|
||||
#include "openmc/error.h"
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/capi.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
namespace cmfd {
|
||||
|
||||
//==============================================================================
|
||||
// Global variables
|
||||
//==============================================================================
|
||||
|
|
@ -22,14 +27,30 @@ int nx, ny, nz, ng;
|
|||
|
||||
xt::xtensor<int, 2> indexmap;
|
||||
|
||||
} // namespace cmfd
|
||||
|
||||
//==============================================================================
|
||||
// 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 % cmfd::ng;
|
||||
i = cmfd::indexmap(irow/cmfd::ng, 0);
|
||||
j = cmfd::indexmap(irow/cmfd::ng, 1);
|
||||
k = cmfd::indexmap(irow/cmfd::ng, 2);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// 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)
|
||||
int get_diagonal_index(int row)
|
||||
{
|
||||
for (int j = cmfd::indptr[row]; j < cmfd::indptr[row+1]; j++) {
|
||||
if (cmfd::indices[j] == row)
|
||||
return j;
|
||||
}
|
||||
|
||||
|
|
@ -41,15 +62,16 @@ int get_diagonal_index(int row) {
|
|||
// 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;
|
||||
void set_indexmap(const int* coremap)
|
||||
{
|
||||
for (int z = 0; z < cmfd::nz; z++) {
|
||||
for (int y = 0; y < cmfd::ny; y++) {
|
||||
for (int x = 0; x < cmfd::nx; x++) {
|
||||
if (coremap[(z*cmfd::ny*cmfd::nx) + (y*cmfd::nx) + x] != CMFD_NOACCEL) {
|
||||
int counter = coremap[(z*cmfd::ny*cmfd::nx) + (y*cmfd::nx) + x];
|
||||
cmfd::indexmap(counter, 0) = x;
|
||||
cmfd::indexmap(counter, 1) = y;
|
||||
cmfd::indexmap(counter, 2) = z;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -60,23 +82,24 @@ void set_indexmap(int* coremap) {
|
|||
// CMFD_LINSOLVER_1G solves a one group CMFD linear system
|
||||
//==============================================================================
|
||||
|
||||
int cmfd_linsolver_1g(double* A_data, double* b, double* x, double tol) {
|
||||
int cmfd_linsolver_1g(const double* A_data, const 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);
|
||||
std::vector<double> tmpx {x, x+cmfd::dim};
|
||||
|
||||
// Perform red/black Gauss-Seidel iterations
|
||||
for (int irb = 0; irb < 2; irb++) {
|
||||
|
||||
// Loop around matrix rows
|
||||
for (int irow = 0; irow < dim; irow++) {
|
||||
for (int irow = 0; irow < cmfd::dim; irow++) {
|
||||
int g, i, j, k;
|
||||
matrix_to_indices(irow, g, i, j, k);
|
||||
|
||||
|
|
@ -88,10 +111,10 @@ int cmfd_linsolver_1g(double* A_data, double* b, double* x, double tol) {
|
|||
|
||||
// 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]];
|
||||
for (int icol = cmfd::indptr[irow]; icol < didx; icol++)
|
||||
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
|
||||
for (int icol = didx + 1; icol < cmfd::indptr[irow + 1]; icol++)
|
||||
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
|
||||
|
||||
// Solve for new x
|
||||
double x1 = (b[irow] - tmp1) / A_data[didx];
|
||||
|
|
@ -106,12 +129,12 @@ int cmfd_linsolver_1g(double* A_data, double* b, double* x, double tol) {
|
|||
}
|
||||
|
||||
// Check convergence
|
||||
err = std::sqrt(err / dim);
|
||||
err = std::sqrt(err / cmfd::dim);
|
||||
if (err < tol)
|
||||
return igs;
|
||||
|
||||
// Calculate new overrelaxation parameter
|
||||
w = 1.0/(1.0 - 0.25 * spectral * w);
|
||||
w = 1.0/(1.0 - 0.25 * cmfd::spectral * w);
|
||||
}
|
||||
|
||||
// Throw error, as max iterations met
|
||||
|
|
@ -125,23 +148,24 @@ int cmfd_linsolver_1g(double* A_data, double* b, double* x, double tol) {
|
|||
// CMFD_LINSOLVER_2G solves a two group CMFD linear system
|
||||
//==============================================================================
|
||||
|
||||
int cmfd_linsolver_2g(double* A_data, double* b, double* x, double tol) {
|
||||
int cmfd_linsolver_2g(const double* A_data, const 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);
|
||||
std::vector<double> tmpx {x, x+cmfd::dim};
|
||||
|
||||
// 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) {
|
||||
for (int irow = 0; irow < cmfd::dim; irow+=2) {
|
||||
int g, i, j, k;
|
||||
matrix_to_indices(irow, g, i, j, k);
|
||||
|
||||
|
|
@ -168,14 +192,14 @@ int cmfd_linsolver_2g(double* A_data, double* b, double* x, double tol) {
|
|||
// 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++)
|
||||
tmp2 += A_data[icol] * x[indices[icol]];
|
||||
for (int icol = cmfd::indptr[irow]; icol < d1idx; icol++)
|
||||
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
|
||||
for (int icol = cmfd::indptr[irow+1]; icol < d2idx-1; icol++)
|
||||
tmp2 += A_data[icol] * x[cmfd::indices[icol]];
|
||||
for (int icol = d1idx + 2; icol < cmfd::indptr[irow + 1]; icol++)
|
||||
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
|
||||
for (int icol = d2idx + 1; icol < cmfd::indptr[irow + 2]; icol++)
|
||||
tmp2 += A_data[icol] * x[cmfd::indices[icol]];
|
||||
|
||||
// Adjust with RHS vector
|
||||
tmp1 = b[irow] - tmp1;
|
||||
|
|
@ -196,12 +220,12 @@ int cmfd_linsolver_2g(double* A_data, double* b, double* x, double tol) {
|
|||
}
|
||||
|
||||
// Check convergence
|
||||
err = std::sqrt(err / dim);
|
||||
err = std::sqrt(err / cmfd::dim);
|
||||
if (err < tol)
|
||||
return igs;
|
||||
|
||||
// Calculate new overrelaxation parameter
|
||||
w = 1.0/(1.0 - 0.25 * spectral * w);
|
||||
w = 1.0/(1.0 - 0.25 * cmfd::spectral * w);
|
||||
}
|
||||
|
||||
// Throw error, as max iterations met
|
||||
|
|
@ -215,29 +239,30 @@ int cmfd_linsolver_2g(double* A_data, double* b, double* x, double tol) {
|
|||
// CMFD_LINSOLVER_NG solves a general CMFD linear system
|
||||
//==============================================================================
|
||||
|
||||
int cmfd_linsolver_ng(double* A_data, double* b, double* x, double tol) {
|
||||
int cmfd_linsolver_ng(const double* A_data, const 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);
|
||||
std::vector<double> tmpx {x, x+cmfd::dim};
|
||||
|
||||
// Loop around matrix rows
|
||||
for (int irow = 0; irow < dim; irow++) {
|
||||
for (int irow = 0; irow < cmfd::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]];
|
||||
for (int icol = cmfd::indptr[irow]; icol < didx; icol++)
|
||||
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
|
||||
for (int icol = didx + 1; icol < cmfd::indptr[irow + 1]; icol++)
|
||||
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
|
||||
|
||||
// Solve for new x
|
||||
double x1 = (b[irow] - tmp1) / A_data[didx];
|
||||
|
|
@ -251,12 +276,12 @@ int cmfd_linsolver_ng(double* A_data, double* b, double* x, double tol) {
|
|||
}
|
||||
|
||||
// Check convergence
|
||||
err = std::sqrt(err / dim);
|
||||
err = std::sqrt(err / cmfd::dim);
|
||||
if (err < tol)
|
||||
return igs;
|
||||
|
||||
// Calculate new overrelaxation parameter
|
||||
w = 1.0/(1.0 - 0.25 * spectral * w);
|
||||
w = 1.0/(1.0 - 0.25 * cmfd::spectral * w);
|
||||
}
|
||||
|
||||
// Throw error, as max iterations met
|
||||
|
|
@ -266,50 +291,40 @@ int cmfd_linsolver_ng(double* A_data, double* b, double* x, double tol) {
|
|||
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) {
|
||||
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)
|
||||
{
|
||||
// Store elements of indptr
|
||||
for (int i = 0; i < len_indptr; i++)
|
||||
openmc::indptr.push_back(indptr[i]);
|
||||
cmfd::indptr.push_back(indptr[i]);
|
||||
|
||||
// Store elements of indices
|
||||
for (int i = 0; i < n_elements; i++)
|
||||
openmc::indices.push_back(indices[i]);
|
||||
cmfd::indices.push_back(indices[i]);
|
||||
|
||||
// Set dimenion of CMFD problem and specral radius
|
||||
openmc::dim = dim;
|
||||
openmc::spectral = spectral;
|
||||
cmfd::dim = dim;
|
||||
cmfd::spectral = spectral;
|
||||
|
||||
// Set number of groups
|
||||
openmc::ng = cmfd_indices[3];
|
||||
cmfd::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];
|
||||
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
|
||||
openmc::indexmap.resize({static_cast<size_t>(dim), 3});
|
||||
cmfd::indexmap.resize({static_cast<size_t>(dim), 3});
|
||||
set_indexmap(map);
|
||||
}
|
||||
}
|
||||
|
|
@ -320,8 +335,10 @@ void openmc_initialize_linsolver(int* indptr, int len_indptr, int* indices,
|
|||
//==============================================================================
|
||||
|
||||
extern "C"
|
||||
int openmc_run_linsolver(double* A_data, double* b, double* x, double tol) {
|
||||
switch (ng) {
|
||||
int openmc_run_linsolver(const double* A_data, const double* b, double* x,
|
||||
double tol)
|
||||
{
|
||||
switch (cmfd::ng) {
|
||||
case 1:
|
||||
return cmfd_linsolver_1g(A_data, b, x, tol);
|
||||
case 2:
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue