max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-13 17:38:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

refactor: inline Coeff Functions

Browse Source
This commit is contained in:
Max Lübke 2023-09-14 15:02:50 +02:00
parent 0eba63f875
commit a0d835e243
2 changed files with 123 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
src/BTCS.cpp
View File

@ -20,42 +20,35 @@
#endif #endif
// calculates coefficient for left boundary in constant case // calculates coefficient for left boundary in constant case
static std::tuple<double, double> static inline std::tuple<double, double>
calcLeftBoundaryCoeffConstant(Eigen::MatrixXd &alpha, int rowIndex, double sx) { calcLeftBoundaryCoeffConstant(Eigen::MatrixXd &alpha, int rowIndex, double sx) {
double centerCoeff; const double centerCoeff =
double rightCoeff;
centerCoeff =
1 + sx * (calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)) + 1 + sx * (calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)) +
2 * alpha(rowIndex, 0)); 2 * alpha(rowIndex, 0));
rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)); const double 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 std::tuple<double, double> static inline std::tuple<double, double>
calcLeftBoundaryCoeffClosed(Eigen::MatrixXd &alpha, int rowIndex, double sx) { calcLeftBoundaryCoeffClosed(Eigen::MatrixXd &alpha, int rowIndex, double sx) {
double centerCoeff; const double centerCoeff =
double rightCoeff;
centerCoeff =
1 + 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)); const double 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 std::tuple<double, double> static inline std::tuple<double, double>
calcRightBoundaryCoeffConstant(Eigen::MatrixXd &alpha, int rowIndex, int n, calcRightBoundaryCoeffConstant(Eigen::MatrixXd &alpha, int rowIndex, int n,
double sx) { double sx) {
double leftCoeff; const double leftCoeff =
double centerCoeff;
leftCoeff =
-sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)); -sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
centerCoeff = const double centerCoeff =
1 + sx * (calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)) + 1 + sx * (calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)) +
2 * alpha(rowIndex, n)); 2 * alpha(rowIndex, n));
@ -63,15 +56,12 @@ calcRightBoundaryCoeffConstant(Eigen::MatrixXd &alpha, int rowIndex, int n,
} }
// calculates coefficient for right boundary in closed case // calculates coefficient for right boundary in closed case
static std::tuple<double, double> static inline std::tuple<double, double>
calcRightBoundaryCoeffClosed(Eigen::MatrixXd &alpha, int rowIndex, int n, calcRightBoundaryCoeffClosed(Eigen::MatrixXd &alpha, int rowIndex, int n,
double sx) { double sx) {
double leftCoeff; const double leftCoeff =
double centerCoeff;
leftCoeff =
-sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)); -sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
centerCoeff = const double centerCoeff =
1 + sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)); 1 + sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
return {leftCoeff, centerCoeff}; return {leftCoeff, centerCoeff};

91
src/FTCS.cpp
View File

@ -19,9 +19,9 @@
#endif #endif
// calculates horizontal change on one cell independent of boundary type // calculates horizontal change on one cell independent of boundary type
static double calcHorizontalChange(Grid &grid, int &row, int &col) { static inline double calcHorizontalChange(Grid &grid, int &row, int &col) {
double result = calcAlphaIntercell(grid.getAlphaX()(row, col + 1), return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) * grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col + 1) - grid.getConcentrations()(row, col + 1) -
(calcAlphaIntercell(grid.getAlphaX()(row, col + 1), (calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
@ -32,14 +32,12 @@ static double calcHorizontalChange(Grid &grid, int &row, int &col) {
calcAlphaIntercell(grid.getAlphaX()(row, col - 1), calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) * grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col - 1); grid.getConcentrations()(row, col - 1);
return result;
} }
// calculates vertical change on one cell independent of boundary type // calculates vertical change on one cell independent of boundary type
static double calcVerticalChange(Grid &grid, int &row, int &col) { static inline double calcVerticalChange(Grid &grid, int &row, int &col) {
double result = calcAlphaIntercell(grid.getAlphaY()(row + 1, col), return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) * grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row + 1, col) - grid.getConcentrations()(row + 1, col) -
(calcAlphaIntercell(grid.getAlphaY()(row + 1, col), (calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
@ -50,15 +48,15 @@ static double calcVerticalChange(Grid &grid, int &row, int &col) {
calcAlphaIntercell(grid.getAlphaY()(row - 1, col), calcAlphaIntercell(grid.getAlphaY()(row - 1, col),
grid.getAlphaY()(row, col)) * grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row - 1, col); grid.getConcentrations()(row - 1, col);
return result;
} }
// calculates horizontal change on one cell with a constant left boundary // calculates horizontal change on one cell with a constant left boundary
static double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, Boundary &bc, static inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid,
int &row, int &col) { Boundary &bc,
int &row,
int &col) {
double result = calcAlphaIntercell(grid.getAlphaX()(row, col + 1), return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) * grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col + 1) - grid.getConcentrations()(row, col + 1) -
(calcAlphaIntercell(grid.getAlphaX()(row, col + 1), (calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
@ -67,24 +65,20 @@ static double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, Boundary &bc,
grid.getConcentrations()(row, col) + grid.getConcentrations()(row, col) +
2 * grid.getAlphaX()(row, col) * 2 * grid.getAlphaX()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_LEFT, row); bc.getBoundaryElementValue(BC_SIDE_LEFT, row);
return result;
} }
// calculates horizontal change on one cell with a closed left boundary // calculates horizontal change on one cell with a closed left boundary
static double calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, static inline double
int &col) { calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) {
double result = calcAlphaIntercell(grid.getAlphaX()(row, col + 1), return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) * grid.getAlphaX()(row, col)) *
(grid.getConcentrations()(row, col + 1) - (grid.getConcentrations()(row, col + 1) -
grid.getConcentrations()(row, col)); grid.getConcentrations()(row, col));
return result;
} }
// checks boundary condition type for a cell on the left edge of grid // checks boundary condition type for a cell on the left edge of grid
static double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, static inline double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc,
int &row, int &col) { int &row, int &col) {
if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) { if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) {
return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col); return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col);
@ -96,11 +90,12 @@ static double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc,
} }
// calculates horizontal change on one cell with a constant right boundary // calculates horizontal change on one cell with a constant right boundary
static double calcHorizontalChangeRightBoundaryConstant(Grid &grid, static inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid,
Boundary &bc, int &row, Boundary &bc,
int &row,
int &col) { int &col) {
double result = 2 * grid.getAlphaX()(row, col) * return 2 * grid.getAlphaX()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) - bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) -
(calcAlphaIntercell(grid.getAlphaX()(row, col - 1), (calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) + grid.getAlphaX()(row, col)) +
@ -109,24 +104,20 @@ static double calcHorizontalChangeRightBoundaryConstant(Grid &grid,
calcAlphaIntercell(grid.getAlphaX()(row, col - 1), calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) * grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col - 1); grid.getConcentrations()(row, col - 1);
return result;
} }
// calculates horizontal change on one cell with a closed right boundary // calculates horizontal change on one cell with a closed right boundary
static double calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, static inline double
int &col) { calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) {
double result = -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1), return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) * grid.getAlphaX()(row, col)) *
(grid.getConcentrations()(row, col) - (grid.getConcentrations()(row, col) -
grid.getConcentrations()(row, col - 1))); grid.getConcentrations()(row, col - 1)));
return result;
} }
// checks boundary condition type for a cell on the right edge of grid // checks boundary condition type for a cell on the right edge of grid
static double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, static inline double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc,
int &row, int &col) { int &row, int &col) {
if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) { if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) {
return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col); return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col);
@ -138,10 +129,11 @@ static double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc,
} }
// calculates vertical change on one cell with a constant top boundary // calculates vertical change on one cell with a constant top boundary
static double calcVerticalChangeTopBoundaryConstant(Grid &grid, Boundary &bc, static inline double calcVerticalChangeTopBoundaryConstant(Grid &grid,
Boundary &bc,
int &row, int &col) { int &row, int &col) {
double result = calcAlphaIntercell(grid.getAlphaY()(row + 1, col), return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) * grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row + 1, col) - grid.getConcentrations()(row + 1, col) -
(calcAlphaIntercell(grid.getAlphaY()(row + 1, col), (calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
@ -150,25 +142,21 @@ static double calcVerticalChangeTopBoundaryConstant(Grid &grid, Boundary &bc,
grid.getConcentrations()(row, col) + grid.getConcentrations()(row, col) +
2 * grid.getAlphaY()(row, col) * 2 * grid.getAlphaY()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_TOP, col); bc.getBoundaryElementValue(BC_SIDE_TOP, col);
return result;
} }
// calculates vertical change on one cell with a closed top boundary // calculates vertical change on one cell with a closed top boundary
static double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, static inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row,
int &col) { int &col) {
double result = calcAlphaIntercell(grid.getAlphaY()(row + 1, col), return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getConcentrations()(row, col)) * grid.getConcentrations()(row, col)) *
(grid.getConcentrations()(row + 1, col) - (grid.getConcentrations()(row + 1, col) -
grid.getConcentrations()(row, col)); grid.getConcentrations()(row, col));
return result;
} }
// checks boundary condition type for a cell on the top edge of grid // checks boundary condition type for a cell on the top edge of grid
static double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, int &row, static inline double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc,
int &col) { int &row, int &col) {
if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) { if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) {
return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col); return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) { } else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) {
@ -179,10 +167,12 @@ static double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, int &row,
} }
// calculates vertical change on one cell with a constant bottom boundary // calculates vertical change on one cell with a constant bottom boundary
static double calcVerticalChangeBottomBoundaryConstant(Grid &grid, Boundary &bc, static inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid,
int &row, int &col) { Boundary &bc,
int &row,
int &col) {
double result = 2 * grid.getAlphaY()(row, col) * return 2 * grid.getAlphaY()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) - bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) -
(calcAlphaIntercell(grid.getAlphaY()(row, col), (calcAlphaIntercell(grid.getAlphaY()(row, col),
grid.getAlphaY()(row - 1, col)) + grid.getAlphaY()(row - 1, col)) +
@ -191,24 +181,20 @@ static double calcVerticalChangeBottomBoundaryConstant(Grid &grid, Boundary &bc,
calcAlphaIntercell(grid.getAlphaY()(row, col), calcAlphaIntercell(grid.getAlphaY()(row, col),
grid.getAlphaY()(row - 1, col)) * grid.getAlphaY()(row - 1, col)) *
grid.getConcentrations()(row - 1, col); grid.getConcentrations()(row - 1, col);
return result;
} }
// calculates vertical change on one cell with a closed bottom boundary // calculates vertical change on one cell with a closed bottom boundary
static double calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, static inline double
int &col) { calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) {
double result = -(calcAlphaIntercell(grid.getAlphaY()(row, col), return -(calcAlphaIntercell(grid.getAlphaY()(row, col),
grid.getAlphaY()(row - 1, col)) * grid.getAlphaY()(row - 1, col)) *
(grid.getConcentrations()(row, col) - (grid.getConcentrations()(row, col) -
grid.getConcentrations()(row - 1, col))); grid.getConcentrations()(row - 1, col)));
return result;
} }
// checks boundary condition type for a cell on the bottom edge of grid // checks boundary condition type for a cell on the bottom edge of grid
static double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, static inline double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc,
int &row, int &col) { int &row, int &col) {
if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) { if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) {
return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col); return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col);
@ -265,7 +251,8 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep,
double deltaCol = grid.getDeltaCol(); double deltaCol = grid.getDeltaCol();
// matrix for concentrations at time t+1 // matrix for concentrations at time t+1
Eigen::MatrixXd concentrations_t1 = Eigen::MatrixXd::Constant(rowMax, colMax, 0); Eigen::MatrixXd concentrations_t1 =
Eigen::MatrixXd::Constant(rowMax, colMax, 0);
// inner cells // inner cells
// these are independent of the boundary condition type // these are independent of the boundary condition type