From f1b5138bccff6f67f3465f3ccdbabf9944f80cdc Mon Sep 17 00:00:00 2001 From: philippun Date: Fri, 11 Aug 2023 12:19:33 +0200 Subject: [PATCH] implementation of closed 1D and 2D cases --- examples/BTCS_2D_proto_example.cpp | 12 +- src/BTCSv2.cpp | 235 +++++++++++++++++++++++------ 2 files changed, 195 insertions(+), 52 deletions(-) diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index 0129c36..1114596 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -35,10 +35,14 @@ int main(int argc, char *argv[]) { // create a boundary with constant values Boundary bc = Boundary(grid); - bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); - bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); - bc.setBoundarySideConstant(BC_SIDE_TOP, 0); - bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); + bc.setBoundarySideClosed(BC_SIDE_LEFT); + bc.setBoundarySideClosed(BC_SIDE_RIGHT); + bc.setBoundarySideClosed(BC_SIDE_TOP); + bc.setBoundarySideClosed(BC_SIDE_BOTTOM); + // bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); + // bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + // bc.setBoundarySideConstant(BC_SIDE_TOP, 0); + // bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); // (optional) set boundary condition values for one side, e.g.: diff --git a/src/BTCSv2.cpp b/src/BTCSv2.cpp index e63ab60..3fc5c69 100644 --- a/src/BTCSv2.cpp +++ b/src/BTCSv2.cpp @@ -4,16 +4,71 @@ using namespace Eigen; -static SparseMatrix createCoeffMatrix(MatrixXd &alpha, int numCols, int rowIndex, double sx) { +static tuple calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, double &sx) { + double centerCoeff; + double rightCoeff; + + centerCoeff = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)) + + 2 * alpha(rowIndex,0)); + rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)); + + return {centerCoeff, rightCoeff}; +} + + +static tuple calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, double &sx) { + double centerCoeff; + double rightCoeff; + + centerCoeff = 1 + sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)); + rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)); + + return {centerCoeff, rightCoeff}; +} + + +static tuple calcRightBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, int &n, double &sx) { + double leftCoeff; + double centerCoeff; + + leftCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)); + centerCoeff = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)) + + 2 * alpha(rowIndex,n)); + + return {leftCoeff, centerCoeff}; +} + + +static tuple calcRightBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, int &n, double &sx) { + double leftCoeff; + double centerCoeff; + + leftCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)); + centerCoeff = 1 + sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)); + + return {leftCoeff, centerCoeff}; +} + + +static SparseMatrix createCoeffMatrix(MatrixXd &alpha, vector bcLeft, vector bcRight, int numCols, int rowIndex, double sx) { // square matrix of column^2 dimension for the coefficients SparseMatrix cm(numCols, numCols); cm.reserve(VectorXi::Constant(numCols, 3)); // left column - cm.insert(0,0) = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)) - + 2 * alpha(rowIndex,0)); - cm.insert(0,1) = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)); + BC_TYPE type = bcLeft[rowIndex].getType(); + if (type == BC_TYPE_CONSTANT) { + auto [centerCoeffTop, rightCoeffTop] = calcLeftBoundaryCoeffConstant(alpha, rowIndex, sx); + cm.insert(0,0) = centerCoeffTop; + cm.insert(0,1) = rightCoeffTop; + } else if (type == BC_TYPE_CLOSED) { + auto [centerCoeffTop, rightCoeffTop] = calcLeftBoundaryCoeffClosed(alpha, rowIndex, sx); + cm.insert(0,0) = centerCoeffTop; + cm.insert(0,1) = rightCoeffTop; + } else { + throw_invalid_argument("Undefined Boundary Condition Type somewhere on Left or Top!"); + } // inner columns int n = numCols-1; @@ -28,22 +83,101 @@ static SparseMatrix createCoeffMatrix(MatrixXd &alpha, int numCols, int } // right column - cm.insert(n,n-1) = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)); - cm.insert(n,n) = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)) - + 2 * alpha(rowIndex,n)); + type = bcRight[rowIndex].getType(); + if (type == BC_TYPE_CONSTANT) { + auto [leftCoeffBottom, centerCoeffBottom] = calcRightBoundaryCoeffConstant(alpha, rowIndex, n, sx); + cm.insert(n,n-1) = leftCoeffBottom; + cm.insert(n,n) = centerCoeffBottom; + } else if (type == BC_TYPE_CLOSED) { + auto [leftCoeffBottom, centerCoeffBottom] = calcRightBoundaryCoeffClosed(alpha, rowIndex, n, sx); + cm.insert(n,n-1) = leftCoeffBottom; + cm.insert(n,n) = centerCoeffBottom; + } else { + throw_invalid_argument("Undefined Boundary Condition Type somewhere on Right or Bottom!"); + } - cm.makeCompressed(); + cm.makeCompressed(); // important for Eigen solver return cm; } +static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrations, + MatrixXd &alpha, vector &bcTop, int &rowIndex, int &i, double &sy) { + double c; + + c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) + * concentrations(rowIndex,i) + + ( + 1 - sy * ( + calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) + + 2 * alpha(rowIndex,i) + ) + ) * concentrations(rowIndex,i) + + sy * alpha(rowIndex,i) * bcTop[i].getValue(); + + return c; +} + + +static double calcExplicitConcentrationsTopBoundaryClosed(MatrixXd &concentrations, + MatrixXd &alpha, int &rowIndex, int &i, double &sy) { + double c; + + c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) + * concentrations(rowIndex,i) + + ( + 1 - sy * ( + calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) + ) + ) * concentrations(rowIndex,i); + + return c; +} + + +static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concentrations, + MatrixXd &alpha, vector &bcBottom, int &rowIndex, int &i, double &sy) { + double c; + + c = sy * alpha(rowIndex,i) * bcBottom[i].getValue() + + ( + 1 - sy * ( + 2 * 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-1,i); + + return c; +} + + +static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentrations, + MatrixXd &alpha, int &rowIndex, int &i, double &sy) { + double c; + + c = ( + 1 - sy * ( + + calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i)) + ) + ) * concentrations(rowIndex,i) + + sy * calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i)) + * concentrations(rowIndex-1,i); + + return c; +} + + static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY, - VectorXd &bcLeft, VectorXd &bcRight, VectorXd &bcTop, - VectorXd &bcBottom, int length, int rowIndex, double sx, double sy) { + vector &bcLeft, vector &bcRight, + vector &bcTop, vector &bcBottom, + int length, int rowIndex, double sx, double sy) { VectorXd sv(length); int numRows = concentrations.rows(); + BC_TYPE type; // inner rows if (rowIndex > 0 && rowIndex < numRows-1) { @@ -65,40 +199,40 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, // first row if (rowIndex == 0) { for (int i = 0; i < length; i++) { - sv(i) = sy * calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i)) - * concentrations(rowIndex,i) - + ( - 1 - sy * ( - calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i)) - + 2 * alphaY(rowIndex,i) - ) - ) * concentrations(rowIndex,i) - + sy * alphaY(rowIndex,i) * bcTop(i) - ; + type = bcTop[i].getType(); + if (type == BC_TYPE_CONSTANT) { + sv(i) = calcExplicitConcentrationsTopBoundaryConstant(concentrations, alphaY, bcTop, rowIndex, i, sy); + } else if (type == BC_TYPE_CLOSED) { + sv(i) = calcExplicitConcentrationsTopBoundaryClosed(concentrations, alphaY, rowIndex, i, sy); + } else { + throw_invalid_argument("Undefined Boundary Condition Type somewhere on Left or Top!"); + } } } // last row if (rowIndex == numRows-1) { for (int i = 0; i < length; i++) { - sv(i) = sy * alphaY(rowIndex,i) * bcBottom(i) - + ( - 1 - sy * ( - 2 * alphaY(rowIndex,i) - + calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i)) - ) - ) * concentrations(rowIndex,i) - + sy * calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i)) - * concentrations(rowIndex-1,i) - ; + type = bcBottom[i].getType(); + if (type == BC_TYPE_CONSTANT) { + sv(i) = calcExplicitConcentrationsBottomBoundaryConstant(concentrations, alphaY, bcBottom, rowIndex, i, sy); + } else if (type == BC_TYPE_CLOSED) { + sv(i) = calcExplicitConcentrationsBottomBoundaryClosed(concentrations, alphaY, rowIndex, i, sy); + } else { + throw_invalid_argument("Undefined Boundary Condition Type somewhere on Right or Bottom!"); + } } } - // first column -> additional fixed concentration change from perpendicular dimension - sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft(rowIndex); + // first column -> additional fixed concentration change from perpendicular dimension in constant bc case + if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) { + sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft[rowIndex].getValue(); + } - // last column -> additional fixed concentration change from perpendicular dimension - sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight(rowIndex); + // last column -> additional fixed concentration change from perpendicular dimension in constant bc case + if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) { + sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight[rowIndex].getValue(); + } return sv; } @@ -124,16 +258,20 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double ×tep) { VectorXd b(length); MatrixXd alpha = grid.getAlpha(); - VectorXd bcLeft = bc.getBoundarySideValues(BC_SIDE_LEFT); - VectorXd bcRight = bc.getBoundarySideValues(BC_SIDE_RIGHT); + vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); MatrixXd concentrations = grid.getConcentrations(); - A = createCoeffMatrix(alpha, length, 0, sx); // this is exactly same as in 2D + A = createCoeffMatrix(alpha, bcLeft, bcRight, length, 0, sx); // this is exactly same as in 2D for (int i = 0; i < length; i++) { - b(i) = concentrations(0,i); // TODO check if this is really the only thing on right hand side of equation in 1D? + b(i) = concentrations(0,i); + } + if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) { + b(0) += 2 * sx * alpha(0,0) * bcLeft[0].getValue(); + } + if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) { + b(length-1) += 2 * sx * alpha(0,length-1) * bcRight[0].getValue(); } - b(0) += 2 * sx * alpha(0,0) * bcLeft(0); - b(length-1) += 2 * sx * alpha(0,length-1) * bcRight(0); concentrations_t1 = solve(A, b); @@ -160,21 +298,21 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double ×tep) { MatrixXd alphaX = grid.getAlphaX(); MatrixXd alphaY = grid.getAlphaY(); - VectorXd bcLeft = bc.getBoundarySideValues(BC_SIDE_LEFT); - VectorXd bcRight = bc.getBoundarySideValues(BC_SIDE_RIGHT); - VectorXd bcTop = bc.getBoundarySideValues(BC_SIDE_TOP); - VectorXd bcBottom = bc.getBoundarySideValues(BC_SIDE_BOTTOM); + vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + vector bcTop = bc.getBoundarySide(BC_SIDE_TOP); + vector bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); MatrixXd concentrations = grid.getConcentrations(); for (int i = 0; i < rowMax; i++) { - A = createCoeffMatrix(alphaX, colMax, i, sx); + A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx); b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, colMax, i, sx, sy); row_t1 = solve(A, b); for (int j = 0; j < colMax; j++) { - concentrations_t1(i,j) = row_t1(j); + concentrations_t1(i,j) = row_t1(j); // TODO Eigen seq ?? } } @@ -184,13 +322,14 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double ×tep) { alphaY.transposeInPlace(); for (int i = 0; i < colMax; i++) { - A = createCoeffMatrix(alphaY, rowMax, i, sy); + // swap alphas, boundary conditions and sx/sy for column-wise calculation + A = createCoeffMatrix(alphaY, bcLeft, bcRight, rowMax, i, sy); b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, bcLeft, bcRight, rowMax, i, sy, sx); row_t1 = solve(A, b); for (int j = 0; j < rowMax; j++) { - concentrations(i,j) = row_t1(j); + concentrations(i,j) = row_t1(j); // TODO Eigen seq ?? } } concentrations.transposeInPlace();