From ba627b66245e37ba0681a9412c5f102eebcfe426 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 15 Sep 2023 10:39:51 +0200 Subject: [PATCH] feat: rewrite library as template library --- examples/BTCS_1D_proto_example.cpp | 2 +- examples/BTCS_2D_proto_example.cpp | 2 +- examples/CRNI_2D_proto_example.cpp | 2 +- examples/FTCS_1D_proto_example.cpp | 2 +- examples/FTCS_2D_proto_closed_mdl.cpp | 2 +- examples/FTCS_2D_proto_example.cpp | 2 +- examples/FTCS_2D_proto_example_mdl.cpp | 2 +- examples/profiling_openmp.cpp | 2 +- examples/profiling_speedup.cpp | 2 +- examples/reference-FTCS_2D_closed.cpp | 2 +- include/tug/Boundary.hpp | 172 +++++++++++++++--- include/tug/Grid.hpp | 236 ++++++++++++++++++++----- include/tug/Simulation.hpp | 145 +++++++++++++-- naaice/BTCS_2D_NAAICE.cpp | 8 +- src/BTCS.cpp | 132 ++++++++------ src/Boundary.cpp | 149 ---------------- src/CMakeLists.txt | 2 +- src/FTCS.cpp | 105 ++++++----- src/Grid.cpp | 167 ----------------- src/Schemes.hpp | 9 +- src/Simulation.cpp | 166 +++-------------- test/testBoundary.cpp | 11 +- test/testGrid.cpp | 20 +-- test/testSimulation.cpp | 8 +- 24 files changed, 664 insertions(+), 686 deletions(-) delete mode 100644 src/Boundary.cpp delete mode 100644 src/Grid.cpp diff --git a/examples/BTCS_1D_proto_example.cpp b/examples/BTCS_1D_proto_example.cpp index 1a53738..4bfd1b4 100644 --- a/examples/BTCS_1D_proto_example.cpp +++ b/examples/BTCS_1D_proto_example.cpp @@ -10,7 +10,7 @@ int main(int argc, char *argv[]) { // create a linear grid with 20 cells int cells = 20; - Grid grid = Grid(cells); + Grid64 grid(cells); MatrixXd concentrations = MatrixXd::Constant(1, 20, 0); concentrations(0, 0) = 2000; diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index 9a103cb..f14e0f5 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -13,7 +13,7 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int row = 40; int col = 50; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); diff --git a/examples/CRNI_2D_proto_example.cpp b/examples/CRNI_2D_proto_example.cpp index 39262e7..d605038 100644 --- a/examples/CRNI_2D_proto_example.cpp +++ b/examples/CRNI_2D_proto_example.cpp @@ -6,7 +6,7 @@ using namespace Eigen; int main(int argc, char *argv[]) { int row = 20; int col = 20; - Grid grid(row, col); + Grid64 grid(row, col); MatrixXd concentrations = MatrixXd::Constant(row, col, 0); concentrations(10, 10) = 2000; diff --git a/examples/FTCS_1D_proto_example.cpp b/examples/FTCS_1D_proto_example.cpp index e97f0be..74b6f3e 100644 --- a/examples/FTCS_1D_proto_example.cpp +++ b/examples/FTCS_1D_proto_example.cpp @@ -11,7 +11,7 @@ int main(int argc, char *argv[]) { // create a linear grid with 20 cells int cells = 20; - Grid grid = Grid(cells); + Grid64 grid(cells); MatrixXd concentrations = MatrixXd::Constant(1, 20, 20); grid.setConcentrations(concentrations); diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index 845d36f..e17f92b 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -30,7 +30,7 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int n2 = row / 2 - 1; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index 5f5e520..89f821d 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -23,7 +23,7 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int row = 20; int col = 20; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp index 2c63c0c..d77710d 100644 --- a/examples/FTCS_2D_proto_example_mdl.cpp +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -21,7 +21,7 @@ int main(int argc, char *argv[]) { int row = 64; int col = 64; int n2 = row / 2 - 1; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); diff --git a/examples/profiling_openmp.cpp b/examples/profiling_openmp.cpp index 5dc13ef..c781900 100644 --- a/examples/profiling_openmp.cpp +++ b/examples/profiling_openmp.cpp @@ -29,7 +29,7 @@ int main(int argc, char *argv[]) { // myfile << "Iterations: " << iterations[j] << endl; for (int k = 0; k < repetition; k++) { cout << "Wiederholung: " << k << endl; - Grid grid = Grid(n[i], n[i]); + Grid64 grid(n[i], n[i]); grid.setDomain(1, 1); MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); diff --git a/examples/profiling_speedup.cpp b/examples/profiling_speedup.cpp index 5dc13ef..c781900 100644 --- a/examples/profiling_speedup.cpp +++ b/examples/profiling_speedup.cpp @@ -29,7 +29,7 @@ int main(int argc, char *argv[]) { // myfile << "Iterations: " << iterations[j] << endl; for (int k = 0; k < repetition; k++) { cout << "Wiederholung: " << k << endl; - Grid grid = Grid(n[i], n[i]); + Grid64 grid(n[i], n[i]); grid.setDomain(1, 1); MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); diff --git a/examples/reference-FTCS_2D_closed.cpp b/examples/reference-FTCS_2D_closed.cpp index 340a6fd..fa16063 100644 --- a/examples/reference-FTCS_2D_closed.cpp +++ b/examples/reference-FTCS_2D_closed.cpp @@ -12,7 +12,7 @@ int main(int argc, char *argv[]) { int domain_col = 10; // Grid - Grid grid = Grid(row, col); + Grid64 grid(row, col); grid.setDomain(domain_row, domain_col); MatrixXd concentrations = MatrixXd::Constant(row, col, 0); diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index b803e75..b51df24 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -9,6 +9,8 @@ #include "Grid.hpp" #include +#include +#include /** * @brief Enum defining the two implemented boundary conditions. @@ -27,7 +29,7 @@ enum BC_SIDE { BC_SIDE_LEFT, BC_SIDE_RIGHT, BC_SIDE_TOP, BC_SIDE_BOTTOM }; * These can be flexibly used and combined later in other classes. * The class serves as an auxiliary class for structuring the Boundary class. */ -class BoundaryElement { +template class BoundaryElement { public: /** * @brief Construct a new Boundary Element object for the closed case. @@ -35,7 +37,7 @@ public: * BC_TYPE_CLOSED, where the value takes -1 and does not hold any * physical meaning. */ - BoundaryElement(); + BoundaryElement(){}; /** * @brief Construct a new Boundary Element object for the constant case. @@ -45,7 +47,7 @@ public: * @param value Value of the constant concentration to be assumed at the * corresponding boundary element. */ - BoundaryElement(double value); + BoundaryElement(T _value) : value(_value), type(BC_TYPE_CONSTANT) {} /** * @brief Allows changing the boundary type of a corresponding @@ -54,7 +56,7 @@ public: * @param type Type of boundary condition. Either BC_TYPE_CONSTANT or BC_TYPE_CLOSED. */ - void setType(BC_TYPE type); + void setType(BC_TYPE type) { this->type = type; }; /** * @brief Sets the value of a boundary condition for the constant case. @@ -62,7 +64,17 @@ public: * @param value Concentration to be considered constant for the * corresponding boundary element. */ - void setValue(double value); + void setValue(double value) { + if (value < 0) { + throw std::invalid_argument("No negative concentration allowed."); + } + if (type == BC_TYPE_CLOSED) { + throw std::invalid_argument( + "No constant boundary concentrations can be set for closed " + "boundaries. Please change type first."); + } + this->value = value; + } /** * @brief Return the type of the boundary condition, i.e. whether the @@ -71,25 +83,25 @@ public: * @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or BC_TYPE_CONSTANT. */ - BC_TYPE getType(); + BC_TYPE getType() const { return this->type; } /** * @brief Return the concentration value for the constant boundary condition. * * @return double Value of the concentration. */ - double getValue(); + T getValue() const { return this->value; } private: BC_TYPE type{BC_TYPE_CLOSED}; - double value{-1}; + T value{-1}; }; /** * This class implements the functionality and management of the boundary * conditions in the grid to be simulated. */ -class Boundary { +template class Boundary { public: /** * @brief Creates a boundary object based on the passed grid object and @@ -98,7 +110,27 @@ public: * @param grid Grid object on the basis of which the simulation takes place * and from which the dimensions (in 2D case) are taken. */ - Boundary(Grid grid); + Boundary(const Grid &grid) + : dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) { + if (this->dim == 1) { + this->boundaries = std::vector>>( + 2); // in 1D only left and right boundary + + this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); + } else if (this->dim == 2) { + this->boundaries = std::vector>>(4); + + this->boundaries[BC_SIDE_LEFT] = + std::vector>(this->rows, BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT] = + std::vector>(this->rows, BoundaryElement()); + this->boundaries[BC_SIDE_TOP] = + std::vector>(this->cols, BoundaryElement()); + this->boundaries[BC_SIDE_BOTTOM] = + std::vector>(this->cols, BoundaryElement()); + } + } /** * @brief Sets all elements of the specified boundary side to the boundary @@ -106,7 +138,21 @@ public: * * @param side Side to be set to closed, e.g. BC_SIDE_LEFT. */ - void setBoundarySideClosed(BC_SIDE side); + void setBoundarySideClosed(BC_SIDE side) { + if (this->dim == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw std::invalid_argument( + "For the one-dimensional case, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } + + const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; + const int n = is_vertical ? this->rows : this->cols; + + this->boundaries[side] = + std::vector>(n, BoundaryElement()); + } /** * @brief Sets all elements of the specified boundary side to the boundary @@ -117,7 +163,21 @@ public: * @param value Concentration to be set for all elements of the specified * page. */ - void setBoundarySideConstant(BC_SIDE side, double value); + void setBoundarySideConstant(BC_SIDE side, double value) { + if (this->dim == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw std::invalid_argument( + "For the one-dimensional case, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } + + const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; + const int n = is_vertical ? this->rows : this->cols; + + this->boundaries[side] = + std::vector>(n, BoundaryElement(value)); + } /** * @brief Specifically sets the boundary element of the specified side @@ -128,7 +188,13 @@ public: * boundary side. Must index an element of the corresponding * side. */ - void setBoundaryElementClosed(BC_SIDE side, int index); + void setBoundaryElemenClosed(BC_SIDE side, int index) { + // tests whether the index really points to an element of the boundary side. + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument("Index is selected either too large or too small."); + } + this->boundaries[side][index].setType(BC_TYPE_CLOSED); + } /** * @brief Specifically sets the boundary element of the specified side @@ -142,7 +208,14 @@ public: * @param value Concentration value to which the boundary element should be set. */ - void setBoundaryElementConstant(BC_SIDE side, int index, double value); + void setBoundaryElementConstant(BC_SIDE side, int index, double value) { + // tests whether the index really points to an element of the boundary side. + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument("Index is selected either too large or too small."); + } + this->boundaries[side][index].setType(BC_TYPE_CONSTANT); + this->boundaries[side][index].setValue(value); + } /** * @brief Returns the boundary condition of a specified side as a vector @@ -150,19 +223,41 @@ public: * * @param side Boundary side from which the boundary conditions are to be * returned. - * @return vector Contains the boundary conditions as - * BoundaryElement objects. + * @return vector> Contains the boundary conditions as + * BoundaryElement objects. */ - const std::vector getBoundarySide(BC_SIDE side); + const std::vector> &getBoundarySide(BC_SIDE side) const { + if (this->dim == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw std::invalid_argument( + "For the one-dimensional trap, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } + return this->boundaries[side]; + } /** * @brief Get thes Boundary Side Values as a vector. Value is -1 in case some specific boundary is closed. * * @param side Boundary side for which the values are to be returned. - * @return VectorXd Vector with values as doubles. + * @return VectorX Vector with values as T. */ - Eigen::VectorXd getBoundarySideValues(BC_SIDE side); + Eigen::VectorX getBoundarySideValues(BC_SIDE side) const { + const std::size_t length = boundaries[side].size(); + Eigen::VectorX values(length); + + for (int i = 0; i < length; i++) { + if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { + values(i) = -1; + continue; + } + values(i) = getBoundaryElementValue(side, i); + } + + return values; + } /** * @brief Returns the boundary condition of a specified element on a given @@ -172,9 +267,15 @@ public: * @param index Index of the boundary element on the corresponding * boundary side. Must index an element of the corresponding * side. - * @return BoundaryElement Boundary condition as a BoundaryElement object. + * @return BoundaryElement Boundary condition as a BoundaryElement + * object. */ - BoundaryElement getBoundaryElement(BC_SIDE side, int index); + BoundaryElement getBoundaryElement(BC_SIDE side, int index) const { + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument("Index is selected either too large or too small."); + } + return this->boundaries[side][index]; + } /** * @brief Returns the type of a boundary condition, i.e. either BC_TYPE_CLOSED @@ -186,25 +287,42 @@ public: side. * @return BC_TYPE Boundary Type of the corresponding boundary condition. */ - BC_TYPE getBoundaryElementType(BC_SIDE side, int index); + BC_TYPE getBoundaryElementType(BC_SIDE side, int index) const { + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument("Index is selected either too large or too small."); + } + return this->boundaries[side][index].getType(); + } /** * @brief Returns the concentration value of a corresponding - * BoundaryElement object if it is a constant boundary condition. + * BoundaryElement object if it is a constant boundary condition. * * @param side Boundary side in which the boundary condition value is * located. * @param index Index of the boundary element on the corresponding * boundary side. Must index an element of the corresponding * side. - * @return double Concentration of the corresponding BoundaryElement object. + * @return double Concentration of the corresponding BoundaryElement + * object. */ - double getBoundaryElementValue(BC_SIDE side, int index); + T getBoundaryElementValue(BC_SIDE side, int index) const { + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument("Index is selected either too large or too small."); + } + if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) { + throw std::invalid_argument( + "A value can only be output if it is a constant boundary condition."); + } + return this->boundaries[side][index].getValue(); + } private: - Grid grid; // Boundary is directly dependent on the dimensions of a predefined + const std::uint8_t dim; + const std::uint32_t cols; + const std::uint32_t rows; - std::vector> + std::vector>> boundaries; // Vector with Boundary Element information }; diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index 6a6f973..3d3638d 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -10,8 +10,9 @@ #include #include +#include -class Grid { +template class Grid { public: /** * @brief Constructs a new 1D-Grid object of a given length, which holds a @@ -21,7 +22,18 @@ public: * * @param length Length of the 1D-Grid. Must be greater than 3. */ - Grid(int length); + Grid(int length) : col(length), domainCol(length) { + if (length <= 3) { + throw std::invalid_argument( + "Given grid length too small. Must be greater than 3."); + } + + this->dim = 1; + this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 + + this->concentrations = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); + this->alphaX = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); + } /** * @brief Constructs a new 2D-Grid object of given dimensions, which holds a @@ -34,146 +46,280 @@ public: * @param row Length of the 2D-Grid in y-direction. Must be greater than 3. * @param col Length of the 2D-Grid in x-direction. Must be greater than 3. */ - Grid(int row, int col); + Grid(int _row, int _col) + : row(_row), col(_col), domainRow(_row), domainCol(_col) { + if (row <= 3 || col <= 3) { + throw std::invalid_argument( + "Given grid dimensions too small. Must each be greater than 3."); + } + + this->dim = 2; + this->deltaRow = double(this->domainRow) / double(this->row); // -> 1 + this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 + + this->concentrations = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); + this->alphaX = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); + this->alphaY = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); + } /** * @brief Sets the concentrations matrix for a 1D or 2D-Grid. * - * @param concentrations An Eigen3 MatrixXd holding the concentrations. Matrix - * must have correct dimensions as defined in row and col. (Or length, in 1D - * case). + * @param concentrations An Eigen3 MatrixX holding the concentrations. + * Matrix must have correct dimensions as defined in row and col. (Or length, + * in 1D case). */ - void setConcentrations(const Eigen::MatrixXd &concentrations); + void setConcentrations(const Eigen::MatrixX &concentrations) { + if (concentrations.rows() != this->row || + concentrations.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of concentrations mismatch with Grid dimensions!"); + } + + this->concentrations = concentrations; + } /** * @brief Gets the concentrations matrix for a Grid. * - * @return MatrixXd An Eigen3 matrix holding the concentrations and having the - * same dimensions as the grid. + * @return MatrixX An Eigen3 matrix holding the concentrations and having + * the same dimensions as the grid. */ - const Eigen::MatrixXd &getConcentrations() { return this->concentrations; } + const Eigen::MatrixX &getConcentrations() { return this->concentrations; } /** * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one * dimensional. * - * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients. - * Matrix columns must have same size as length of grid. + * @param alpha An Eigen3 MatrixX with 1 row holding the alpha + * coefficients. Matrix columns must have same size as length of grid. */ - void setAlpha(const Eigen::MatrixXd &alpha); + void setAlpha(const Eigen::MatrixX &alpha) { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use 2D setter function!"); + } + if (alpha.rows() != 1 || alpha.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of alpha coefficients mismatch with Grid dimensions!"); + } + + this->alphaX = alpha; + } /** * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two * dimensional. * - * @param alphaX An Eigen3 MatrixXd holding the alpha coefficients in + * @param alphaX An Eigen3 MatrixX holding the alpha coefficients in * x-direction. Matrix must be of same size as the grid. - * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in + * @param alphaY An Eigen3 MatrixX holding the alpha coefficients in * y-direction. Matrix must be of same size as the grid. */ - void setAlpha(const Eigen::MatrixXd &alphaX, const Eigen::MatrixXd &alphaY); + void setAlpha(const Eigen::MatrixX &alphaX, + const Eigen::MatrixX &alphaY) { + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably " + "use 1D setter function!"); + } + if (alphaX.rows() != this->row || alphaX.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of alpha coefficients in x-direction " + "mismatch with GRid dimensions!"); + } + if (alphaY.rows() != this->row || alphaY.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of alpha coefficients in y-direction " + "mismatch with GRid dimensions!"); + } + + this->alphaX = alphaX; + this->alphaY = alphaY; + } /** * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one * dimensional. * - * @return MatrixXd A matrix with 1 row holding the alpha coefficients. + * @return MatrixX A matrix with 1 row holding the alpha coefficients. */ - const Eigen::MatrixXd &getAlpha(); + const Eigen::MatrixX &getAlpha() const { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use either getAlphaX() or getAlphaY()!"); + } + + return this->alphaX; + } /** * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. * Grid must be two dimensional. * - * @return MatrixXd A matrix holding the alpha coefficients in x-direction. + * @return MatrixX A matrix holding the alpha coefficients in x-direction. */ - const Eigen::MatrixXd &getAlphaX(); + const Eigen::MatrixX &getAlphaX() const { + + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably use getAlpha()!"); + } + + return this->alphaX; + } /** * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. * Grid must be two dimensional. * - * @return MatrixXd A matrix holding the alpha coefficients in y-direction. + * @return MatrixX A matrix holding the alpha coefficients in y-direction. */ - const Eigen::MatrixXd &getAlphaY(); + const Eigen::MatrixX &getAlphaY() const { + + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably use getAlpha()!"); + } + + return this->alphaY; + } /** * @brief Gets the dimensions of the grid. * * @return int Dimensions, either 1 or 2. */ - int getDim(); + int getDim() const { return this->dim; } /** * @brief Gets length of 1D grid. Must be one dimensional grid. * * @return int Length of 1D grid. */ - int getLength(); + int getLength() const { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use getRow() or getCol()!"); + } + + return col; + } /** * @brief Gets the number of rows of the grid. * * @return int Number of rows. */ - int getRow(); + int getRow() const { return this->row; } /** * @brief Gets the number of columns of the grid. * * @return int Number of columns. */ - int getCol(); + int getCol() const { return this->col; } /** * @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional. * * @param domainLength A double value of the domain length. Must be positive. */ - void setDomain(double domainLength); + void setDomain(double domainLength) { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probaly " + "use the 2D domain setter!"); + } + if (domainLength <= 0) { + throw std::invalid_argument("Given domain length is not positive!"); + } + + this->domainCol = domainLength; + this->deltaCol = double(this->domainCol) / double(this->col); + } /** * @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional. * - * @param domainRow A double value of the domain size in y-direction. Must be - * positive. - * @param domainCol A double value of the domain size in x-direction. Must be - * positive. + * @param domainRow A double value of the domain size in y-direction. Must + * be positive. + * @param domainCol A double value of the domain size in x-direction. Must + * be positive. */ - void setDomain(double domainRow, double domainCol); + void setDomain(double domainRow, double domainCol) { + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably " + "use the 1D domain setter!"); + } + if (domainRow <= 0 || domainCol <= 0) { + throw std::invalid_argument("Given domain size is not positive!"); + } + + this->domainRow = domainRow; + this->domainCol = domainCol; + this->deltaRow = double(this->domainRow) / double(this->row); + this->deltaCol = double(this->domainCol) / double(this->col); + } /** * @brief Gets the delta value for 1D-Grid. Grid must be one dimensional. * * @return double Delta value. */ - double getDelta(); + T getDelta() const { + + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use the 2D delta getters"); + } + + return this->deltaCol; + } /** * @brief Gets the delta value in x-direction. * * @return double Delta value in x-direction. */ - double getDeltaCol(); + T getDeltaCol() const { return this->deltaCol; } /** * @brief Gets the delta value in y-direction. Must be two dimensional grid. * * @return double Delta value in y-direction. */ - double getDeltaRow(); + T getDeltaRow() const { + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, meaning there is no " + "delta in y-direction!"); + } + + return this->deltaRow; + } private: - int col; // number of grid columns - int row{1}; // number of grid rows - int dim; // 1D or 2D - double domainCol; // number of domain columns - double domainRow{0}; // number of domain rows - double deltaCol; // delta in x-direction (between columns) - double deltaRow{0}; // delta in y-direction (between rows) - Eigen::MatrixXd concentrations; // Matrix holding grid concentrations - Eigen::MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction - Eigen::MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction + int col; // number of grid columns + int row{1}; // number of grid rows + int dim; // 1D or 2D + T domainCol; // number of domain columns + T domainRow{0}; // number of domain rows + T deltaCol; // delta in x-direction (between columns) + T deltaRow{0}; // delta in y-direction (between rows) + Eigen::MatrixX concentrations; // Matrix holding grid concentrations + Eigen::MatrixX alphaX; // Matrix holding alpha coefficients in x-direction + Eigen::MatrixX alphaY; // Matrix holding alpha coefficients in y-direction + + static constexpr double MAT_INIT_VAL = 0; }; +using Grid64 = Grid; +using Grid32 = Grid; + #endif // GRID_H_ diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 0c09a7e..7ed50c4 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -11,6 +11,11 @@ #include "Boundary.hpp" #include "Grid.hpp" +#include +#include +#include +#include +#include #include #include @@ -77,7 +82,8 @@ enum TIME_MEASURE { * time step, number of simulations, etc. * */ -class Simulation { + +template class Simulation { public: /** * @brief Set up a simulation environment. The timestep and number of @@ -91,7 +97,8 @@ public: * @param bc Valid boundary condition object * @param approach Approach to solving the problem. Either FTCS or BTCS. */ - Simulation(Grid &grid, Boundary &bc, APPROACH approach); + Simulation(Grid &_grid, Boundary &_bc, APPROACH _approach) + : grid(_grid), bc(_bc), approach(_approach){}; /** * @brief Set the option to output the results to a CSV file. Off by default. @@ -107,7 +114,13 @@ public: * - CSV_OUTPUT_XTREME: produce csv output with all * concentration matrices and simulation environment */ - void setOutputCSV(CSV_OUTPUT csv_output); + void setOutputCSV(CSV_OUTPUT csv_output) { + if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { + throw std::invalid_argument("Invalid CSV output option given!"); + } + + this->csv_output = csv_output; + } /** * @brief Set the options for outputting information to the console. Off by @@ -122,7 +135,14 @@ public: * - CONSOLE_OUTPUT_VERBOSE: print all concentration * matrices to console */ - void setOutputConsole(CONSOLE_OUTPUT console_output); + void setOutputConsole(CONSOLE_OUTPUT console_output) { + if (console_output < CONSOLE_OUTPUT_OFF && + console_output > CONSOLE_OUTPUT_VERBOSE) { + throw std::invalid_argument("Invalid console output option given!"); + } + + this->console_output = console_output; + } /** * @brief Set the Time Measure option. Off by default. @@ -133,7 +153,13 @@ public: * - TIME_MEASURE_ON: Time of simulation run is printed to * console */ - void setTimeMeasure(TIME_MEASURE time_measure); + void setTimeMeasure(TIME_MEASURE time_measure) { + if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) { + throw std::invalid_argument("Invalid time measure option given!"); + } + + this->time_measure = time_measure; + } /** * @brief Setting the time step for each iteration step. Time step must be @@ -141,14 +167,14 @@ public: * * @param timestep Valid timestep greater than zero. */ - void setTimestep(double timestep); + void setTimestep(T timestep); /** * @brief Currently set time step is returned. * * @return double timestep */ - double getTimestep(); + T getTimestep() const { return this->timestep; } /** * @brief Set the desired iterations to be calculated. A value greater @@ -156,7 +182,13 @@ public: * * @param iterations Number of iterations to be simulated. */ - void setIterations(int iterations); + void setIterations(int iterations) { + if (iterations <= 0) { + throw std::invalid_argument( + "Number of iterations must be greater than zero."); + } + this->iterations = iterations; + } /** * @brief Set the desired linear equation solver to be used for BTCS approach. @@ -165,7 +197,17 @@ public: * @param solver Solver to be used. Default is Thomas Algorithm as it is more * efficient for tridiagonal Matrices. */ - void setSolver(SOLVER solver); + void setSolver(SOLVER solver) { + if (this->approach == FTCS_APPROACH) { + std::cerr + << "Warning: Solver was set, but FTCS approach initialized. Setting " + "the solver " + "is thus without effect." + << std::endl; + } + + this->solver = solver; + } /** * @brief Set the number of desired openMP Threads. @@ -175,20 +217,32 @@ public: * maximum number of processors is set as the default case during Simulation * construction. */ - void setNumberThreads(int num_threads); + void setNumberThreads(int num_threads) { + if (numThreads > 0 && numThreads <= omp_get_num_procs()) { + this->numThreads = numThreads; + } else { + int maxThreadNumber = omp_get_num_procs(); + throw std::invalid_argument( + "Number of threads exceeds the number of processor cores (" + + std::to_string(maxThreadNumber) + ") or is less than 1."); + } + } /** * @brief Return the currently set iterations to be calculated. * * @return int Number of iterations. */ - int getIterations(); + int getIterations() const { return this->iterations; } /** * @brief Outputs the current concentrations of the grid on the console. * */ - void printConcentrationsConsole(); + inline void printConcentrationsConsole() const { + std::cout << grid.getConcentrations() << std::endl; + std::cout << std::endl; + } /** * @brief Creates a CSV file with a name containing the current simulation @@ -199,7 +253,52 @@ public: * * @return string Filename with configured simulation parameters. */ - std::string createCSVfile(); + std::string createCSVfile() const { + std::ofstream file; + int appendIdent = 0; + std::string appendIdentString; + + // string approachString = (approach == 0) ? "FTCS" : "BTCS"; + const std::string &approachString = this->approach_names[approach]; + std::string row = std::to_string(grid.getRow()); + std::string col = std::to_string(grid.getCol()); + std::string numIterations = std::to_string(iterations); + + std::string filename = + approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; + + while (std::filesystem::exists(filename)) { + appendIdent += 1; + appendIdentString = std::to_string(appendIdent); + filename = approachString + "_" + row + "_" + col + "_" + numIterations + + "-" + appendIdentString + ".csv"; + } + + file.open(filename); + if (!file) { + exit(1); + } + + // adds lines at the beginning of verbose output csv that represent the + // boundary conditions and their values -1 in case of closed boundary + if (csv_output == CSV_OUTPUT_XTREME) { + Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "", + " "); + file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) + << std::endl; // boundary left + file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) + << std::endl; // boundary right + file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) + << std::endl; // boundary top + file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) + << std::endl; // boundary bottom + file << std::endl << std::endl; + } + + file.close(); + + return filename; + } /** * @brief Writes the currently calculated concentration values of the grid @@ -208,7 +307,19 @@ public: * @param filename Name of the file to which the concentration values are * to be written. */ - void printConcentrationsCSV(const std::string &filename); + void printConcentrationsCSV(const std::string &filename) const { + std::ofstream file; + + file.open(filename, std::ios_base::app); + if (!file) { + exit(1); + } + + Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); + file << grid.getConcentrations().format(do_not_align) << std::endl; + file << std::endl << std::endl; + file.close(); + } /** * @brief Method starts the simulation process with the previously set @@ -217,7 +328,7 @@ public: void run(); private: - double timestep{-1}; + T timestep{-1}; int iterations{-1}; int innerIterations{1}; int numThreads{omp_get_num_procs()}; @@ -225,8 +336,8 @@ private: CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT_OFF}; TIME_MEASURE time_measure{TIME_MEASURE_OFF}; - Grid &grid; - Boundary &bc; + Grid &grid; + Boundary &bc; APPROACH approach; SOLVER solver{THOMAS_ALGORITHM_SOLVER}; diff --git a/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp index fe15455..bc80a8b 100644 --- a/naaice/BTCS_2D_NAAICE.cpp +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -103,23 +103,23 @@ int main(int argc, char *argv[]) { // create a grid with a 5 x 10 field constexpr int row = 5; constexpr int col = 10; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: grid.setDomain(0.005, 0.01); const auto init_values_vec = CSVToVector(INPUT_CONC_FILE); - MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col); + Eigen::MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col); grid.setConcentrations(concentrations); const double sum_init = concentrations.sum(); // // (optional) set alphas of the grid, e.g.: const auto alphax_vec = CSVToVector(INPUT_ALPHAX_FILE); - MatrixXd alphax = rmVecTocmMatrix(alphax_vec, row, col); + Eigen::MatrixXd alphax = rmVecTocmMatrix(alphax_vec, row, col); constexpr double alphay_val = 5e-10; - MatrixXd alphay = MatrixXd::Constant(row, col, alphay_val); // row,col,value + Eigen::MatrixXd alphay = Eigen::MatrixXd::Constant(row, col, alphay_val); // row,col,value grid.setAlpha(alphax, alphay); // // ****************** diff --git a/src/BTCS.cpp b/src/BTCS.cpp index 3b314e7..7a9dc7e 100644 --- a/src/BTCS.cpp +++ b/src/BTCS.cpp @@ -45,13 +45,15 @@ constexpr std::pair calcBoundaryCoeffClosed(T alpha_center, T alpha_side, } // creates coefficient matrix for next time step from alphas in x-direction -static Eigen::SparseMatrix -createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector &bcLeft, - std::vector &bcRight, int numCols, - int rowIndex, double sx) { +template +static Eigen::SparseMatrix +createCoeffMatrix(const Eigen::MatrixX &alpha, + const std::vector> &bcLeft, + const std::vector> &bcRight, int numCols, + int rowIndex, T sx) { // square matrix of column^2 dimension for the coefficients - Eigen::SparseMatrix cm(numCols, numCols); + Eigen::SparseMatrix cm(numCols, numCols); cm.reserve(Eigen::VectorXi::Constant(numCols, 3)); // left column @@ -142,14 +144,18 @@ constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc, // creates a solution vector for next time step from the current state of // concentrations -static Eigen::VectorXd createSolutionVector( - Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alphaX, - Eigen::MatrixXd &alphaY, std::vector &bcLeft, - std::vector &bcRight, std::vector &bcTop, - std::vector &bcBottom, int length, int rowIndex, double sx, - double sy) { +template +static Eigen::VectorX +createSolutionVector(const Eigen::MatrixX &concentrations, + const Eigen::MatrixX &alphaX, + const Eigen::MatrixX &alphaY, + const std::vector> &bcLeft, + const std::vector> &bcRight, + const std::vector> &bcTop, + const std::vector> &bcBottom, + int length, int rowIndex, T sx, T sy) { - Eigen::VectorXd sv(length); + Eigen::VectorX sv(length); const std::size_t numRows = concentrations.rows(); // inner rows @@ -235,10 +241,11 @@ static Eigen::VectorXd createSolutionVector( // solver for linear equation system; A corresponds to coefficient matrix, // b to the solution vector // use of EigenLU solver -static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix &A, - Eigen::VectorXd &b) { +template +static Eigen::VectorX EigenLUAlgorithm(Eigen::SparseMatrix &A, + Eigen::VectorX &b) { - Eigen::SparseLU> solver; + Eigen::SparseLU> solver; solver.analyzePattern(A); solver.factorize(A); @@ -248,14 +255,15 @@ static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix &A, // solver for linear equation system; A corresponds to coefficient matrix, // b to the solution vector // implementation of Thomas Algorithm -static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix &A, - Eigen::VectorXd &b) { +template +static Eigen::VectorX ThomasAlgorithm(Eigen::SparseMatrix &A, + Eigen::VectorX &b) { Eigen::Index n = b.size(); - Eigen::VectorXd a_diag(n); - Eigen::VectorXd b_diag(n); - Eigen::VectorXd c_diag(n); - Eigen::VectorXd x_vec = b; + Eigen::VectorX a_diag(n); + Eigen::VectorX b_diag(n); + Eigen::VectorX c_diag(n); + Eigen::VectorX x_vec = b; // Fill diagonals vectors b_diag[0] = A.coeff(0, 0); @@ -291,23 +299,24 @@ static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix &A, } // BTCS solution for 1D grid -static void -BTCS_1D(Grid &grid, Boundary &bc, double timestep, - Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix &A, - Eigen::VectorXd &b)) { +template +static void BTCS_1D(Grid &grid, Boundary &bc, T timestep, + Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, + Eigen::VectorX &b)) { int length = grid.getLength(); - double sx = timestep / (grid.getDelta() * grid.getDelta()); + T sx = timestep / (grid.getDelta() * grid.getDelta()); - Eigen::VectorXd concentrations_t1(length); + Eigen::VectorX concentrations_t1(length); - Eigen::SparseMatrix A; - Eigen::VectorXd b(length); + Eigen::SparseMatrix A; + Eigen::VectorX b(length); - Eigen::MatrixXd alpha = grid.getAlpha(); - std::vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - std::vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + const auto &alpha = grid.getAlpha(); - Eigen::MatrixXd concentrations = grid.getConcentrations(); + const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + + Eigen::MatrixX concentrations = grid.getConcentrations(); int rowIndex = 0; A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx); // this is exactly same as in 2D @@ -331,31 +340,32 @@ BTCS_1D(Grid &grid, Boundary &bc, double timestep, } // BTCS solution for 2D grid -static void -BTCS_2D(Grid &grid, Boundary &bc, double timestep, - Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix &A, - Eigen::VectorXd &b), - int numThreads) { +template +static void BTCS_2D(Grid &grid, Boundary &bc, T timestep, + Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, + Eigen::VectorX &b), + int numThreads) { int rowMax = grid.getRow(); int colMax = grid.getCol(); - double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); - double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); + T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); + T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); - Eigen::MatrixXd concentrations_t1 = - Eigen::MatrixXd::Constant(rowMax, colMax, 0); - Eigen::VectorXd row_t1(colMax); + Eigen::MatrixX concentrations_t1 = + Eigen::MatrixX::Constant(rowMax, colMax, 0); + Eigen::VectorX row_t1(colMax); - Eigen::SparseMatrix A; - Eigen::VectorXd b; + Eigen::SparseMatrix A; + Eigen::VectorX b; - Eigen::MatrixXd alphaX = grid.getAlphaX(); - Eigen::MatrixXd alphaY = grid.getAlphaY(); - std::vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - std::vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); - std::vector bcTop = bc.getBoundarySide(BC_SIDE_TOP); - std::vector bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); + auto alphaX = grid.getAlphaX(); + auto alphaY = grid.getAlphaY(); - Eigen::MatrixXd concentrations = grid.getConcentrations(); + const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + const auto &bcTop = bc.getBoundarySide(BC_SIDE_TOP); + const auto &bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); + + Eigen::MatrixX concentrations = grid.getConcentrations(); #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) for (int i = 0; i < rowMax; i++) { @@ -364,8 +374,6 @@ BTCS_2D(Grid &grid, Boundary &bc, double timestep, b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, colMax, i, sx, sy); - Eigen::SparseLU> solver; - row_t1 = solverFunc(A, b); concentrations_t1.row(i) = row_t1; @@ -394,7 +402,8 @@ BTCS_2D(Grid &grid, Boundary &bc, double timestep, } // entry point for EigenLU solver; differentiate between 1D and 2D grid -void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { +template +void BTCS_LU(Grid &grid, Boundary &bc, T timestep, int numThreads) { if (grid.getDim() == 1) { BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); } else if (grid.getDim() == 2) { @@ -406,7 +415,8 @@ void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { } // entry point for Thomas algorithm solver; differentiate 1D and 2D grid -void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) { +template +void BTCS_Thomas(Grid &grid, Boundary &bc, T timestep, int numThreads) { if (grid.getDim() == 1) { BTCS_1D(grid, bc, timestep, ThomasAlgorithm); } else if (grid.getDim() == 2) { @@ -416,3 +426,13 @@ void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) { "Error: Only 1- and 2-dimensional grids are defined!"); } } + +template void BTCS_Thomas(Grid &grid, Boundary &bc, + double timestep, int numThreads); +template void BTCS_Thomas(Grid &grid, Boundary &bc, + float timestep, int numThreads); + +template void BTCS_LU(Grid &grid, Boundary &bc, double timestep, + int numThreads); +template void BTCS_LU(Grid &grid, Boundary &bc, float timestep, + int numThreads); diff --git a/src/Boundary.cpp b/src/Boundary.cpp deleted file mode 100644 index def9bbd..0000000 --- a/src/Boundary.cpp +++ /dev/null @@ -1,149 +0,0 @@ -#include "TugUtils.hpp" - -#include -#include -#include - -BoundaryElement::BoundaryElement() {} - -BoundaryElement::BoundaryElement(double _value) - : value(_value), type(BC_TYPE_CONSTANT) {} - -void BoundaryElement::setType(BC_TYPE type) { this->type = type; } - -void BoundaryElement::setValue(double value) { - if (value < 0) { - throw_invalid_argument("No negative concentration allowed."); - } - if (type == BC_TYPE_CLOSED) { - throw_invalid_argument( - "No constant boundary concentrations can be set for closed " - "boundaries. Please change type first."); - } - this->value = value; -} - -BC_TYPE BoundaryElement::getType() { return this->type; } - -double BoundaryElement::getValue() { return this->value; } - -Boundary::Boundary(Grid grid) : grid(grid) { - if (grid.getDim() == 1) { - this->boundaries = std::vector>( - 2); // in 1D only left and right boundary - - this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); - this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); - } else if (grid.getDim() == 2) { - this->boundaries = std::vector>(4); - - this->boundaries[BC_SIDE_LEFT] = - std::vector(grid.getRow(), BoundaryElement()); - this->boundaries[BC_SIDE_RIGHT] = - std::vector(grid.getRow(), BoundaryElement()); - this->boundaries[BC_SIDE_TOP] = - std::vector(grid.getCol(), BoundaryElement()); - this->boundaries[BC_SIDE_BOTTOM] = - std::vector(grid.getCol(), BoundaryElement()); - } -} - -void Boundary::setBoundarySideClosed(BC_SIDE side) { - if (grid.getDim() == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw_invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } - - const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; - const int n = is_vertical ? grid.getRow() : grid.getCol(); - - this->boundaries[side] = std::vector(n, BoundaryElement()); -} - -void Boundary::setBoundarySideConstant(BC_SIDE side, double value) { - if (grid.getDim() == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw_invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } - - const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; - const int n = is_vertical ? grid.getRow() : grid.getCol(); - - this->boundaries[side] = - std::vector(n, BoundaryElement(value)); -} - -void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) { - // tests whether the index really points to an element of the boundary side. - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - this->boundaries[side][index].setType(BC_TYPE_CLOSED); -} - -void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, - double value) { - // tests whether the index really points to an element of the boundary side. - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - this->boundaries[side][index].setType(BC_TYPE_CONSTANT); - this->boundaries[side][index].setValue(value); -} - -const std::vector Boundary::getBoundarySide(BC_SIDE side) { - if (grid.getDim() == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw_invalid_argument( - "For the one-dimensional trap, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } - return this->boundaries[side]; -} - -Eigen::VectorXd Boundary::getBoundarySideValues(BC_SIDE side) { - const std::size_t length = boundaries[side].size(); - Eigen::VectorXd values(length); - - for (int i = 0; i < length; i++) { - if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { - values(i) = -1; - continue; - } - values(i) = getBoundaryElementValue(side, i); - } - - return values; -} - -BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) { - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - return this->boundaries[side][index]; -} - -BC_TYPE Boundary::getBoundaryElementType(BC_SIDE side, int index) { - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - return this->boundaries[side][index].getType(); -} - -double Boundary::getBoundaryElementValue(BC_SIDE side, int index) { - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) { - throw_invalid_argument( - "A value can only be output if it is a constant boundary condition."); - } - return this->boundaries[side][index].getValue(); -} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5527d09..91d141d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,4 @@ -add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCS.cpp) +add_library(tug Simulation.cpp FTCS.cpp BTCS.cpp) IF(TUG_NAAICE_EXAMPLE) target_compile_definitions(tug PRIVATE WRITE_THOMAS_CSV) diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 375fe50..5ee4657 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -19,7 +19,8 @@ #endif // calculates horizontal change on one cell independent of boundary type -static inline double calcHorizontalChange(Grid &grid, int &row, int &col) { +template +static inline T calcHorizontalChange(Grid &grid, int &row, int &col) { return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), grid.getAlphaX()(row, col)) * @@ -35,7 +36,8 @@ static inline double calcHorizontalChange(Grid &grid, int &row, int &col) { } // calculates vertical change on one cell independent of boundary type -static inline double calcVerticalChange(Grid &grid, int &row, int &col) { +template +static inline T calcVerticalChange(Grid &grid, int &row, int &col) { return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), grid.getAlphaY()(row, col)) * @@ -51,10 +53,10 @@ static inline double calcVerticalChange(Grid &grid, int &row, int &col) { } // calculates horizontal change on one cell with a constant left boundary -static inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, - int &col) { +template +static inline T calcHorizontalChangeLeftBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, int &col) { return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), grid.getAlphaX()(row, col)) * @@ -68,8 +70,9 @@ static inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, } // calculates horizontal change on one cell with a closed left boundary -static inline double -calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) { +template +static inline T calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, + int &col) { return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), grid.getAlphaX()(row, col)) * @@ -78,8 +81,9 @@ calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) { } // checks boundary condition type for a cell on the left edge of grid -static inline double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T 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) { @@ -90,10 +94,10 @@ static inline double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, } // calculates horizontal change on one cell with a constant right boundary -static inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, - int &col) { +template +static inline T calcHorizontalChangeRightBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, int &col) { return 2 * grid.getAlphaX()(row, col) * bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) - @@ -107,8 +111,9 @@ static inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid, } // calculates horizontal change on one cell with a closed right boundary -static inline double -calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) { +template +static inline T calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, + int &col) { return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1), grid.getAlphaX()(row, col)) * @@ -117,8 +122,10 @@ calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) { } // checks boundary condition type for a cell on the right edge of grid -static inline double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T 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) { @@ -129,9 +136,10 @@ static inline double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, } // calculates vertical change on one cell with a constant top boundary -static inline double calcVerticalChangeTopBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, int &col) { +template +static inline T calcVerticalChangeTopBoundaryConstant(Grid &grid, + Boundary &bc, int &row, + int &col) { return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), grid.getAlphaY()(row, col)) * @@ -145,8 +153,9 @@ static inline double calcVerticalChangeTopBoundaryConstant(Grid &grid, } // calculates vertical change on one cell with a closed top boundary -static inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, - int &col) { +template +static inline T calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, + int &col) { return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), grid.getConcentrations()(row, col)) * @@ -155,8 +164,9 @@ static inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, } // checks boundary condition type for a cell on the top edge of grid -static inline double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T 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) { @@ -167,10 +177,10 @@ static inline double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, } // calculates vertical change on one cell with a constant bottom boundary -static inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, - int &col) { +template +static inline T calcVerticalChangeBottomBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, int &col) { return 2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) - @@ -184,8 +194,9 @@ static inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid, } // calculates vertical change on one cell with a closed bottom boundary -static inline double -calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) { +template +static inline T calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, + int &col) { return -(calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row - 1, col)) * @@ -194,8 +205,9 @@ calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) { } // checks boundary condition type for a cell on the bottom edge of grid -static inline double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T 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) { @@ -206,12 +218,14 @@ static inline double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, } // FTCS solution for 1D grid -static void FTCS_1D(Grid &grid, Boundary &bc, double ×tep) { +template +static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { int colMax = grid.getCol(); - double deltaCol = grid.getDeltaCol(); + T deltaCol = grid.getDeltaCol(); // matrix for concentrations at time t+1 - Eigen::MatrixXd concentrations_t1 = Eigen::MatrixXd::Constant(1, colMax, 0); + Eigen::MatrixX concentrations_t1 = + Eigen::MatrixX::Constant(1, colMax, 0); // only one row in 1D case -> row constant at index 0 int row = 0; @@ -243,16 +257,17 @@ static void FTCS_1D(Grid &grid, Boundary &bc, double ×tep) { } // FTCS solution for 2D grid -static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep, +template +static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, int numThreads) { int rowMax = grid.getRow(); int colMax = grid.getCol(); - double deltaRow = grid.getDeltaRow(); - double deltaCol = grid.getDeltaCol(); + T deltaRow = grid.getDeltaRow(); + T deltaCol = grid.getDeltaCol(); // matrix for concentrations at time t+1 - Eigen::MatrixXd concentrations_t1 = - Eigen::MatrixXd::Constant(rowMax, colMax, 0); + Eigen::MatrixX concentrations_t1 = + Eigen::MatrixX::Constant(rowMax, colMax, 0); // inner cells // these are independent of the boundary condition type @@ -369,7 +384,8 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep, } // entry point; differentiate between 1D and 2D grid -void FTCS(Grid &grid, Boundary &bc, double ×tep, int &numThreads) { +template +void FTCS(Grid &grid, Boundary &bc, T timestep, int &numThreads) { if (grid.getDim() == 1) { FTCS_1D(grid, bc, timestep); } else if (grid.getDim() == 2) { @@ -379,3 +395,8 @@ void FTCS(Grid &grid, Boundary &bc, double ×tep, int &numThreads) { "Error: Only 1- and 2-dimensional grids are defined!"); } } + +template void FTCS(Grid &grid, Boundary &bc, double timestep, + int &numThreads); +template void FTCS(Grid &grid, Boundary &bc, float timestep, + int &numThreads); diff --git a/src/Grid.cpp b/src/Grid.cpp deleted file mode 100644 index c6573e4..0000000 --- a/src/Grid.cpp +++ /dev/null @@ -1,167 +0,0 @@ -#include "TugUtils.hpp" - -#include -#include - -constexpr double MAT_INIT_VAL = 0; - -Grid::Grid(int length) : col(length), domainCol(length) { - if (length <= 3) { - throw_invalid_argument( - "Given grid length too small. Must be greater than 3."); - } - - this->dim = 1; - this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 - - this->concentrations = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); - this->alphaX = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); -} - -Grid::Grid(int _row, int _col) - : row(_row), col(_col), domainRow(_row), domainCol(_col) { - if (row <= 3 || col <= 3) { - throw_invalid_argument( - "Given grid dimensions too small. Must each be greater than 3."); - } - - this->dim = 2; - this->deltaRow = double(this->domainRow) / double(this->row); // -> 1 - this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 - - this->concentrations = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); - this->alphaX = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); - this->alphaY = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); -} - -void Grid::setConcentrations(const Eigen::MatrixXd &concentrations) { - if (concentrations.rows() != this->row || - concentrations.cols() != this->col) { - throw_invalid_argument( - "Given matrix of concentrations mismatch with Grid dimensions!"); - } - - this->concentrations = concentrations; -} - -void Grid::setAlpha(const Eigen::MatrixXd &alpha) { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use 2D setter function!"); - } - if (alpha.rows() != 1 || alpha.cols() != this->col) { - throw_invalid_argument( - "Given matrix of alpha coefficients mismatch with Grid dimensions!"); - } - - this->alphaX = alpha; -} - -void Grid::setAlpha(const Eigen::MatrixXd &alphaX, - const Eigen::MatrixXd &alphaY) { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably " - "use 1D setter function!"); - } - if (alphaX.rows() != this->row || alphaX.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients in x-direction " - "mismatch with GRid dimensions!"); - } - if (alphaY.rows() != this->row || alphaY.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients in y-direction " - "mismatch with GRid dimensions!"); - } - - this->alphaX = alphaX; - this->alphaY = alphaY; -} - -const Eigen::MatrixXd &Grid::getAlpha() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use either getAlphaX() or getAlphaY()!"); - } - - return this->alphaX; -} - -const Eigen::MatrixXd &Grid::getAlphaX() { - if (dim != 2) { - throw_invalid_argument( - "Grid is not two dimensional, you should probably use getAlpha()!"); - } - - return this->alphaX; -} - -const Eigen::MatrixXd &Grid::getAlphaY() { - if (dim != 2) { - throw_invalid_argument( - "Grid is not two dimensional, you should probably use getAlpha()!"); - } - - return this->alphaY; -} - -int Grid::getDim() { return dim; } - -int Grid::getLength() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use getRow() or getCol()!"); - } - - return col; -} - -int Grid::getRow() { return row; } - -int Grid::getCol() { return col; } - -void Grid::setDomain(double domainLength) { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probaly " - "use the 2D domain setter!"); - } - if (domainLength <= 0) { - throw_invalid_argument("Given domain length is not positive!"); - } - - this->domainCol = domainLength; - this->deltaCol = double(this->domainCol) / double(this->col); -} - -void Grid::setDomain(double domainRow, double domainCol) { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably " - "use the 1D domain setter!"); - } - if (domainRow <= 0 || domainCol <= 0) { - throw_invalid_argument("Given domain size is not positive!"); - } - - this->domainRow = domainRow; - this->domainCol = domainCol; - this->deltaRow = double(this->domainRow) / double(this->row); - this->deltaCol = double(this->domainCol) / double(this->col); -} - -double Grid::getDelta() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use the 2D delta getters"); - } - - return this->deltaCol; -} - -double Grid::getDeltaCol() { return this->deltaCol; } - -double Grid::getDeltaRow() { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, meaning there is no " - "delta in y-direction!"); - } - - return this->deltaRow; -} diff --git a/src/Schemes.hpp b/src/Schemes.hpp index d10e0c5..771a2a4 100644 --- a/src/Schemes.hpp +++ b/src/Schemes.hpp @@ -16,13 +16,16 @@ #include // entry point; differentiate between 1D and 2D grid -extern void FTCS(Grid &grid, Boundary &bc, double ×tep, int &numThreads); +template +extern void FTCS(Grid &grid, Boundary &bc, T timestep, int &numThreads); // entry point for EigenLU solver; differentiate between 1D and 2D grid -extern void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads); +template +extern void BTCS_LU(Grid &grid, Boundary &bc, T timestep, int numThreads); // entry point for Thomas algorithm solver; differentiate 1D and 2D grid -extern void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, +template +extern void BTCS_Thomas(Grid &grid, Boundary &bc, T timestep, int numThreads); #endif // SCHEMES_H_ diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 0a5a961..90a1459 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -12,56 +12,28 @@ #include "Schemes.hpp" #include "TugUtils.hpp" -Simulation::Simulation(Grid &_grid, Boundary &_bc, APPROACH _approach) - : grid(_grid), bc(_bc), approach(_approach) {} - -void Simulation::setOutputCSV(CSV_OUTPUT csv_output) { - if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { - throw_invalid_argument("Invalid CSV output option given!"); - } - - this->csv_output = csv_output; -} - -void Simulation::setOutputConsole(CONSOLE_OUTPUT console_output) { - if (console_output < CONSOLE_OUTPUT_OFF && - console_output > CONSOLE_OUTPUT_VERBOSE) { - throw_invalid_argument("Invalid console output option given!"); - } - - this->console_output = console_output; -} - -void Simulation::setTimeMeasure(TIME_MEASURE time_measure) { - if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) { - throw_invalid_argument("Invalid time measure option given!"); - } - - this->time_measure = time_measure; -} - -void Simulation::setTimestep(double timestep) { +template void Simulation::setTimestep(T timestep) { if (timestep <= 0) { throw_invalid_argument("Timestep has to be greater than zero."); } if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { - const double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); // will be 0 if 1D, else ... - const double deltaRowSquare = grid.getDim() != 1 - ? grid.getDeltaRow() * grid.getDeltaRow() - : deltaColSquare; - const double minDeltaSquare = + const T deltaRowSquare = grid.getDim() != 1 + ? grid.getDeltaRow() * grid.getDeltaRow() + : deltaColSquare; + const T minDeltaSquare = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - double maxAlpha = std::numeric_limits::quiet_NaN(); + T maxAlpha = std::numeric_limits::quiet_NaN(); // determine maximum alpha if (grid.getDim() == 2) { - const double maxAlphaX = grid.getAlphaX().maxCoeff(); - const double maxAlphaY = grid.getAlphaY().maxCoeff(); + const T maxAlphaX = grid.getAlphaX().maxCoeff(); + const T maxAlphaY = grid.getAlphaY().maxCoeff(); maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; } else if (grid.getDim() == 1) { @@ -74,7 +46,7 @@ void Simulation::setTimestep(double timestep) { const std::string dim = std::to_string(grid.getDim()) + "D"; // Courant-Friedrichs-Lewy condition - double cfl = minDeltaSquare / (4 * maxAlpha); + T cfl = minDeltaSquare / (4 * maxAlpha); // stability equation from Wikipedia; might be useful if applied cfl does // not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha * @@ -112,109 +84,7 @@ void Simulation::setTimestep(double timestep) { } } -double Simulation::getTimestep() { return this->timestep; } - -void Simulation::setIterations(int iterations) { - if (iterations <= 0) { - throw_invalid_argument("Number of iterations must be greater than zero."); - } - this->iterations = iterations; -} - -void Simulation::setSolver(SOLVER solver) { - if (this->approach == FTCS_APPROACH) { - std::cerr - << "Warning: Solver was set, but FTCS approach initialized. Setting " - "the solver " - "is thus without effect." - << std::endl; - } - - this->solver = solver; -} - -void Simulation::setNumberThreads(int numThreads) { - if (numThreads > 0 && numThreads <= omp_get_num_procs()) { - this->numThreads = numThreads; - } else { - int maxThreadNumber = omp_get_num_procs(); - std::string outputMessage = - "Number of threads exceeds the number of processor cores (" + - std::to_string(maxThreadNumber) + ") or is less than 1."; - - throw_invalid_argument(outputMessage); - } -} - -int Simulation::getIterations() { return this->iterations; } - -void Simulation::printConcentrationsConsole() { - std::cout << grid.getConcentrations() << std::endl; - std::cout << std::endl; -} - -std::string Simulation::createCSVfile() { - std::ofstream file; - int appendIdent = 0; - std::string appendIdentString; - - // string approachString = (approach == 0) ? "FTCS" : "BTCS"; - const std::string &approachString = this->approach_names[approach]; - std::string row = std::to_string(grid.getRow()); - std::string col = std::to_string(grid.getCol()); - std::string numIterations = std::to_string(iterations); - - std::string filename = - approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; - - while (std::filesystem::exists(filename)) { - appendIdent += 1; - appendIdentString = std::to_string(appendIdent); - filename = approachString + "_" + row + "_" + col + "_" + numIterations + - "-" + appendIdentString + ".csv"; - } - - file.open(filename); - if (!file) { - exit(1); - } - - // adds lines at the beginning of verbose output csv that represent the - // boundary conditions and their values -1 in case of closed boundary - if (csv_output == CSV_OUTPUT_XTREME) { - Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "", - " "); - file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) - << std::endl; // boundary left - file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) - << std::endl; // boundary right - file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) - << std::endl; // boundary top - file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) - << std::endl; // boundary bottom - file << std::endl << std::endl; - } - - file.close(); - - return filename; -} - -void Simulation::printConcentrationsCSV(const std::string &filename) { - std::ofstream file; - - file.open(filename, std::ios_base::app); - if (!file) { - exit(1); - } - - Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); - file << grid.getConcentrations().format(do_not_align) << std::endl; - file << std::endl << std::endl; - file.close(); -} - -void Simulation::run() { +template void Simulation::run() { if (this->timestep == -1) { throw_invalid_argument("Timestep is not set!"); } @@ -280,14 +150,14 @@ void Simulation::run() { } else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case - double beta = 0.5; + constexpr T beta = 0.5; // TODO this implementation is very inefficient! // a separate implementation that sets up a specific tridiagonal matrix for // Crank-Nicolson would be better - Eigen::MatrixXd concentrations; - Eigen::MatrixXd concentrationsFTCS; - Eigen::MatrixXd concentrationsResult; + Eigen::MatrixX concentrations; + Eigen::MatrixX concentrationsFTCS; + Eigen::MatrixX concentrationsResult; for (int i = 0; i < iterations * innerIterations; i++) { if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { printConcentrationsConsole(); @@ -324,3 +194,9 @@ void Simulation::run() { << milliseconds.count() << "ms" << std::endl; } } + +template void Simulation::setTimestep(double timestep); +template void Simulation::setTimestep(float timestep); + +template void Simulation::run(); +template void Simulation::run(); diff --git a/test/testBoundary.cpp b/test/testBoundary.cpp index 1d2bff2..be91b1e 100644 --- a/test/testBoundary.cpp +++ b/test/testBoundary.cpp @@ -6,12 +6,11 @@ #include using namespace std; - TEST_CASE("BoundaryElement") { SUBCASE("Closed case") { - BoundaryElement boundaryElementClosed = BoundaryElement(); - CHECK_NOTHROW(BoundaryElement()); + BoundaryElement boundaryElementClosed = BoundaryElement(); + CHECK_NOTHROW(BoundaryElement()); CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED); CHECK_EQ(boundaryElementClosed.getValue(), -1); CHECK_THROWS(boundaryElementClosed.setValue(0.2)); @@ -28,11 +27,11 @@ TEST_CASE("BoundaryElement") { } TEST_CASE("Boundary Class") { - Grid grid1D = Grid(10); - Grid grid2D = Grid(10, 12); + Grid grid1D = Grid64(10); + Grid grid2D = Grid64(10, 12); Boundary boundary1D = Boundary(grid1D); Boundary boundary2D = Boundary(grid2D); - vector boundary1DVector(1, BoundaryElement(1.0)); + vector> boundary1DVector(1, BoundaryElement(1.0)); SUBCASE("Boundaries 1D case") { CHECK_NOTHROW(Boundary boundary(grid1D)); diff --git a/test/testGrid.cpp b/test/testGrid.cpp index a1f2886..ba831d3 100644 --- a/test/testGrid.cpp +++ b/test/testGrid.cpp @@ -8,22 +8,22 @@ using namespace std; TEST_CASE("1D Grid, too small length") { int l = 2; - CHECK_THROWS(Grid(l)); + CHECK_THROWS(Grid64(l)); } -TEST_CASE("2D Grid, too small side") { +TEST_CASE("2D Grid64, too small side") { int r = 2; int c = 4; - CHECK_THROWS(Grid(r, c)); + CHECK_THROWS(Grid64(r, c)); r = 4; c = 2; - CHECK_THROWS(Grid(r, c)); + CHECK_THROWS(Grid64(r, c)); } -TEST_CASE("1D Grid") { +TEST_CASE("1D Grid64") { int l = 12; - Grid grid(l); + Grid64 grid(l); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 1); @@ -89,9 +89,9 @@ TEST_CASE("1D Grid") { } } -TEST_CASE("2D Grid quadratic") { +TEST_CASE("2D Grid64 quadratic") { int rc = 12; - Grid grid(rc, rc); + Grid64 grid(rc, rc); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 2); @@ -170,10 +170,10 @@ TEST_CASE("2D Grid quadratic") { } } -TEST_CASE("2D Grid non-quadratic") { +TEST_CASE("2D Grid64 non-quadratic") { int r = 12; int c = 15; - Grid grid(r, c); + Grid64 grid(r, c); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 2); diff --git a/test/testSimulation.cpp b/test/testSimulation.cpp index 004799a..e1b0cf3 100644 --- a/test/testSimulation.cpp +++ b/test/testSimulation.cpp @@ -11,7 +11,7 @@ using namespace Eigen; using namespace std; -static Grid setupSimulation(APPROACH approach, double timestep, +static Grid64 setupSimulation(APPROACH approach, double timestep, int iterations) { int row = 11; int col = 11; @@ -19,7 +19,7 @@ static Grid setupSimulation(APPROACH approach, double timestep, int domain_col = 10; // Grid - Grid grid = Grid(row, col); + Grid grid = Grid64(row, col); grid.setDomain(domain_row, domain_col); MatrixXd concentrations = MatrixXd::Constant(row, col, 0); @@ -80,7 +80,7 @@ TEST_CASE("equality to reference matrix with BTCS") { TEST_CASE("Initialize environment") { int rc = 12; - Grid grid(rc, rc); + Grid64 grid(rc, rc); Boundary boundary(grid); CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH)); @@ -88,7 +88,7 @@ TEST_CASE("Initialize environment") { TEST_CASE("Simulation environment") { int rc = 12; - Grid grid(rc, rc); + Grid64 grid(rc, rc); Boundary boundary(grid); Simulation sim(grid, boundary, FTCS_APPROACH);