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

implemented Simulation with FTCS and tried a first example

Browse Source
This commit is contained in:
philippun 2023-07-18 17:14:16 +02:00
parent 0a9b58e8ff
commit 542601fdcd
8 changed files with 94 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
examples/CMakeLists.txt
View File

@ -1,7 +1,9 @@
add_executable(first_example first_example.cpp)
add_executable(second_example second_example.cpp)
add_executable(boundary_example1D boundary_example1D.cpp)
add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp)
target_link_libraries(first_example tug)
target_link_libraries(second_example tug)
target_link_libraries(boundary_example1D tug)
target_link_libraries(boundary_example1D tug)
target_link_libraries(FTCS_2D_proto_example tug)

11
examples/FTCS_2D_proto_example.cpp Normal file
View File

@ -0,0 +1,11 @@
#include <tug/Simulation.hpp>
int main(int argc, char *argv[]) {
Grid grid = Grid(20,20);
Boundary bc = Boundary(grid, BC_TYPE_CONSTANT);
Simulation simulation = Simulation(grid, bc, FTCS_APPROACH);
simulation.run();
}

2
include/tug/Boundary.hpp
View File

@ -50,7 +50,7 @@ class Boundary {
* @param side
* @return auto
*/
auto getBoundaryConditionValue(BC_SIDE side);
VectorXd getBoundaryConditionValue(BC_SIDE side);
private:

2
include/tug/Simulation.hpp
View File

@ -62,7 +62,7 @@ class Simulation {
*
* @return auto
*/
auto run();
void run();
private:

2
src/Boundary.cpp
View File

@ -51,7 +51,7 @@ void Boundary::setBoundaryConditionValue(BC_SIDE side, VectorXd &values) {
}
}
auto Boundary::getBoundaryConditionValue(BC_SIDE side) {
VectorXd Boundary::getBoundaryConditionValue(BC_SIDE side) {
switch (side) {
case BC_SIDE_LEFT:
return this->left;

2
src/CMakeLists.txt
View File

@ -1,4 +1,4 @@
add_library(tug BTCS.cpp BoundaryCondition.cpp Solver.cpp)
add_library(tug BTCS.cpp BoundaryCondition.cpp Solver.cpp Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp)
target_link_libraries(tug Eigen3::Eigen)

80
src/FTCS.cpp
View File

@ -23,7 +23,7 @@ auto FTCS_constant(Grid &grid, Boundary &bc, double timestep) {
for (int row = 1; row < rowMax-1; row++) {
for (int col = 1; col < colMax-1; col++) {
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep / deltaRow*deltaRow * (
+ timestep / (deltaRow*deltaRow) * (
calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
* grid.getConcentrations()(row+1,col)
- (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
@ -32,7 +32,7 @@ auto FTCS_constant(Grid &grid, Boundary &bc, double timestep) {
+ calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
* grid.getConcentrations()(row-1,col)
)
- timestep / deltaCol*deltaCol * (
- timestep / (deltaCol*deltaCol) * (
calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
* grid.getConcentrations()(row,col+1)
- (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
@ -46,17 +46,87 @@ auto FTCS_constant(Grid &grid, Boundary &bc, double timestep) {
// boundary conditions
// left without corners / looping over rows
int col = 0;
for (int row = 1; row < rowMax-1; row++) {
concentrations_t1(row, col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol)
* (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
* grid.getConcentrations()(row,col+1)
- (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
+ 2 * grid.getAlphaX()(row,col)) * grid.getConcentrations()(row,col)
+ 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row, 1))
+ timestep / (deltaRow*deltaRow)
* (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
* grid.getConcentrations()(row+1,col)
- (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+ calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
* grid.getConcentrations()(row,col)
+ calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getConcentrations()(row,col))
* grid.getConcentrations()(row-1,col));
}
return 0;
// right without corners / looping over columns
col = colMax-1;
for (int row = 1; row < rowMax-1; row++) {
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol)
* (2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row, 1)
- (calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+ 2 * grid.getAlphaX()(row,col)) + 2 * grid.getAlphaX()(row,col)
* grid.getConcentrations()(row,col)
+ calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
* grid.getConcentrations()(row,col-1))
+ timestep / (deltaRow*deltaRow)
* (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
* grid.getConcentrations()(row+1,col)
- (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+ calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
* grid.getConcentrations()(row,col)
+ calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
* grid.getConcentrations()(row-1,col));
}
// top without corners / looping over cols
for(int col=1; col<colMax-1;col++){
int row = 0;
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep/(grid.getDeltaRow()*grid.getDeltaRow()) * (calc_alpha_intercell(grid.getAlphaY()(1, col), grid.getAlphaY()(0, col)) * grid.getConcentrations()(1,col)
- (calc_alpha_intercell(grid.getAlphaY()(1, col), grid.getAlphaY()(0, col)) + 2 * grid.getAlphaY()(0, col)) * grid.getConcentrations()(0, col)
+ 2 * grid.getAlphaY()(0, col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(1, col))
+ timestep/(grid.getDeltaCol()*grid.getDeltaCol()) * (calc_alpha_intercell(grid.getAlphaX()(0, col+1), grid.getAlphaX()(0, col)) * grid.getConcentrations()(0, col+1)
- (calc_alpha_intercell(grid.getAlphaX()(0, col+1), grid.getAlphaX()(0, col)) + calc_alpha_intercell(grid.getAlphaX()(0, col-1), grid.getAlphaX()(0, col))) * grid.getConcentrations()(0, col)
+ calc_alpha_intercell(grid.getAlphaX()(0, col-1), grid.getAlphaX()(0, col)) * grid.getConcentrations()(0, col-1));
}
// bottom without corners / looping over cols
int row = rowMax-1;
for(int col=1; row<colMax-1;col++){
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep/(grid.getDeltaRow()*grid.getDeltaRow()) * (2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(1, col)
- (calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) + 2 * grid.getAlphaY()(row, col)) * grid.getConcentrations()(row, col)
+ calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)))
+ timestep/(grid.getDeltaCol()*grid.getDeltaCol()) * (calc_alpha_intercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) * grid.getConcentrations()(row, col+1)
- (calc_alpha_intercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) + calc_alpha_intercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col))) * grid.getConcentrations()(row, col)
+ calc_alpha_intercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col-1)) * grid.getConcentrations()(row, col-1));
}
concentrations_t1(0,0) = 0;
concentrations_t1(rowMax-1,0) = 0;
concentrations_t1(0,colMax-1) = 0;
concentrations_t1(rowMax-1,colMax-1) = 0;
grid.setConcentrations(concentrations_t1);
}
auto FTCS_closed(Grid &grid, Boundary &bc, double timestep) {
return 0;
}
auto FTCS(Grid &grid, Boundary &bc, double timestep) {
void FTCS(Grid &grid, Boundary &bc, double timestep) {
if (bc.getBoundaryConditionType() == BC_TYPE_CONSTANT) {
int test = FTCS_constant(grid, bc, timestep);
FTCS_constant(grid, bc, timestep);
} else if (bc.getBoundaryConditionType() == BC_TYPE_CLOSED) {
FTCS_closed(grid, bc, timestep);
}

2
src/Simulation.cpp
View File

@ -46,7 +46,7 @@ auto Simulation::getIterations() {
}
auto Simulation::run() {
void Simulation::run() {
if (approach == FTCS_APPROACH) {
for (int i = 0; i < iterations; i++) {
FTCS(grid, bc, timestep);