From d457c2b9a7dc42c772f3127696eff0b6e6e81fb9 Mon Sep 17 00:00:00 2001 From: philippun Date: Wed, 19 Jul 2023 11:19:00 +0200 Subject: [PATCH] fixed Eigen implementation bugs --- examples/FTCS_2D_proto_example.cpp | 1 + include/tug/Boundary.hpp | 2 +- include/tug/Simulation.hpp | 2 +- src/Boundary.cpp | 2 +- src/FTCS.cpp | 131 +++++++++++++++-------------- src/Grid.cpp | 6 +- src/Simulation.cpp | 11 ++- 7 files changed, 85 insertions(+), 70 deletions(-) diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index fba5876..9586454 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -8,6 +8,7 @@ int main(int argc, char *argv[]) { Boundary bc = Boundary(grid, BC_TYPE_CONSTANT); Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); + simulation.setIterations(2); simulation.run(); diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index c49d456..fc7543a 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -27,7 +27,7 @@ class Boundary { * @param grid * @param type */ - Boundary(Grid &grid, BC_TYPE type); + Boundary(Grid grid, BC_TYPE type); /** * @brief Get the Boundary Condition Type object diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 8d0d809..8e27308 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -21,7 +21,7 @@ class Simulation { * @param bc * @param aproach */ - Simulation(Grid &grid, Boundary &bc, APPROACH approach); + Simulation(Grid grid, Boundary bc, APPROACH approach); /** * @brief diff --git a/src/Boundary.cpp b/src/Boundary.cpp index 311a439..b464f0d 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -4,7 +4,7 @@ using namespace std; -Boundary::Boundary(Grid &grid, BC_TYPE type) : grid(grid) { +Boundary::Boundary(Grid grid, BC_TYPE type) : grid(grid) { //probably to DEBUG assignment grid this->type = type; diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 695cd1d..65cce82 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -1,6 +1,9 @@ +#include #include #include +using namespace std; + auto calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = false) { if (useHarmonic) { return 2 / ((1/alpha1) + (1/alpha2)); @@ -9,7 +12,7 @@ auto calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = false } } -auto FTCS_constant(Grid &grid, Boundary &bc, double timestep) { +MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) { int rowMax = grid.getRow(); int colMax = grid.getCol(); double deltaRow = grid.getDeltaRow(); @@ -17,9 +20,13 @@ auto FTCS_constant(Grid &grid, Boundary &bc, double timestep) { // Matrix with concentrations at time t+1 // TODO profiler / only use 2 matrices - MatrixXd concentrations_t1 = MatrixXd(rowMax, colMax); + MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 1); // inner cells + cout << "Concentration 5,5: " << grid.getConcentrations()(5,5) << endl; + cout << "Alpha Y 5,5: " << grid.getAlphaY()(5,5) << endl; + cout << "calc alpha Y 5,5; 5,6: " << calc_alpha_intercell(grid.getAlphaY()(5,5), grid.getAlphaY()(5,6)) << endl; + cout << "t1 Concentrations 5,5: " << concentrations_t1(5,5) << endl; for (int row = 1; row < rowMax-1; row++) { for (int col = 1; col < colMax-1; col++) { concentrations_t1(row, col) = grid.getConcentrations()(row, col) @@ -32,7 +39,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)) @@ -47,69 +54,69 @@ 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)); - } + // 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)) + // + 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)); + // } // 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)); - } + // 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) + // - (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; coldomain_col = domain_col; - this->delta_col = this->domain_col/this->col; + this->delta_col = double(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; + this->delta_row = double(this->domain_row)/this->row; + this->delta_col = double(this->domain_col)/this->col; } double Grid::getDeltaCol() { diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 92a042a..44d1736 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -9,7 +9,7 @@ using namespace std; // } -Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid), bc(bc) { +Simulation::Simulation(Grid grid, Boundary bc, APPROACH approach) : grid(grid), bc(bc) { //probably to DEBUG assignment of grid and bc this->grid = grid; this->approach = approach; @@ -48,9 +48,16 @@ auto Simulation::getIterations() { void Simulation::run() { if (approach == FTCS_APPROACH) { + cout << grid.getConcentrations() << endl; + cout << endl; for (int i = 0; i < iterations; i++) { - FTCS(grid, bc, timestep); + grid.setConcentrations(FTCS(grid, bc, timestep)); + if (i == 10) { + cout << grid.getConcentrations() << endl; + cout << endl; + } } + cout << grid.getConcentrations() << endl; } else if (approach == BTCS_APPROACH) { for (int i = 0; i < iterations; i++) { //TODO