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

feat: implementation of dynamic boundary conditions in FTCS

Browse Source
This commit is contained in:
philippun 2023-07-28 15:54:03 +02:00
parent 9703dba718
commit 0872639c54
4 changed files with 79 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
examples/FTCS_1D_proto_example.cpp
View File

@ -1,3 +1,4 @@
#include "tug/Boundary.hpp"
#include <tug/Simulation.hpp>
int main(int argc, char *argv[]) {
@ -18,7 +19,9 @@ int main(int argc, char *argv[]) {
// ******************
// create a boundary with constant values
Boundary bc = Boundary(grid, BC_TYPE_CONSTANT);
Boundary bc = Boundary(grid);
bc.setBoundarySideConstant(BC_SIDE_LEFT, 1);
bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1);
// ************************

9
examples/FTCS_2D_proto_example.cpp
View File

@ -44,7 +44,12 @@ int main(int argc, char *argv[]) {
// ******************
// create a boundary with constant values
Boundary bc = Boundary(grid, BC_TYPE_CONSTANT);
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);
// (optional) set boundary condition values for one side, e.g.:
// VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value
@ -68,7 +73,7 @@ int main(int argc, char *argv[]) {
simulation.setTimestep(0.1); // timestep
// (optional) set the number of iterations
simulation.setIterations(2);
simulation.setIterations(10000);
// (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);

5
src/Boundary.cpp
View File

@ -34,12 +34,15 @@ double BoundaryElement::getValue() {
Boundary::Boundary(Grid grid) : grid(grid) {
//probably to DEBUG assignment grid
if (grid.getDim() == 1) {
this->boundaries = vector<vector<BoundaryElement>>(2);
this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement());
this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement());
} else if (grid.getDim() == 2) {
this->boundaries = vector<vector<BoundaryElement>>(4);
this->boundaries[BC_SIDE_LEFT] = vector<BoundaryElement>(grid.getRow(), BoundaryElement());
this->boundaries[BC_SIDE_RIGHT] = vector<BoundaryElement>(grid.getRow(), BoundaryElement());
this->boundaries[BC_SIDE_TOP] = vector<BoundaryElement>(grid.getCol(), BoundaryElement());

82
src/FTCS.cpp
View File

@ -1,9 +1,11 @@
#include "TugUtils.hpp"
#include <cstddef>
#include <tug/Boundary.hpp>
#include <iostream>
using namespace std;
double calcAlphaIntercell(double alpha1, double alpha2, bool useHarmonic = false) {
if (useHarmonic) {
return 2 / ((1/alpha1) + (1/alpha2));
@ -57,7 +59,7 @@ double calcHorizontalChangeLeftBoundaryConstant(Grid grid, Boundary bc, int row,
+ 2 * grid.getAlphaX()(row,col)
)
* grid.getConcentrations()(row,col)
+ 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row);
+ 2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_LEFT, row);
return result;
}
@ -73,10 +75,21 @@ double calcHorizontalChangeLeftBoundaryClosed(Grid grid, int row, int col) {
}
double calcHorizontalChangeLeftBoundary(Grid grid, Boundary bc, int row, int col) {
if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) {
return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CLOSED) {
return calcHorizontalChangeLeftBoundaryClosed(grid, row, col);
} else {
throw_invalid_argument("Undefined Boundary Condition Type!");
}
}
double calcHorizontalChangeRightBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
double result =
2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_RIGHT, row)
- (
calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+ 2 * grid.getAlphaX()(row,col)
@ -99,6 +112,17 @@ double calcHorizontalChangeRightBoundaryClosed(Grid grid, int row, int col) {
}
double calcHorizontalChangeRightBoundary(Grid grid, Boundary bc, int row, int col) {
if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) {
return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CLOSED) {
return calcHorizontalChangeRightBoundaryClosed(grid, row, col);
} else {
throw_invalid_argument("Undefined Boundary Condition Type!");
}
}
double calcVerticalChangeTopBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
double result =
@ -109,7 +133,7 @@ double calcVerticalChangeTopBoundaryConstant(Grid grid, Boundary bc, int row, in
+ 2 * grid.getAlphaY()(row, col)
)
* grid.getConcentrations()(row, col)
+ 2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col);
+ 2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_TOP, col);
return result;
}
@ -125,10 +149,21 @@ double calcVerticalChangeTopBoundaryClosed(Grid grid, int row, int col) {
}
double calcVerticalChangeTopBoundary(Grid grid, Boundary bc, int row, int col) {
if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) {
return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) {
return calcVerticalChangeTopBoundaryClosed(grid, row, col);
} else {
throw_invalid_argument("Undefined Boundary Condition Type!");
}
}
double calcVerticalChangeBottomBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
double result =
2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col)
- (
calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col))
+ 2 * grid.getAlphaY()(row, col)
@ -151,6 +186,17 @@ double calcVerticalChangeBottomBoundaryClosed(Grid grid, int row, int col) {
}
double calcVerticalChangeBottomBoundary(Grid grid, Boundary bc, int row, int col) {
if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) {
return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) {
return calcVerticalChangeBottomBoundaryClosed(grid, row, col);
} else {
throw_invalid_argument("Undefined Boundary Condition Type!");
}
}
MatrixXd FTCS_1D(Grid grid, Boundary bc, double timestep) {
int colMax = grid.getCol();
double deltaCol = grid.getDeltaCol();
@ -175,7 +221,7 @@ MatrixXd FTCS_1D(Grid grid, Boundary bc, double timestep) {
concentrations_t1(row, col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
)
;
@ -185,7 +231,7 @@ MatrixXd FTCS_1D(Grid grid, Boundary bc, double timestep) {
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
calcHorizontalChangeRightBoundary(grid, bc, row, col)
)
;
@ -240,7 +286,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
concentrations_t1(row, col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
)
+ timestep / (deltaRow*deltaRow)
* (
@ -255,7 +301,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
calcHorizontalChangeRightBoundary(grid, bc, row, col)
)
+ timestep / (deltaRow*deltaRow)
* (
@ -271,7 +317,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep / (deltaRow*deltaRow)
* (
calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
calcVerticalChangeTopBoundary(grid, bc, row, col)
)
+ timestep / (deltaCol*deltaCol)
* (
@ -286,7 +332,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep / (deltaRow*deltaRow)
* (
calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
calcVerticalChangeBottomBoundary(grid, bc, row, col)
)
+ timestep / (deltaCol*deltaCol)
* (
@ -301,11 +347,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep/(deltaCol*deltaCol)
* (
calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
)
+ timestep/(deltaRow*deltaRow)
* (
calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
calcVerticalChangeTopBoundary(grid, bc, row, col)
)
;
@ -315,11 +361,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep/(deltaCol*deltaCol)
* (
calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
calcHorizontalChangeRightBoundary(grid, bc, row, col)
)
+ timestep/(deltaRow*deltaRow)
* (
calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
calcVerticalChangeTopBoundary(grid, bc, row, col)
)
;
@ -329,11 +375,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep/(deltaCol*deltaCol)
* (
calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
)
+ timestep/(deltaRow*deltaRow)
* (
calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
calcVerticalChangeBottomBoundary(grid, bc, row, col)
)
;
@ -343,11 +389,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep/(deltaCol*deltaCol)
* (
calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
calcHorizontalChangeRightBoundary(grid, bc, row, col)
)
+ timestep/(deltaRow*deltaRow)
* (
calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
calcVerticalChangeBottomBoundary(grid, bc, row, col)
)
;