From 7cf4ad4772dc3e1d0306b89bf44b281a78885cd7 Mon Sep 17 00:00:00 2001 From: philippun Date: Wed, 26 Jul 2023 16:02:23 +0200 Subject: [PATCH 01/12] optimized example --- examples/FTCS_2D_proto_example.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index 5f159bb..27dd892 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -29,7 +29,7 @@ int main(int argc, char *argv[]) { // (optional) set the concentrations, e.g.: // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value // grid.setConcentrations(concentrations); - MatrixXd concentrations = MatrixXd::Constant(row, col,1); + MatrixXd concentrations = MatrixXd::Constant(row,col,0); concentrations(0,0) = 2000; grid.setConcentrations(concentrations); From 049fc319db9386f77884bd5a748e136b9786a870 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Wed, 26 Jul 2023 22:17:37 +0200 Subject: [PATCH 02/12] feat: add boundary conditons for closed cases in independent functions --- src/FTCS.cpp | 40 ++++++++++++++++++++++++++++++---------- 1 file changed, 30 insertions(+), 10 deletions(-) diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 58812a9..04194d0 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -62,9 +62,14 @@ double calcHorizontalChangeLeftBoundaryConstant(Grid grid, Boundary bc, int row, return result; } - -double calcHorizontalChangeLeftBoundaryClosed() { - return 0; +// TODO: REVIEW +double calcHorizontalChangeLeftBoundaryClosed(Grid grid, int row, int col) { + + double result = + calcAlphaIntercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) + * (grid.getConcentrations()(row, col+1) - grid.getConcentrations()(row, col)); + + return result; } @@ -83,9 +88,14 @@ double calcHorizontalChangeRightBoundaryConstant(Grid grid, Boundary bc, int row return result; } - -double calcHorizontalChangeRightBoundaryClosed() { - return 0; +// TODO: REVIEW +double calcHorizontalChangeRightBoundaryClosed(Grid grid, int row, int col) { + + double result = + - (calcAlphaIntercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col)) + * (grid.getConcentrations()(row, col) - grid.getConcentrations()(row, col-1))); + + return result; } @@ -104,9 +114,14 @@ double calcVerticalChangeTopBoundaryConstant(Grid grid, Boundary bc, int row, in return result; } +// TODO: REVIEW +double calcVerticalChangeTopBoundaryClosed(Grid grid, int row, int col) { + + double result = + calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getConcentrations()(row, col)) + * (grid.getConcentrations()(row+1, col) - grid.getConcentrations()(row, col)); -double calcVerticalChangeTopBoundaryClosed() { - return 0; + return result; } @@ -125,9 +140,14 @@ double calcVerticalChangeBottomBoundaryConstant(Grid grid, Boundary bc, int row, return result; } +// TODO: REVIEW +double calcVerticalChangeBottomBoundaryClosed(Grid grid, int row, int col) { -double calcVerticalChangeBottomBoundaryClosed() { - return 0; + double result = + - (calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) + * (grid.getConcentrations()(row, col) - grid.getConcentrations()(row-1, col))); + + return result; } From 8f48830bdab673cf86d9178ac805f135d7ab0531 Mon Sep 17 00:00:00 2001 From: philippun Date: Fri, 28 Jul 2023 10:46:05 +0200 Subject: [PATCH 03/12] reviewed closed case functions --- src/FTCS.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 04194d0..bb91844 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -62,7 +62,7 @@ double calcHorizontalChangeLeftBoundaryConstant(Grid grid, Boundary bc, int row, return result; } -// TODO: REVIEW + double calcHorizontalChangeLeftBoundaryClosed(Grid grid, int row, int col) { double result = @@ -88,7 +88,7 @@ double calcHorizontalChangeRightBoundaryConstant(Grid grid, Boundary bc, int row return result; } -// TODO: REVIEW + double calcHorizontalChangeRightBoundaryClosed(Grid grid, int row, int col) { double result = @@ -114,7 +114,7 @@ double calcVerticalChangeTopBoundaryConstant(Grid grid, Boundary bc, int row, in return result; } -// TODO: REVIEW + double calcVerticalChangeTopBoundaryClosed(Grid grid, int row, int col) { double result = @@ -140,7 +140,7 @@ double calcVerticalChangeBottomBoundaryConstant(Grid grid, Boundary bc, int row, return result; } -// TODO: REVIEW + double calcVerticalChangeBottomBoundaryClosed(Grid grid, int row, int col) { double result = From 88b00931f8b47415b12452703e8ad0230c13e626 Mon Sep 17 00:00:00 2001 From: philippun Date: Fri, 28 Jul 2023 10:48:07 +0200 Subject: [PATCH 04/12] test commit --- examples/FTCS_2D_proto_example.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index 27dd892..17bfaf8 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -30,7 +30,7 @@ int main(int argc, char *argv[]) { // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value // grid.setConcentrations(concentrations); MatrixXd concentrations = MatrixXd::Constant(row,col,0); - concentrations(0,0) = 2000; + concentrations(0,0) = 1999; grid.setConcentrations(concentrations); // (optional) set alphas of the grid, e.g.: From 3364042af714e7fe8df41a43988083bfa80f5ad7 Mon Sep 17 00:00:00 2001 From: philippun Date: Fri, 28 Jul 2023 10:55:28 +0200 Subject: [PATCH 05/12] unset easyprofiler as required --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a980a35..9cd2625 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ set(CMAKE_CXX_STANDARD 17) find_package(Eigen3 REQUIRED NO_MODULE) find_package(OpenMP) -find_package(easy_profiler REQUIRED) +find_package(easy_profiler) ## SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -mfma") option(TUG_USE_OPENMP "Compile with OpenMP support" ON) From c9b714b241f8b435f8207fe0f79d7e51a2377181 Mon Sep 17 00:00:00 2001 From: philippun Date: Fri, 28 Jul 2023 10:59:08 +0200 Subject: [PATCH 06/12] commented out easy profiler --- examples/FTCS_2D_proto_example.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index 27dd892..99502cc 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -8,16 +8,16 @@ #include #include -#define EASY_PROFILER_ENABLE ::profiler::setEnabled(true); +// #define EASY_PROFILER_ENABLE ::profiler::setEnabled(true); int main(int argc, char *argv[]) { - EASY_PROFILER_ENABLE; - profiler::startListen(); + // EASY_PROFILER_ENABLE; + // profiler::startListen(); // ************** // **** GRID **** // ************** - profiler::startListen(); + // profiler::startListen(); // create a grid with a 20 x 20 field int row = 20; int col = 20; @@ -77,9 +77,9 @@ int main(int argc, char *argv[]) { // run the simulation - EASY_BLOCK("SIMULATION") + // EASY_BLOCK("SIMULATION") simulation.run(); - EASY_END_BLOCK; - profiler::dumpBlocksToFile("test_profile.prof"); - profiler::stopListen(); + // EASY_END_BLOCK; + // profiler::dumpBlocksToFile("test_profile.prof"); + // profiler::stopListen(); } \ No newline at end of file From 605000de5e1dfc15ec27e6b5171038839acdc957 Mon Sep 17 00:00:00 2001 From: philippun Date: Fri, 28 Jul 2023 11:00:21 +0200 Subject: [PATCH 07/12] commented out another easy profiler include --- examples/FTCS_2D_proto_example.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index 99502cc..ac27d25 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -7,7 +7,7 @@ */ #include -#include +// #include // #define EASY_PROFILER_ENABLE ::profiler::setEnabled(true); From 16d3663909f7141eb6997725612f583e7abd04dd Mon Sep 17 00:00:00 2001 From: philippun Date: Fri, 28 Jul 2023 11:02:56 +0200 Subject: [PATCH 08/12] commented out the probably last easy profiler reference --- examples/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 47042c2..b7c0d71 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -9,4 +9,4 @@ target_link_libraries(second_example tug) target_link_libraries(boundary_example1D tug) target_link_libraries(FTCS_2D_proto_example tug) target_link_libraries(FTCS_1D_proto_example tug) -target_link_libraries(FTCS_2D_proto_example easy_profiler) \ No newline at end of file +# target_link_libraries(FTCS_2D_proto_example easy_profiler) \ No newline at end of file From 4680e9823fa5302eef6c4932e967d544a99c7a0a Mon Sep 17 00:00:00 2001 From: philippun Date: Fri, 28 Jul 2023 12:57:09 +0200 Subject: [PATCH 09/12] feat: reimplented Boundary to incorporate different types for each boundary element --- include/tug/Boundary.hpp | 82 ++++++++++++----------------- src/Boundary.cpp | 109 ++++++++++++++++++++------------------- src/FTCS.cpp | 13 +++++ 3 files changed, 104 insertions(+), 100 deletions(-) diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index df3ad0f..c4ad577 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -2,6 +2,7 @@ #define BOUNDARY_H_ #include +#include #include "Grid.hpp" using namespace std; @@ -19,26 +20,21 @@ enum BC_SIDE { BC_SIDE_BOTTOM }; -/******************************************* -class WallElement { +class BoundaryElement { public: - WallElement() { - this->type = BC_TYPE_CLOSED; - this->value = 0; - } + // bc type closed + BoundaryElement(); - WallElement(double value) { - this->type = BC_TYPE_CONSTANT; - this->value = value; - } + // bc type constant + BoundaryElement(double value); - BC_TYPE getType() { - return this->type; - } + void setType(BC_TYPE type); - double getValue() { - return this->value; - } + void setValue(double value); + + BC_TYPE getType(); + + double getValue(); private: BC_TYPE type; @@ -47,16 +43,22 @@ class WallElement { class BoundaryWall { public: - BoundaryWall(int length) { - // create array with length many wall elements - } + BoundaryWall(int length); + + void setWall(BC_TYPE type, double value = NAN); + + vector getWall(); + + void setBoundaryElement(int index, BC_TYPE type, double value = NAN); + + BoundaryElement getBoundaryElement(); private: BC_SIDE side; int length; - vector wall; + vector wall; + }; -***********************/ class Boundary { public: @@ -67,40 +69,24 @@ class Boundary { * @param grid * @param type */ - Boundary(Grid grid, BC_TYPE type); + Boundary(Grid grid); - /** - * @brief Get the Boundary Condition Type object - * - * @return auto - */ - BC_TYPE getBoundaryConditionType(); + void setBoundarySideClosed(BC_SIDE side); - /** - * @brief Set the Boundary Condition Value object - * - * @param side - * @param values - */ - void setBoundaryConditionValue(BC_SIDE side, VectorXd values); + void setBoundarySideConstant(BC_SIDE side, double value); - /** - * @brief Get the Boundary Condition Value object - * - * @param side - * @return auto - */ - VectorXd getBoundaryConditionValue(BC_SIDE side); + void setBoundaryElementClosed(BC_SIDE side, int index); + void setBoundaryElementConstant(BC_SIDE side, int index, double value); + + vector getBoundarySide(BC_SIDE side); + + BoundaryElement getBoundaryElement(BC_SIDE side, int index); private: Grid grid; - - // need a way to save the bc type and value for each single 'boundary cell' - // perhaps an array for each side with structs containing the bc type as well as a value - // or another object that contains one boundary side - BC_TYPE type; - VectorXd left, right, top, bottom; + + vector> boundaries; }; #endif diff --git a/src/Boundary.cpp b/src/Boundary.cpp index 516c422..3c3ac04 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -1,68 +1,73 @@ +#include "tug/BoundaryCondition.hpp" #include +#include #include #include using namespace std; -Boundary::Boundary(Grid grid, BC_TYPE type) : grid(grid) { - //probably to DEBUG assignment grid +BoundaryElement::BoundaryElement() { + this->type = BC_TYPE_CLOSED; + this->value = NAN; +} + +BoundaryElement::BoundaryElement(double value) { + this->type = BC_TYPE_CONSTANT; + this->value = value; +} + +void BoundaryElement::setType(BC_TYPE type) { this->type = type; - - if (type == BC_TYPE_CONSTANT) { - if (grid.getDim() == 1) { - this->left = VectorXd::Constant(1, 1); - this->right = VectorXd::Constant(1, 1); - } else if (grid.getDim() == 2) { - this->left = VectorXd::Constant(grid.getRow(), 1); - this->right = VectorXd::Constant(grid.getRow(), 1); - this->top = VectorXd::Constant(grid.getCol(), 1); - this->bottom = VectorXd::Constant(grid.getCol(), 1); - - } else { - throw invalid_argument("Dimension must be 1 or 2!"); - } - } } -BC_TYPE Boundary::getBoundaryConditionType() { - return this->type; +void BoundaryElement::setValue(double value) { + this->value = value; } -void Boundary::setBoundaryConditionValue(BC_SIDE side, VectorXd values) { - if (type != BC_TYPE_CONSTANT) { - // TODO check if correct way for handling warning - cerr << "Values will not be used, wrong BC_TYPE!"; - } +BC_TYPE BoundaryElement::getType() { + return this->type; +} - switch (side) { - case BC_SIDE_LEFT: - this->left = values; - break; - case BC_SIDE_RIGHT: - this->right = values; - break; - case BC_SIDE_TOP: - this->top = values; - break; - case BC_SIDE_BOTTOM: - this->bottom = values; - break; - default: - throw invalid_argument("Invalid side given!"); +double BoundaryElement::getValue() { + return this->value; +} + +Boundary::Boundary(Grid grid) : grid(grid) { + //probably to DEBUG assignment grid + + + if (grid.getDim() == 1) { + this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); + } else if (grid.getDim() == 2) { + this->boundaries[BC_SIDE_LEFT] = vector(grid.getRow(), BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT] = vector(grid.getRow(), BoundaryElement()); + this->boundaries[BC_SIDE_TOP] = vector(grid.getCol(), BoundaryElement()); + this->boundaries[BC_SIDE_BOTTOM] = vector(grid.getCol(), BoundaryElement()); } } -VectorXd Boundary::getBoundaryConditionValue(BC_SIDE side) { - switch (side) { - case BC_SIDE_LEFT: - return this->left; - case BC_SIDE_RIGHT: - return this->right; - case BC_SIDE_TOP: - return this->top; - case BC_SIDE_BOTTOM: - return this->bottom; - default: - throw invalid_argument("Invalid side given!"); - } +void Boundary::setBoundarySideClosed(BC_SIDE side) { + this->boundaries[side] = vector(grid.getRow(), BoundaryElement()); +} + +void Boundary::setBoundarySideConstant(BC_SIDE side, double value) { + this->boundaries[side] = vector(grid.getRow(), BoundaryElement(value)); +} + +void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) { + this->boundaries[side][index].setType(BC_TYPE_CLOSED); +} + +void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, double value) { + this->boundaries[side][index].setType(BC_TYPE_CONSTANT); + this->boundaries[side][index].setValue(value); +} + +vector Boundary::getBoundarySide(BC_SIDE side) { + return this->boundaries[side]; +} + +BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) { + return this->boundaries[side][index]; } \ No newline at end of file diff --git a/src/FTCS.cpp b/src/FTCS.cpp index bb91844..1b24e5a 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -224,6 +224,19 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) { // 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); + // if (bc.getBoundaryConditionType(BC_SIDE_LEFT, row)) { + // concentrations_t1(row, col) += timestep / (deltaCol*deltaCol) + // + calcHorizontalChangeLeftBoundaryClosed(grid, row, col); + // } else { + // concentrations_t1(row, col) += timestep / (deltaCol*deltaCol) + // + calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col); + // } + // concentrations_t1(row, col) += timestep / (deltaRow*deltaRow) + // * ( + // calcVerticalChange(grid, row, col) + // ); + concentrations_t1(row, col) = grid.getConcentrations()(row,col) + timestep / (deltaCol*deltaCol) * ( From 0a0e16bb56a3a39d67a3b808ff40f4450286968f Mon Sep 17 00:00:00 2001 From: philippun Date: Fri, 28 Jul 2023 13:02:45 +0200 Subject: [PATCH 10/12] feat: added getters to Boundary --- include/tug/Boundary.hpp | 5 +++++ src/Boundary.cpp | 11 ++++++++++- 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index c4ad577..bf1d6a9 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -83,6 +83,10 @@ class Boundary { BoundaryElement getBoundaryElement(BC_SIDE side, int index); + BC_TYPE getBoundaryElementType(BC_SIDE side, int index); + + double getBoundaryElementValue(BC_SIDE side, int index); + private: Grid grid; @@ -90,3 +94,4 @@ class Boundary { }; #endif + diff --git a/src/Boundary.cpp b/src/Boundary.cpp index 3c3ac04..9a1012a 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -70,4 +70,13 @@ vector Boundary::getBoundarySide(BC_SIDE side) { BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) { return this->boundaries[side][index]; -} \ No newline at end of file +} + +BC_TYPE Boundary::getBoundaryElementType(BC_SIDE side, int index) { + return this->boundaries[side][index].getType(); +} + +double Boundary::getBoundaryElementValue(BC_SIDE side, int index) { + return this->boundaries[side][index].getValue(); +} + From 0872639c540f9f995dbcf4ee8f2c861567e6d3c6 Mon Sep 17 00:00:00 2001 From: philippun Date: Fri, 28 Jul 2023 15:54:03 +0200 Subject: [PATCH 11/12] feat: implementation of dynamic boundary conditions in FTCS --- examples/FTCS_1D_proto_example.cpp | 5 +- examples/FTCS_2D_proto_example.cpp | 9 +++- src/Boundary.cpp | 5 +- src/FTCS.cpp | 82 +++++++++++++++++++++++------- 4 files changed, 79 insertions(+), 22 deletions(-) diff --git a/examples/FTCS_1D_proto_example.cpp b/examples/FTCS_1D_proto_example.cpp index c0936c4..3ac0df6 100644 --- a/examples/FTCS_1D_proto_example.cpp +++ b/examples/FTCS_1D_proto_example.cpp @@ -1,3 +1,4 @@ +#include "tug/Boundary.hpp" #include 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); // ************************ diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index 7171f62..deb778d 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -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); diff --git a/src/Boundary.cpp b/src/Boundary.cpp index 9a1012a..87752ea 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -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>(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>(4); + this->boundaries[BC_SIDE_LEFT] = vector(grid.getRow(), BoundaryElement()); this->boundaries[BC_SIDE_RIGHT] = vector(grid.getRow(), BoundaryElement()); this->boundaries[BC_SIDE_TOP] = vector(grid.getCol(), BoundaryElement()); diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 1b24e5a..f00891a 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -1,9 +1,11 @@ +#include "TugUtils.hpp" #include #include #include 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) ) ; From ffa48347b81f859364290983d98f7f467fff0bcd Mon Sep 17 00:00:00 2001 From: philippun Date: Fri, 28 Jul 2023 15:55:59 +0200 Subject: [PATCH 12/12] removed some old commentary --- src/FTCS.cpp | 33 +-------------------------------- 1 file changed, 1 insertion(+), 32 deletions(-) diff --git a/src/FTCS.cpp b/src/FTCS.cpp index f00891a..277468f 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -270,19 +270,6 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) { // 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); - // if (bc.getBoundaryConditionType(BC_SIDE_LEFT, row)) { - // concentrations_t1(row, col) += timestep / (deltaCol*deltaCol) - // + calcHorizontalChangeLeftBoundaryClosed(grid, row, col); - // } else { - // concentrations_t1(row, col) += timestep / (deltaCol*deltaCol) - // + calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col); - // } - // concentrations_t1(row, col) += timestep / (deltaRow*deltaRow) - // * ( - // calcVerticalChange(grid, row, col) - // ); - concentrations_t1(row, col) = grid.getConcentrations()(row,col) + timestep / (deltaCol*deltaCol) * ( @@ -403,29 +390,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) { MatrixXd FTCS(Grid grid, Boundary bc, double timestep) { - // inner cells - // TODO only the boundary cells are different in constant and closed case - - // if 1D: - // do inner cells - // do left boundary according to bc type - // do right boundary according to bc type - // if 2D: - // do inner cells - // do left boundaries according to bc type - // do right boundaries according to bc type - // ... - + if (grid.getDim() == 1) { return FTCS_1D(grid, bc, timestep); } else { return FTCS_2D(grid, bc, timestep); } - // checking the boundary condition type first does not work - // if the boundary condition types change dynamically for a grid - // meaning: - // - boundary condition type needs to be checked for every single boundary cell - // -> this check is last in order - // - }