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) 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 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 7673bc2..dd7cc37 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -7,17 +7,17 @@ */ #include -#include -#define EASY_PROFILER_ENABLE ::profiler::setEnabled(true); +// #include +// #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; @@ -29,8 +29,8 @@ 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); - concentrations(0,0) = 2000; + MatrixXd concentrations = MatrixXd::Constant(row,col,0); + concentrations(0,0) = 1999; grid.setConcentrations(concentrations); @@ -45,7 +45,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 @@ -69,7 +74,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); @@ -78,9 +83,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 diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index df3ad0f..bf1d6a9 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,29 @@ 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); + + BC_TYPE getBoundaryElementType(BC_SIDE side, int index); + + double getBoundaryElementValue(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..87752ea 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -1,68 +1,85 @@ +#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 = 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()); + 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!"); - } -} \ No newline at end of file +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]; +} + +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(); +} + diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 58812a9..277468f 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,21 +59,37 @@ 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; } -double calcHorizontalChangeLeftBoundaryClosed() { - return 0; +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; +} + + +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) @@ -84,8 +102,24 @@ double calcHorizontalChangeRightBoundaryConstant(Grid grid, Boundary bc, int row } -double calcHorizontalChangeRightBoundaryClosed() { - return 0; +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; +} + + +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!"); + } } @@ -99,21 +133,37 @@ 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; } -double calcVerticalChangeTopBoundaryClosed() { - return 0; +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)); + + return result; +} + + +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) @@ -126,8 +176,24 @@ double calcVerticalChangeBottomBoundaryConstant(Grid grid, Boundary bc, int row, } -double calcVerticalChangeBottomBoundaryClosed() { - return 0; +double calcVerticalChangeBottomBoundaryClosed(Grid grid, int row, int col) { + + double result = + - (calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) + * (grid.getConcentrations()(row, col) - grid.getConcentrations()(row-1, col))); + + return result; +} + + +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!"); + } } @@ -155,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) ) ; @@ -165,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) ) ; @@ -207,7 +273,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) * ( @@ -222,7 +288,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) * ( @@ -238,7 +304,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) * ( @@ -253,7 +319,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) * ( @@ -268,11 +334,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) ) ; @@ -282,11 +348,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) ) ; @@ -296,11 +362,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) ) ; @@ -310,11 +376,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) ) ; @@ -324,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 - // - }