From 0a9b58e8ff42b61e8b108e639098cb31c85cbae0 Mon Sep 17 00:00:00 2001 From: philippun Date: Tue, 18 Jul 2023 15:27:17 +0200 Subject: [PATCH] implementing FTCS --- include/tug/Boundary.hpp | 7 ++++- include/tug/Grid.hpp | 20 +++++++++++-- src/Boundary.cpp | 2 +- src/FTCS.cpp | 62 ++++++++++++++++++++++++++++++++++++++-- src/Grid.cpp | 33 +++++++++++++++++++-- src/Simulation.cpp | 7 ++++- 6 files changed, 122 insertions(+), 9 deletions(-) diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index fb59f47..383a6b9 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -1,3 +1,6 @@ +#ifndef BOUNDARY_H_ +#define BOUNDARY_H_ + #include #include "Grid.hpp" @@ -31,7 +34,7 @@ class Boundary { * * @return auto */ - auto getBoundaryConditionType(); + BC_TYPE getBoundaryConditionType(); /** * @brief Set the Boundary Condition Value object @@ -55,3 +58,5 @@ class Boundary { BC_TYPE type; VectorXd left, right, top, bottom; }; + +#endif diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index 7736b73..0f828f3 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -32,7 +32,7 @@ class Grid { * * @return auto */ - auto getConcentrations(); + Matrix2d getConcentrations(); /** * @brief Set the Alpha object @@ -49,18 +49,34 @@ class Grid { */ void setAlpha(Matrix2d alpha_x, Matrix2d alpha_y); + Matrix2d getAlphaX(); + + Matrix2d getAlphaY(); + int getDim(); int getRow(); int getCol(); + void setDomain(int domain_col); + + void setDomain(int domain_row, int domain_col); + + double getDeltaCol(); + + double getDeltaRow(); + private: int dim; - int row; int col; + int row; + int domain_col; + int domain_row; + double delta_col; + double delta_row; Matrix2d concentrations; Matrix2d alpha_x; Matrix2d alpha_y; diff --git a/src/Boundary.cpp b/src/Boundary.cpp index 5b95f72..2c25848 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -23,7 +23,7 @@ Boundary::Boundary(Grid &grid, BC_TYPE type) : grid(grid) { } } -auto Boundary::getBoundaryConditionType() { +BC_TYPE Boundary::getBoundaryConditionType() { return this->type; } diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 94c2595..f7aca0d 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -1,5 +1,63 @@ #include -auto FTCS(Grid &grid, Boundary &bc, double timestep) { -} \ No newline at end of file +auto calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = false) { + if (useHarmonic) { + return 2 / ((1/alpha1) + (1/alpha2)); + } else { + return 0.5 * (alpha1 + alpha2); + } +} + +auto FTCS_constant(Grid &grid, Boundary &bc, double timestep) { + int rowMax = grid.getRow(); + int colMax = grid.getCol(); + double deltaRow = grid.getDeltaRow(); + double deltaCol = grid.getDeltaCol(); + + // Matrix with concentrations at time t+1 + // TODO profiler / only use 2 matrices + Matrix2d concentrations_t1 = Matrix2d(rowMax, colMax); + + // inner cells + 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 * ( + 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) + ) + - 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)) + + 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)) + * grid.getConcentrations()(row,col-1) + ); + } + } + + // boundary conditions + // left without corners / looping over rows + + return 0; +} + +auto FTCS_closed(Grid &grid, Boundary &bc, double timestep) { + return 0; +} + +auto FTCS(Grid &grid, Boundary &bc, double timestep) { + if (bc.getBoundaryConditionType() == BC_TYPE_CONSTANT) { + int test = FTCS_constant(grid, bc, timestep); + } else if (bc.getBoundaryConditionType() == BC_TYPE_CLOSED) { + FTCS_closed(grid, bc, timestep); + } +} diff --git a/src/Grid.cpp b/src/Grid.cpp index 65fe6cc..55ea4f9 100644 --- a/src/Grid.cpp +++ b/src/Grid.cpp @@ -22,7 +22,7 @@ void Grid::setConcentrations(Matrix2d concentrations) { this->concentrations = concentrations; } -auto Grid::getConcentrations() { +Matrix2d Grid::getConcentrations() { return this->concentrations; } @@ -35,6 +35,14 @@ void Grid::setAlpha(Matrix2d alpha_x, Matrix2d alpha_y) { this->alpha_y = alpha_y; } +Matrix2d Grid::getAlphaX() { + return this->alpha_x; +} + +Matrix2d Grid::getAlphaY() { + return this->alpha_y; +} + int Grid::getDim() { return dim; } @@ -45,4 +53,25 @@ int Grid::getRow() { int Grid::getCol() { return col; -} \ No newline at end of file +} + +void Grid::setDomain(int domain_col) { + this->domain_col = domain_col; + this->delta_col = this->domain_col/this->col; +} + +void Grid::setDomain(int domain_row, int domain_col) { + this->domain_row = domain_row; + this->domain_col = domain_col; + + this->domain_row = this->domain_row/this->row; + this->domain_col = this->domain_col/this->col; +} + +double Grid::getDeltaCol() { + return this->delta_col; +} + +double Grid::getDeltaRow() { + return this->delta_row; +} diff --git a/src/Simulation.cpp b/src/Simulation.cpp index b4cd947..0e44d32 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -5,6 +5,10 @@ using namespace std; +// auto FTCS(Grid &grid, Boundary &bc, double timestep) { + +// } + Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid), bc(bc) { //probably to DEBUG assignment of grid and bc this->grid = grid; @@ -41,6 +45,7 @@ auto Simulation::getIterations() { return this->iterations; } + auto Simulation::run() { if (approach == FTCS_APPROACH) { for (int i = 0; i < iterations; i++) { @@ -52,4 +57,4 @@ auto Simulation::run() { break; } } -} \ No newline at end of file +}