max/TugJulia
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

apply format changes (LLVM)

Browse Source
This commit is contained in:
Max Lübke 2023-09-05 16:42:05 +02:00
parent 6e388f3d99
commit 7d05320f24
1 changed files with 358 additions and 335 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

301
src/BTCSv2.cpp
View File

@ -1,78 +1,88 @@
/** /**
* @file BTCSv2.cpp * @file BTCSv2.cpp
* @brief Implementation of heterogenous BTCS (backward time-centered space) solution * @brief Implementation of heterogenous BTCS (backward time-centered space)
* of diffusion equation in 1D and 2D space. Internally the alternating-direction * solution of diffusion equation in 1D and 2D space. Internally the
* implicit (ADI) method is used. Version 2, because Version 1 was an * alternating-direction implicit (ADI) method is used. Version 2, because
* implementation for the homogeneous BTCS solution. * Version 1 was an implementation for the homogeneous BTCS solution.
* *
*/ */
#include "FTCS.cpp" #include "FTCS.cpp"
#include <tug/Boundary.hpp>
#include <omp.h> #include <omp.h>
#include <tug/Boundary.hpp>
#ifdef WRITE_THOMAS_CSV #ifdef WRITE_THOMAS_CSV
#include<fstream> #include <fstream>
#include<string> #include <string>
#endif #endif
#define NUM_THREADS_BTCS 10 #define NUM_THREADS_BTCS 10
using namespace Eigen; using namespace Eigen;
// calculates coefficient for left boundary in constant case // calculates coefficient for left boundary in constant case
static tuple<double, double> calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) { static tuple<double, double>
calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) {
double centerCoeff; double centerCoeff;
double rightCoeff; double rightCoeff;
centerCoeff = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)) centerCoeff =
+ 2 * alpha(rowIndex,0)); 1 + sx * (calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)) +
rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)); 2 * alpha(rowIndex, 0));
rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1));
return {centerCoeff, rightCoeff}; return {centerCoeff, rightCoeff};
} }
// calculates coefficient for left boundary in closed case // calculates coefficient for left boundary in closed case
static tuple<double, double> calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) { static tuple<double, double>
calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) {
double centerCoeff; double centerCoeff;
double rightCoeff; double rightCoeff;
centerCoeff = 1 + sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)); centerCoeff =
rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)); 1 + sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1));
rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1));
return {centerCoeff, rightCoeff}; return {centerCoeff, rightCoeff};
} }
// calculates coefficient for right boundary in constant case // calculates coefficient for right boundary in constant case
static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, int n, double sx) { static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha,
int rowIndex, int n,
double sx) {
double leftCoeff; double leftCoeff;
double centerCoeff; double centerCoeff;
leftCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)); leftCoeff =
centerCoeff = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)) -sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
+ 2 * alpha(rowIndex,n)); centerCoeff =
1 + sx * (calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)) +
2 * alpha(rowIndex, n));
return {leftCoeff, centerCoeff}; return {leftCoeff, centerCoeff};
} }
// calculates coefficient for right boundary in closed case // calculates coefficient for right boundary in closed case
static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) { static tuple<double, double>
calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) {
double leftCoeff; double leftCoeff;
double centerCoeff; double centerCoeff;
leftCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)); leftCoeff =
centerCoeff = 1 + sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)); -sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
centerCoeff =
1 + sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
return {leftCoeff, centerCoeff}; return {leftCoeff, centerCoeff};
} }
// creates coefficient matrix for next time step from alphas in x-direction // creates coefficient matrix for next time step from alphas in x-direction
static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight, int numCols, int rowIndex, double sx) { static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha,
vector<BoundaryElement> &bcLeft,
vector<BoundaryElement> &bcRight,
int numCols, int rowIndex,
double sx) {
// square matrix of column^2 dimension for the coefficients // square matrix of column^2 dimension for the coefficients
SparseMatrix<double> cm(numCols, numCols); SparseMatrix<double> cm(numCols, numCols);
@ -81,41 +91,48 @@ static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryEl
// left column // left column
BC_TYPE type = bcLeft[rowIndex].getType(); BC_TYPE type = bcLeft[rowIndex].getType();
if (type == BC_TYPE_CONSTANT) { if (type == BC_TYPE_CONSTANT) {
auto [centerCoeffTop, rightCoeffTop] = calcLeftBoundaryCoeffConstant(alpha, rowIndex, sx); auto [centerCoeffTop, rightCoeffTop] =
cm.insert(0,0) = centerCoeffTop; calcLeftBoundaryCoeffConstant(alpha, rowIndex, sx);
cm.insert(0,1) = rightCoeffTop; cm.insert(0, 0) = centerCoeffTop;
cm.insert(0, 1) = rightCoeffTop;
} else if (type == BC_TYPE_CLOSED) { } else if (type == BC_TYPE_CLOSED) {
auto [centerCoeffTop, rightCoeffTop] = calcLeftBoundaryCoeffClosed(alpha, rowIndex, sx); auto [centerCoeffTop, rightCoeffTop] =
cm.insert(0,0) = centerCoeffTop; calcLeftBoundaryCoeffClosed(alpha, rowIndex, sx);
cm.insert(0,1) = rightCoeffTop; cm.insert(0, 0) = centerCoeffTop;
cm.insert(0, 1) = rightCoeffTop;
} else { } else {
throw_invalid_argument("Undefined Boundary Condition Type somewhere on Left or Top!"); throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Left or Top!");
} }
// inner columns // inner columns
int n = numCols-1; int n = numCols - 1;
for (int i = 1; i < n; i++) { for (int i = 1; i < n; i++) {
cm.insert(i,i-1) = -sx * calcAlphaIntercell(alpha(rowIndex,i-1), alpha(rowIndex,i)); cm.insert(i, i - 1) =
cm.insert(i,i) = 1 + sx * ( -sx * calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i));
calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex,i+1)) cm.insert(i, i) =
+ calcAlphaIntercell(alpha(rowIndex,i-1), alpha(rowIndex,i)) 1 +
) sx * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)) +
; calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i)));
cm.insert(i,i+1) = -sx * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex,i+1)); cm.insert(i, i + 1) =
-sx * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1));
} }
// right column // right column
type = bcRight[rowIndex].getType(); type = bcRight[rowIndex].getType();
if (type == BC_TYPE_CONSTANT) { if (type == BC_TYPE_CONSTANT) {
auto [leftCoeffBottom, centerCoeffBottom] = calcRightBoundaryCoeffConstant(alpha, rowIndex, n, sx); auto [leftCoeffBottom, centerCoeffBottom] =
cm.insert(n,n-1) = leftCoeffBottom; calcRightBoundaryCoeffConstant(alpha, rowIndex, n, sx);
cm.insert(n,n) = centerCoeffBottom; cm.insert(n, n - 1) = leftCoeffBottom;
cm.insert(n, n) = centerCoeffBottom;
} else if (type == BC_TYPE_CLOSED) { } else if (type == BC_TYPE_CLOSED) {
auto [leftCoeffBottom, centerCoeffBottom] = calcRightBoundaryCoeffClosed(alpha, rowIndex, n, sx); auto [leftCoeffBottom, centerCoeffBottom] =
cm.insert(n,n-1) = leftCoeffBottom; calcRightBoundaryCoeffClosed(alpha, rowIndex, n, sx);
cm.insert(n,n) = centerCoeffBottom; cm.insert(n, n - 1) = leftCoeffBottom;
cm.insert(n, n) = centerCoeffBottom;
} else { } else {
throw_invalid_argument("Undefined Boundary Condition Type somewhere on Right or Bottom!"); throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Right or Bottom!");
} }
cm.makeCompressed(); // important for Eigen solver cm.makeCompressed(); // important for Eigen solver
@ -123,81 +140,72 @@ static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryEl
return cm; return cm;
} }
// calculates explicity concentration at top boundary in constant case // calculates explicity concentration at top boundary in constant case
static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrations, static double calcExplicitConcentrationsTopBoundaryConstant(
MatrixXd &alpha, vector<BoundaryElement> &bcTop, int rowIndex, int i, double sy) { MatrixXd &concentrations, MatrixXd &alpha, vector<BoundaryElement> &bcTop,
int rowIndex, int i, double sy) {
double c; double c;
c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) *
* concentrations(rowIndex,i) concentrations(rowIndex, i) +
+ ( (1 -
1 - sy * ( sy * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) +
calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) 2 * alpha(rowIndex, i))) *
+ 2 * alpha(rowIndex,i) concentrations(rowIndex, i) +
) sy * alpha(rowIndex, i) * bcTop[i].getValue();
) * concentrations(rowIndex,i)
+ sy * alpha(rowIndex,i) * bcTop[i].getValue();
return c; return c;
} }
// calculates explicit concentration at top boundary in closed case // calculates explicit concentration at top boundary in closed case
static double calcExplicitConcentrationsTopBoundaryClosed(MatrixXd &concentrations, static double calcExplicitConcentrationsTopBoundaryClosed(
MatrixXd &alpha, int rowIndex, int i, double sy) { MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) {
double c; double c;
c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) *
* concentrations(rowIndex,i) concentrations(rowIndex, i) +
+ ( (1 -
1 - sy * ( sy * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)))) *
calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) concentrations(rowIndex, i);
)
) * concentrations(rowIndex,i);
return c; return c;
} }
// calculates explicit concentration at bottom boundary in constant case // calculates explicit concentration at bottom boundary in constant case
static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concentrations, static double calcExplicitConcentrationsBottomBoundaryConstant(
MatrixXd &alpha, vector<BoundaryElement> &bcBottom, int rowIndex, int i, double sy) { MatrixXd &concentrations, MatrixXd &alpha,
vector<BoundaryElement> &bcBottom, int rowIndex, int i, double sy) {
double c; double c;
c = sy * alpha(rowIndex,i) * bcBottom[i].getValue() c = sy * alpha(rowIndex, i) * bcBottom[i].getValue() +
+ ( (1 -
1 - sy * ( sy * (2 * alpha(rowIndex, i) +
2 * alpha(rowIndex,i) calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)))) *
+ calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i)) concentrations(rowIndex, i) +
) sy * calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)) *
) * concentrations(rowIndex,i) concentrations(rowIndex - 1, i);
+ sy * calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i))
* concentrations(rowIndex-1,i);
return c; return c;
} }
// calculates explicit concentration at bottom boundary in closed case // calculates explicit concentration at bottom boundary in closed case
static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentrations, static double calcExplicitConcentrationsBottomBoundaryClosed(
MatrixXd &alpha, int rowIndex, int i, double sy) { MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) {
double c; double c;
c = ( c = (1 -
1 - sy * ( sy * (+calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)))) *
+ calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i)) concentrations(rowIndex, i) +
) sy * calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)) *
) * concentrations(rowIndex,i) concentrations(rowIndex - 1, i);
+ sy * calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i))
* concentrations(rowIndex-1,i);
return c; return c;
} }
// creates a solution vector for next time step from the current state of
// creates a solution vector for next time step from the current state of concentrations // concentrations
static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY, static VectorXd createSolutionVector(
MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY,
vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight, vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight,
vector<BoundaryElement> &bcTop, vector<BoundaryElement> &bcBottom, vector<BoundaryElement> &bcTop, vector<BoundaryElement> &bcBottom,
int length, int rowIndex, double sx, double sy) { int length, int rowIndex, double sx, double sy) {
@ -207,19 +215,20 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
BC_TYPE type; BC_TYPE type;
// inner rows // inner rows
if (rowIndex > 0 && rowIndex < numRows-1) { if (rowIndex > 0 && rowIndex < numRows - 1) {
for (int i = 0; i < length; i++) { for (int i = 0; i < length; i++) {
sv(i) = sy * calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i)) sv(i) =
* concentrations(rowIndex+1,i) sy *
+ ( calcAlphaIntercell(alphaY(rowIndex, i), alphaY(rowIndex + 1, i)) *
1 - sy * ( concentrations(rowIndex + 1, i) +
calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i)) (1 - sy * (calcAlphaIntercell(alphaY(rowIndex, i),
+ calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i)) alphaY(rowIndex + 1, i)) +
) calcAlphaIntercell(alphaY(rowIndex - 1, i),
) * concentrations(rowIndex,i) alphaY(rowIndex, i)))) *
+ sy * calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i)) concentrations(rowIndex, i) +
* concentrations(rowIndex-1,i) sy *
; calcAlphaIntercell(alphaY(rowIndex - 1, i), alphaY(rowIndex, i)) *
concentrations(rowIndex - 1, i);
} }
} }
@ -228,43 +237,51 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
for (int i = 0; i < length; i++) { for (int i = 0; i < length; i++) {
type = bcTop[i].getType(); type = bcTop[i].getType();
if (type == BC_TYPE_CONSTANT) { if (type == BC_TYPE_CONSTANT) {
sv(i) = calcExplicitConcentrationsTopBoundaryConstant(concentrations, alphaY, bcTop, rowIndex, i, sy); sv(i) = calcExplicitConcentrationsTopBoundaryConstant(
concentrations, alphaY, bcTop, rowIndex, i, sy);
} else if (type == BC_TYPE_CLOSED) { } else if (type == BC_TYPE_CLOSED) {
sv(i) = calcExplicitConcentrationsTopBoundaryClosed(concentrations, alphaY, rowIndex, i, sy); sv(i) = calcExplicitConcentrationsTopBoundaryClosed(
concentrations, alphaY, rowIndex, i, sy);
} else { } else {
throw_invalid_argument("Undefined Boundary Condition Type somewhere on Left or Top!"); throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Left or Top!");
} }
} }
} }
// last row // last row
else if (rowIndex == numRows-1) { else if (rowIndex == numRows - 1) {
for (int i = 0; i < length; i++) { for (int i = 0; i < length; i++) {
type = bcBottom[i].getType(); type = bcBottom[i].getType();
if (type == BC_TYPE_CONSTANT) { if (type == BC_TYPE_CONSTANT) {
sv(i) = calcExplicitConcentrationsBottomBoundaryConstant(concentrations, alphaY, bcBottom, rowIndex, i, sy); sv(i) = calcExplicitConcentrationsBottomBoundaryConstant(
concentrations, alphaY, bcBottom, rowIndex, i, sy);
} else if (type == BC_TYPE_CLOSED) { } else if (type == BC_TYPE_CLOSED) {
sv(i) = calcExplicitConcentrationsBottomBoundaryClosed(concentrations, alphaY, rowIndex, i, sy); sv(i) = calcExplicitConcentrationsBottomBoundaryClosed(
concentrations, alphaY, rowIndex, i, sy);
} else { } else {
throw_invalid_argument("Undefined Boundary Condition Type somewhere on Right or Bottom!"); throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Right or Bottom!");
} }
} }
} }
// first column -> additional fixed concentration change from perpendicular dimension in constant bc case // first column -> additional fixed concentration change from perpendicular
// dimension in constant bc case
if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) { if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) {
sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft[rowIndex].getValue(); sv(0) += 2 * sx * alphaX(rowIndex, 0) * bcLeft[rowIndex].getValue();
} }
// last column -> additional fixed concentration change from perpendicular dimension in constant bc case // last column -> additional fixed concentration change from perpendicular
// dimension in constant bc case
if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) { if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) {
sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight[rowIndex].getValue(); sv(length - 1) +=
2 * sx * alphaX(rowIndex, length - 1) * bcRight[rowIndex].getValue();
} }
return sv; return sv;
} }
// solver for linear equation system; A corresponds to coefficient matrix, // solver for linear equation system; A corresponds to coefficient matrix,
// b to the solution vector // b to the solution vector
// use of EigenLU solver // use of EigenLU solver
@ -302,8 +319,8 @@ static VectorXd ThomasAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
// HACK: write CSV to file // HACK: write CSV to file
#ifdef WRITE_THOMAS_CSV #ifdef WRITE_THOMAS_CSV
#include<fstream> #include <fstream>
#include<string> #include <string>
static std::uint32_t file_index = 0; static std::uint32_t file_index = 0;
std::string file_name = "Thomas_" + std::to_string(file_index++) + ".csv"; std::string file_name = "Thomas_" + std::to_string(file_index++) + ".csv";
@ -316,7 +333,8 @@ static VectorXd ThomasAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
// iterate through all elements // iterate through all elements
for (std::size_t i = 0; i < n; i++) { for (std::size_t i = 0; i < n; i++) {
out_file << a_diag[i] << ", " << b_diag[i] << ", " << c_diag[i] << ", " << b[i] << "\n"; out_file << a_diag[i] << ", " << b_diag[i] << ", " << c_diag[i] << ", "
<< b[i] << "\n";
} }
out_file.close(); out_file.close();
@ -343,9 +361,10 @@ static VectorXd ThomasAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
return x_vec; return x_vec;
} }
// BTCS solution for 1D grid // BTCS solution for 1D grid
static void BTCS_1D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solverFunc) (SparseMatrix<double> &A, VectorXd &b)) { static void BTCS_1D(Grid &grid, Boundary &bc, double timestep,
VectorXd (*solverFunc)(SparseMatrix<double> &A,
VectorXd &b)) {
int length = grid.getLength(); int length = grid.getLength();
double sx = timestep / (grid.getDelta() * grid.getDelta()); double sx = timestep / (grid.getDelta() * grid.getDelta());
@ -360,29 +379,32 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solver
MatrixXd concentrations = grid.getConcentrations(); MatrixXd concentrations = grid.getConcentrations();
int rowIndex = 0; int rowIndex = 0;
A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx); // this is exactly same as in 2D A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex,
sx); // this is exactly same as in 2D
for (int i = 0; i < length; i++) { for (int i = 0; i < length; i++) {
b(i) = concentrations(0,i); b(i) = concentrations(0, i);
} }
if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) { if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) {
b(0) += 2 * sx * alpha(0,0) * bcLeft[0].getValue(); b(0) += 2 * sx * alpha(0, 0) * bcLeft[0].getValue();
} }
if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) { if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) {
b(length-1) += 2 * sx * alpha(0,length-1) * bcRight[0].getValue(); b(length - 1) += 2 * sx * alpha(0, length - 1) * bcRight[0].getValue();
} }
concentrations_t1 = solverFunc(A, b); concentrations_t1 = solverFunc(A, b);
for (int j = 0; j < length; j++) { for (int j = 0; j < length; j++) {
concentrations(0,j) = concentrations_t1(j); concentrations(0, j) = concentrations_t1(j);
} }
grid.setConcentrations(concentrations); grid.setConcentrations(concentrations);
} }
// BTCS solution for 2D grid // BTCS solution for 2D grid
static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solverFunc) (SparseMatrix<double> &A, VectorXd &b), int numThreads) { static void BTCS_2D(Grid &grid, Boundary &bc, double timestep,
VectorXd (*solverFunc)(SparseMatrix<double> &A,
VectorXd &b),
int numThreads) {
int rowMax = grid.getRow(); int rowMax = grid.getRow();
int colMax = grid.getCol(); int colMax = grid.getCol();
double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
@ -403,10 +425,9 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solver
MatrixXd concentrations = grid.getConcentrations(); MatrixXd concentrations = grid.getConcentrations();
#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
for (int i = 0; i < rowMax; i++) { for (int i = 0; i < rowMax; i++) {
A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx); A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx);
b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
bcTop, bcBottom, colMax, i, sx, sy); bcTop, bcBottom, colMax, i, sx, sy);
@ -423,7 +444,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solver
alphaX.transposeInPlace(); alphaX.transposeInPlace();
alphaY.transposeInPlace(); alphaY.transposeInPlace();
#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
for (int i = 0; i < colMax; i++) { for (int i = 0; i < colMax; i++) {
// swap alphas, boundary conditions and sx/sy for column-wise calculation // swap alphas, boundary conditions and sx/sy for column-wise calculation
A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy); A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy);
@ -440,7 +461,6 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solver
grid.setConcentrations(concentrations); grid.setConcentrations(concentrations);
} }
// entry point for EigenLU solver; differentiate between 1D and 2D grid // entry point for EigenLU solver; differentiate between 1D and 2D grid
static void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { static void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) {
if (grid.getDim() == 1) { if (grid.getDim() == 1) {
@ -448,17 +468,20 @@ static void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) {
} else if (grid.getDim() == 2) { } else if (grid.getDim() == 2) {
BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads); BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads);
} else { } else {
throw_invalid_argument("Error: Only 1- and 2-dimensional grids are defined!"); throw_invalid_argument(
"Error: Only 1- and 2-dimensional grids are defined!");
} }
} }
// entry point for Thomas algorithm solver; differentiate 1D and 2D grid // entry point for Thomas algorithm solver; differentiate 1D and 2D grid
static void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) { static void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep,
int numThreads) {
if (grid.getDim() == 1) { if (grid.getDim() == 1) {
BTCS_1D(grid, bc, timestep, ThomasAlgorithm); BTCS_1D(grid, bc, timestep, ThomasAlgorithm);
} else if (grid.getDim() == 2) { } else if (grid.getDim() == 2) {
BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads); BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads);
} else { } else {
throw_invalid_argument("Error: Only 1- and 2-dimensional grids are defined!"); throw_invalid_argument(
"Error: Only 1- and 2-dimensional grids are defined!");
} }
} }