diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index b7ed776..95421b7 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -12,7 +12,7 @@ test: stage: test script: - mkdir build && cd build - - cmake -DCMAKE_BUILD_TYPE=Release -DTUG_ENABLE_TESTING=ON -G Ninja .. + - cmake -DCMAKE_BUILD_TYPE=Debug -DTUG_ENABLE_TESTING=ON -G Ninja .. - ninja - ctest --output-junit test_results.xml artifacts: @@ -27,7 +27,7 @@ pages: image: python:slim before_script: - apt-get update && apt-get install --no-install-recommends -y graphviz imagemagick doxygen make - - pip install --upgrade pip && pip install Sphinx Pillow breathe sphinx-rtd-theme m2r2 + - pip install --upgrade pip && pip install Sphinx Pillow breathe sphinx-rtd-theme sphinx-mdinclude - mkdir public script: - pushd docs_sphinx diff --git a/docs_sphinx/Simulation.rst b/docs_sphinx/Diffusion.rst similarity index 79% rename from docs_sphinx/Simulation.rst rename to docs_sphinx/Diffusion.rst index a461073..1ab1cc1 100644 --- a/docs_sphinx/Simulation.rst +++ b/docs_sphinx/Diffusion.rst @@ -1,4 +1,4 @@ -Simulation +Diffusion ========== .. doxygenenum:: tug::APPROACH @@ -7,4 +7,4 @@ Simulation .. doxygenenum:: tug::CONSOLE_OUTPUT .. doxygenenum:: tug::TIME_MEASURE -.. doxygenclass:: tug::Simulation +.. doxygenclass:: tug::Diffusion diff --git a/docs_sphinx/Grid.rst b/docs_sphinx/Grid.rst deleted file mode 100644 index e4dc497..0000000 --- a/docs_sphinx/Grid.rst +++ /dev/null @@ -1,4 +0,0 @@ -Grid -==== - -.. doxygenclass:: tug::Grid diff --git a/docs_sphinx/conf.py b/docs_sphinx/conf.py index 9087a8a..75137c9 100644 --- a/docs_sphinx/conf.py +++ b/docs_sphinx/conf.py @@ -42,7 +42,7 @@ extensions = [ 'sphinx.ext.viewcode', 'sphinx.ext.inheritance_diagram', 'breathe', - 'm2r2' + 'sphinx_mdinclude' ] html_baseurl = "/index.html" diff --git a/docs_sphinx/user.rst b/docs_sphinx/user.rst index ddd4941..84352ea 100644 --- a/docs_sphinx/user.rst +++ b/docs_sphinx/user.rst @@ -3,6 +3,5 @@ User API .. toctree:: - Grid Boundary - Simulation + Diffusion diff --git a/examples/BTCS_1D_proto_example.cpp b/examples/BTCS_1D_proto_example.cpp deleted file mode 100644 index af9888e..0000000 --- a/examples/BTCS_1D_proto_example.cpp +++ /dev/null @@ -1,51 +0,0 @@ -#include -#include - -using namespace Eigen; -using namespace tug; - -int main(int argc, char *argv[]) { - // ************** - // **** GRID **** - // ************** - - // create a linear grid with 20 cells - int cells = 20; - Grid64 grid(cells); - - MatrixXd concentrations = MatrixXd::Constant(1, 20, 0); - concentrations(0, 0) = 2000; - // TODO add option to set concentrations with a vector in 1D case - grid.setConcentrations(concentrations); - - // ****************** - // **** BOUNDARY **** - // ****************** - - // create a boundary with constant values - Boundary bc = Boundary(grid); - bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); - bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); - - // ************************ - // **** SIMULATION ENV **** - // ************************ - - // set up a simulation environment - Simulation simulation = Simulation(grid, bc); // grid,boundary - - // set the timestep of the simulation - simulation.setTimestep(0.1); // timestep - - // set the number of iterations - simulation.setIterations(100); - - // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, - // CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - - // **** RUN SIMULATION **** - - // run the simulation - simulation.run(); -} diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index 98acdb9..c09d93d 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -1,5 +1,5 @@ #include -#include +#include using namespace Eigen; using namespace tug; @@ -14,7 +14,6 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int row = 40; int col = 50; - Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); @@ -24,7 +23,7 @@ int main(int argc, char *argv[]) { // #row,#col,value grid.setConcentrations(concentrations); MatrixXd concentrations = MatrixXd::Constant(row, col, 0); concentrations(10, 10) = 2000; - grid.setConcentrations(concentrations); + Grid64 grid(concentrations); // (optional) set alphas of the grid, e.g.: // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value @@ -61,8 +60,7 @@ int main(int argc, char *argv[]) { // ************************ // set up a simulation environment - Simulation simulation = - Simulation(grid, bc); // grid,boundary + Diffusion simulation(grid, bc); // grid,boundary // set the timestep of the simulation simulation.setTimestep(0.1); // timestep diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 0a3173d..3c8bf28 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,20 +1,7 @@ -add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp) -add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp) -add_executable(BTCS_1D_proto_example BTCS_1D_proto_example.cpp) add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp) -add_executable(CRNI_2D_proto_example CRNI_2D_proto_example.cpp) -add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp) -add_executable(profiling_openmp profiling_openmp.cpp) - -target_link_libraries(FTCS_1D_proto_example tug) -target_link_libraries(FTCS_2D_proto_example tug) -target_link_libraries(BTCS_1D_proto_example tug) -target_link_libraries(BTCS_2D_proto_example tug) -target_link_libraries(CRNI_2D_proto_example tug) -target_link_libraries(reference-FTCS_2D_closed tug) -target_link_libraries(profiling_openmp tug) - add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp) add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp) + +target_link_libraries(BTCS_2D_proto_example tug) target_link_libraries(FTCS_2D_proto_closed_mdl tug) target_link_libraries(FTCS_2D_proto_example_mdl tug) diff --git a/examples/CRNI_2D_proto_example.cpp b/examples/CRNI_2D_proto_example.cpp deleted file mode 100644 index fc3f046..0000000 --- a/examples/CRNI_2D_proto_example.cpp +++ /dev/null @@ -1,29 +0,0 @@ -#include -#include - -using namespace Eigen; -using namespace tug; - -int main(int argc, char *argv[]) { - int row = 20; - int col = 20; - Grid64 grid(row, col); - - MatrixXd concentrations = MatrixXd::Constant(row, col, 0); - concentrations(10, 10) = 2000; - grid.setConcentrations(concentrations); - - Boundary bc = Boundary(grid); - bc.setBoundarySideClosed(BC_SIDE_LEFT); - bc.setBoundarySideClosed(BC_SIDE_RIGHT); - bc.setBoundarySideClosed(BC_SIDE_TOP); - bc.setBoundarySideClosed(BC_SIDE_BOTTOM); - - Simulation simulation = - Simulation(grid, bc); - simulation.setTimestep(0.1); - simulation.setIterations(50); - simulation.setOutputCSV(CSV_OUTPUT_XTREME); - - simulation.run(); -} diff --git a/examples/FTCS_1D_proto_example.cpp b/examples/FTCS_1D_proto_example.cpp deleted file mode 100644 index b3dc905..0000000 --- a/examples/FTCS_1D_proto_example.cpp +++ /dev/null @@ -1,51 +0,0 @@ -#include "tug/Boundary.hpp" -#include -#include - -using namespace Eigen; -using namespace tug; - -int main(int argc, char *argv[]) { - // ************** - // **** GRID **** - // ************** - - // create a linear grid with 20 cells - int cells = 20; - Grid64 grid(cells); - - MatrixXd concentrations = MatrixXd::Constant(1, 20, 20); - grid.setConcentrations(concentrations); - - // ****************** - // **** BOUNDARY **** - // ****************** - - // create a boundary with constant values - Boundary bc = Boundary(grid); - bc.setBoundarySideConstant(BC_SIDE_LEFT, 1); - bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1); - - // ************************ - // **** SIMULATION ENV **** - // ************************ - - // set up a simulation environment - Simulation simulation = - Simulation(grid, bc); // grid,boundary,simulation-approach - - // (optional) set the timestep of the simulation - simulation.setTimestep(0.1); // timestep - - // (optional) set the number of iterations - simulation.setIterations(100); - - // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, - // CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_OFF); - - // **** RUN SIMULATION **** - - // run the simulation - simulation.run(); -} diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index c6ed02f..69b31c3 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -9,7 +9,7 @@ #include #include #include -#include +#include using namespace Eigen; using namespace tug; @@ -31,7 +31,6 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int n2 = row / 2 - 1; - Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); @@ -44,7 +43,7 @@ int main(int argc, char *argv[]) { concentrations(n2, n2 + 1) = 1; concentrations(n2 + 1, n2) = 1; concentrations(n2 + 1, n2 + 1) = 1; - grid.setConcentrations(concentrations); + Grid64 grid(concentrations); // (optional) set alphas of the grid, e.g.: MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value @@ -69,8 +68,8 @@ int main(int argc, char *argv[]) { // ************************ // set up a simulation environment - Simulation simulation = - Simulation(grid, bc); // grid,boundary,simulation-approach + Diffusion simulation( + grid, bc); // grid,boundary,simulation-approach // set the timestep of the simulation simulation.setTimestep(10000); // timestep diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp deleted file mode 100644 index 9b8987a..0000000 --- a/examples/FTCS_2D_proto_example.cpp +++ /dev/null @@ -1,92 +0,0 @@ -/** - * @file FTCS_2D_proto_example.cpp - * @author Hannes Signer, Philipp Ungrund - * @brief Creates a prototypical standard TUG simulation in 2D with FTCS - * approach and constant boundary condition - * - */ - -#include -#include - -using namespace Eigen; -using namespace tug; -// #include -// #define EASY_PROFILER_ENABLE ::profiler::setEnabled(true); - -int main(int argc, char *argv[]) { - // EASY_PROFILER_ENABLE; - // profiler::startListen(); - // ************** - // **** GRID **** - // ************** - // profiler::startListen(); - // create a grid with a 20 x 20 field - int row = 20; - int col = 20; - Grid64 grid(row, col); - - // (optional) set the domain, e.g.: - // grid.setDomain(20, 20); - - // (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, 0); - concentrations(0, 0) = 1999; - grid.setConcentrations(concentrations); - - // (optional) set alphas of the grid, e.g.: - // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value - // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value - // grid.setAlpha(alphax, alphay); - - // ****************** - // **** BOUNDARY **** - // ****************** - - // create a boundary with constant values - 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 - // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values - // VectorXd bc_zero_values = VectorXd::Constant(20,0); - // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); - // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); - // VectorXd bc_front_values = VectorXd::Constant(20,2000); - // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); - // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); - - // ************************ - // **** SIMULATION ENV **** - // ************************ - - // set up a simulation environment - Simulation simulation = - Simulation(grid, bc); // grid,boundary,simulation-approach - - // set the timestep of the simulation - simulation.setTimestep(0.1); // timestep - - // set the number of iterations - simulation.setIterations(10000); - - // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, - // CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - - // **** RUN SIMULATION **** - - // run the simulation - - // EASY_BLOCK("SIMULATION") - simulation.run(); - // EASY_END_BLOCK; - // profiler::dumpBlocksToFile("test_profile.prof"); - // profiler::stopListen(); -} diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp index fdf40d1..41ddb91 100644 --- a/examples/FTCS_2D_proto_example_mdl.cpp +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -7,7 +7,7 @@ */ #include -#include +#include using namespace Eigen; using namespace tug; @@ -22,7 +22,6 @@ int main(int argc, char *argv[]) { int row = 64; int col = 64; int n2 = row / 2 - 1; - Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); @@ -35,7 +34,7 @@ int main(int argc, char *argv[]) { concentrations(n2, n2 + 1) = 1; concentrations(n2 + 1, n2) = 1; concentrations(n2 + 1, n2 + 1) = 1; - grid.setConcentrations(concentrations); + Grid64 grid(concentrations); // (optional) set alphas of the grid, e.g.: MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value @@ -60,8 +59,8 @@ int main(int argc, char *argv[]) { // ************************ // set up a simulation environment - Simulation simulation = - Simulation(grid, bc); // grid,boundary,simulation-approach + Diffusion simulation( + grid, bc); // grid,boundary,simulation-approach // (optional) set the timestep of the simulation simulation.setTimestep(1000); // timestep diff --git a/examples/profiling_openmp.cpp b/examples/profiling_openmp.cpp deleted file mode 100644 index 8ad62e0..0000000 --- a/examples/profiling_openmp.cpp +++ /dev/null @@ -1,70 +0,0 @@ -#include -#include -#include -#include -#include -#include - -using namespace Eigen; -using namespace std; -using namespace tug; - -int main(int argc, char *argv[]) { - - int n[] = {2000}; - int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; - int iterations[1] = {1}; - int repetition = 10; - - for (int l = 0; l < size(threads); l++) { - // string filename = "ftcs_openmp_" + to_string(threads[l]) + ".csv"; - ofstream myfile; - myfile.open("speedup_1000.csv", std::ios::app); - myfile << "Number threads: " << threads[l] << endl; - - for (int i = 0; i < size(n); i++) { - cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl; - // myfile << "Grid size: " << n[i] << " x " << n[i] << endl << endl; - for (int j = 0; j < size(iterations); j++) { - cout << "Iterations: " << iterations[j] << endl; - // myfile << "Iterations: " << iterations[j] << endl; - for (int k = 0; k < repetition; k++) { - cout << "Wiederholung: " << k << endl; - Grid64 grid(n[i], n[i]); - grid.setDomain(1, 1); - - MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); - concentrations(n[i] / 2, n[i] / 2) = 1; - grid.setConcentrations(concentrations); - MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5); - - Boundary bc = Boundary(grid); - - Simulation sim = Simulation(grid, bc); - - if (argc == 2) { - int numThreads = atoi(argv[1]); - sim.setNumberThreads(numThreads); - } else { - sim.setNumberThreads(threads[l]); - } - - sim.setTimestep(0.01); - sim.setIterations(iterations[j]); - sim.setOutputCSV(CSV_OUTPUT_OFF); - - auto begin = std::chrono::high_resolution_clock::now(); - sim.run(); - auto end = std::chrono::high_resolution_clock::now(); - auto milliseconds = - std::chrono::duration_cast(end - - begin); - myfile << milliseconds.count() << endl; - } - } - cout << endl; - myfile << endl; - } - myfile.close(); - } -} diff --git a/examples/profiling_speedup.cpp b/examples/profiling_speedup.cpp deleted file mode 100644 index 8ad62e0..0000000 --- a/examples/profiling_speedup.cpp +++ /dev/null @@ -1,70 +0,0 @@ -#include -#include -#include -#include -#include -#include - -using namespace Eigen; -using namespace std; -using namespace tug; - -int main(int argc, char *argv[]) { - - int n[] = {2000}; - int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; - int iterations[1] = {1}; - int repetition = 10; - - for (int l = 0; l < size(threads); l++) { - // string filename = "ftcs_openmp_" + to_string(threads[l]) + ".csv"; - ofstream myfile; - myfile.open("speedup_1000.csv", std::ios::app); - myfile << "Number threads: " << threads[l] << endl; - - for (int i = 0; i < size(n); i++) { - cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl; - // myfile << "Grid size: " << n[i] << " x " << n[i] << endl << endl; - for (int j = 0; j < size(iterations); j++) { - cout << "Iterations: " << iterations[j] << endl; - // myfile << "Iterations: " << iterations[j] << endl; - for (int k = 0; k < repetition; k++) { - cout << "Wiederholung: " << k << endl; - Grid64 grid(n[i], n[i]); - grid.setDomain(1, 1); - - MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); - concentrations(n[i] / 2, n[i] / 2) = 1; - grid.setConcentrations(concentrations); - MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5); - - Boundary bc = Boundary(grid); - - Simulation sim = Simulation(grid, bc); - - if (argc == 2) { - int numThreads = atoi(argv[1]); - sim.setNumberThreads(numThreads); - } else { - sim.setNumberThreads(threads[l]); - } - - sim.setTimestep(0.01); - sim.setIterations(iterations[j]); - sim.setOutputCSV(CSV_OUTPUT_OFF); - - auto begin = std::chrono::high_resolution_clock::now(); - sim.run(); - auto end = std::chrono::high_resolution_clock::now(); - auto milliseconds = - std::chrono::duration_cast(end - - begin); - myfile << milliseconds.count() << endl; - } - } - cout << endl; - myfile << endl; - } - myfile.close(); - } -} diff --git a/examples/reference-FTCS_2D_closed.cpp b/examples/reference-FTCS_2D_closed.cpp deleted file mode 100644 index f08c689..0000000 --- a/examples/reference-FTCS_2D_closed.cpp +++ /dev/null @@ -1,53 +0,0 @@ -#include "Eigen/Core" -#include -#include - -using namespace std; -using namespace Eigen; -using namespace tug; - -int main(int argc, char *argv[]) { - int row = 50; - int col = 50; - int domain_row = 10; - int domain_col = 10; - - // Grid - Grid64 grid(row, col); - grid.setDomain(domain_row, domain_col); - - MatrixXd concentrations = MatrixXd::Constant(row, col, 0); - concentrations(5, 5) = 1; - grid.setConcentrations(concentrations); - - MatrixXd alpha = MatrixXd::Constant(row, col, 1); - for (int i = 0; i < 5; i++) { - for (int j = 0; j < 6; j++) { - alpha(i, j) = 0.01; - } - } - for (int i = 0; i < 5; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.001; - } - } - for (int i = 5; i < 11; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.1; - } - } - grid.setAlpha(alpha, alpha); - - // Boundary - Boundary bc = Boundary(grid); - - // Simulation - Simulation sim = Simulation(grid, bc); - sim.setTimestep(0.001); - sim.setIterations(10000); - sim.setOutputCSV(CSV_OUTPUT_OFF); - sim.setOutputConsole(CONSOLE_OUTPUT_OFF); - - // RUN - sim.run(); -} diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index ebcde98..114edad 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -7,11 +7,13 @@ #ifndef BOUNDARY_H_ #define BOUNDARY_H_ -#include "Grid.hpp" +#include "tug/Core/TugUtils.hpp" +#include #include #include #include +#include #include #include @@ -44,7 +46,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. @@ -114,7 +116,7 @@ public: * * @param length Length of the grid */ - Boundary(std::uint32_t length) : Boundary(Grid(length)){}; + Boundary(std::uint32_t length) : Boundary(1, length) {}; /** * @brief Creates a boundary object for a 2D grid @@ -123,17 +125,7 @@ public: * @param cols Number of columns of the grid */ Boundary(std::uint32_t rows, std::uint32_t cols) - : Boundary(Grid(rows, cols)){}; - - /** - * @brief Creates a boundary object based on the passed grid object and - * initializes the boundaries as closed. - * - * @param grid Grid object on the basis of which the simulation takes place - * and from which the dimensions (in 2D case) are taken. - */ - Boundary(const Grid &grid) - : dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) { + : dim(rows == 1 ? 1 : 2), cols(cols), rows(rows) { if (this->dim == 1) { this->boundaries = std::vector>>( 2); // in 1D only left and right boundary @@ -152,8 +144,37 @@ public: this->boundaries[BC_SIDE_BOTTOM] = std::vector>(this->cols, BoundaryElement()); } - } + }; + /** + * @brief Creates a boundary object based on the passed grid object and + * initializes the boundaries as closed. + * + * @param grid Grid object on the basis of which the simulation takes place + * and from which the dimensions (in 2D case) are taken. + */ + // 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 * condition closed. @@ -161,13 +182,11 @@ public: * @param side Side to be set to closed, e.g. BC_SIDE_LEFT. */ 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."); - } - } + tug_assert((this->dim > 1) || + ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)), + "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; @@ -186,13 +205,11 @@ public: * page. */ 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."); - } - } + tug_assert((this->dim > 1) || + ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)), + "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; @@ -212,10 +229,9 @@ public: */ 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."); - } + tug_assert(boundaries[side].size() > index && index >= 0, + "Index is selected either too large or too small."); + this->boundaries[side][index].setType(BC_TYPE_CLOSED); } @@ -233,10 +249,8 @@ public: */ 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."); - } + tug_assert(boundaries[side].size() > index && index >= 0, + "Index is selected either too large or too small."); this->boundaries[side][index].setType(BC_TYPE_CONSTANT); this->boundaries[side][index].setValue(value); } @@ -251,13 +265,11 @@ public: * BoundaryElement objects. */ 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."); - } - } + tug_assert((this->dim > 1) || + ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)), + "For the " + "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT " + "borders exist."); return this->boundaries[side]; } @@ -295,10 +307,8 @@ public: * object. */ 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."); - } + tug_assert(boundaries[side].size() > index && index >= 0, + "Index is selected either too large or too small."); return this->boundaries[side][index]; } @@ -313,10 +323,8 @@ public: * @return Boundary Type of the corresponding boundary condition. */ 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."); - } + tug_assert(boundaries[side].size() > index && index >= 0, + "Index is selected either too large or too small."); return this->boundaries[side][index].getType(); } @@ -333,14 +341,12 @@ public: * object. */ 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."); - } + tug_assert(boundaries[side].size() > index && index >= 0, + "Index is selected either too large or too small."); + tug_assert( + boundaries[side][index].getType() == BC_TYPE_CONSTANT, + "A value can only be output if it is a constant boundary condition."); + return this->boundaries[side][index].getValue(); } @@ -382,13 +388,8 @@ public: * @param value Value of the inner constant boundary condition */ void setInnerBoundary(std::uint32_t index, T value) { - if (this->dim != 1) { - throw std::invalid_argument( - "This function is only available for 1D grids."); - } - if (index >= this->cols) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(this->dim == 1, "This function is only available for 1D grids."); + tug_assert(index < this->cols, "Index is out of bounds."); this->inner_boundary[std::make_pair(0, index)] = value; } @@ -401,13 +402,8 @@ public: * @param value Value of the inner constant boundary condition */ void setInnerBoundary(std::uint32_t row, std::uint32_t col, T value) { - if (this->dim != 2) { - throw std::invalid_argument( - "This function is only available for 2D grids."); - } - if (row >= this->rows || col >= this->cols) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(this->dim == 2, "This function is only available for 2D grids."); + tug_assert(row < this->rows && col < this->cols, "Index is out of bounds."); this->inner_boundary[std::make_pair(row, col)] = value; } @@ -420,13 +416,8 @@ public: * set or not) and value of the inner constant boundary condition */ std::pair getInnerBoundary(std::uint32_t index) const { - if (this->dim != 1) { - throw std::invalid_argument( - "This function is only available for 1D grids."); - } - if (index >= this->cols) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(this->dim == 1, "This function is only available for 1D grids."); + tug_assert(index < this->cols, "Index is out of bounds."); auto it = this->inner_boundary.find(std::make_pair(0, index)); if (it == this->inner_boundary.end()) { @@ -445,13 +436,8 @@ public: */ std::pair getInnerBoundary(std::uint32_t row, std::uint32_t col) const { - if (this->dim != 2) { - throw std::invalid_argument( - "This function is only available for 2D grids."); - } - if (row >= this->rows || col >= this->cols) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(this->dim == 2, "This function is only available for 2D grids."); + tug_assert(row < this->rows && col < this->cols, "Index is out of bounds."); auto it = this->inner_boundary.find(std::make_pair(row, col)); if (it == this->inner_boundary.end()) { @@ -471,9 +457,7 @@ public: * condition */ std::vector> getInnerBoundaryRow(std::uint32_t row) const { - if (row >= this->rows) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(row < this->rows, "Index is out of bounds."); if (inner_boundary.empty()) { return std::vector>(this->cols, @@ -499,14 +483,8 @@ public: * condition */ std::vector> getInnerBoundaryCol(std::uint32_t col) const { - if (this->dim != 2) { - throw std::invalid_argument( - "This function is only available for 2D grids."); - } - - if (col >= this->cols) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(this->dim == 2, "This function is only available for 2D grids."); + tug_assert(col < this->cols, "Index is out of bounds."); if (inner_boundary.empty()) { return std::vector>(this->rows, diff --git a/include/tug/Core/BaseSimulation.hpp b/include/tug/Core/BaseSimulation.hpp new file mode 100644 index 0000000..82a1480 --- /dev/null +++ b/include/tug/Core/BaseSimulation.hpp @@ -0,0 +1,388 @@ +#pragma once + +#include "tug/Boundary.hpp" +#include +#include +#include +#include + +namespace tug { + +/** + * @brief Enum holding different options for .csv output. + * + */ +enum class CSV_OUTPUT { + OFF, /*!< do not produce csv output */ + ON, /*!< produce csv output with last concentration matrix */ + VERBOSE, /*!< produce csv output with all concentration matrices */ + XTREME /*!< csv output like VERBOSE but additional boundary + conditions at beginning */ +}; + +/** + * @brief Enum holding different options for console output. + * + */ +enum class CONSOLE_OUTPUT { + OFF, /*!< do not print any output to console */ + ON, /*!< print before and after concentrations to console */ + VERBOSE /*!< print all concentration matrices to console */ +}; + +/** + * @brief Enum holding different options for time measurement. + * + */ +enum class TIME_MEASURE { + OFF, /*!< do not print any time measures */ + ON /*!< print time measure after last iteration */ +}; + +/** + * @brief A base class for simulation grids. + * + * This class provides a base implementation for simulation grids, including + * methods for setting and getting grid dimensions, domain sizes, and output + * options. It also includes methods for running simulations, which must be + * implemented by derived classes. + * + * @tparam T The type of the elements in the grid. + */ +template class BaseSimulationGrid { +private: + CSV_OUTPUT csv_output{CSV_OUTPUT::OFF}; + CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT::OFF}; + TIME_MEASURE time_measure{TIME_MEASURE::OFF}; + + int iterations{1}; + RowMajMatMap concentrationMatrix; + Boundary boundaryConditions; + + const std::uint8_t dim; + + T delta_col; + T delta_row; + +public: + /** + * @brief Constructs a BaseSimulationGrid from a given RowMajMat object. + * + * This constructor initializes a BaseSimulationGrid using the data, number of + * rows, and number of columns from the provided RowMajMat object. + * + * @tparam T The type of the elements in the RowMajMat. + * @param origin The RowMajMat object from which to initialize the + * BaseSimulationGrid. + */ + BaseSimulationGrid(RowMajMat &origin) + : BaseSimulationGrid(origin.data(), origin.rows(), origin.cols()) {} + + /** + * @brief Constructs a BaseSimulationGrid object. + * + * @tparam T The type of the data elements. + * @param data Pointer to the data array. + * @param rows Number of rows in the grid. + * @param cols Number of columns in the grid. + * + * Initializes the concentration_matrix with the provided data, rows, and + * columns. Sets delta_col and delta_row to 1. Determines the dimension (dim) + * based on the number of rows: if rows == 1, dim is set to 1; otherwise, dim + * is set to 2. + */ + BaseSimulationGrid(T *data, std::size_t rows, std::size_t cols) + : concentrationMatrix(data, rows, cols), boundaryConditions(rows, cols), + delta_col(1), delta_row(1), dim(rows == 1 ? 1 : 2) {} + + /** + * @brief Constructs a BaseSimulationGrid with a single dimension. + * + * This constructor initializes a BaseSimulationGrid object with the provided + * data and length. It assumes the grid has only one dimension. + * + * @param data Pointer to the data array. + * @param length The length of the data array. + */ + BaseSimulationGrid(T *data, std::size_t length) + : BaseSimulationGrid(data, 1, length) {} + + /** + * @brief Overloaded function call operator to access elements in a + * one-dimensional grid. + * + * This operator provides access to elements in the concentration matrix using + * a single index. It asserts that the grid is one-dimensional before + * accessing the element. + * + * @tparam T The type of elements in the concentration matrix. + * @param index The index of the element to access. + * @return A reference to the element at the specified index in the + * concentration matrix. + */ + constexpr T &operator()(std::size_t index) { + tug_assert(dim == 1, "Grid is not one dimensional, use 2D index operator!"); + + return concentrationMatrix(index); + } + + /** + * @brief Overloaded function call operator to access elements in a 2D + * concentration matrix. + * + * This operator allows accessing elements in the concentration matrix using + * row and column indices. It asserts that the grid is two-dimensional before + * accessing the element. + * + * @param row The row index of the element to access. + * @param col The column index of the element to access. + * @return A reference to the element at the specified row and column in the + * concentration matrix. + */ + constexpr T &operator()(std::size_t row, std::size_t col) { + tug_assert(dim == 2, "Grid is not two dimensional, use 1D index operator!"); + + return concentrationMatrix(row, col); + } + + /** + * @brief Retrieves the concentration matrix. + * + * @tparam T The data type of the elements in the concentration matrix. + * @return RowMajMat& Reference to the concentration matrix. + */ + RowMajMatMap &getConcentrationMatrix() { return concentrationMatrix; } + + const RowMajMatMap &getConcentrationMatrix() const { + return concentrationMatrix; + } + + /** + * @brief Retrieves the boundary conditions for the simulation. + * + * @tparam T The type parameter for the Boundary class. + * @return Boundary& A reference to the boundary conditions. + */ + Boundary &getBoundaryConditions() { return boundaryConditions; } + + const Boundary &getBoundaryConditions() const { + return boundaryConditions; + } + + /** + * @brief Retrieves the dimension value. + * + * @return The dimension value as an 8-bit unsigned integer. + */ + std::uint8_t getDim() const { return dim; } + + /** + * @brief Returns the number of rows in the concentration matrix. + * + * @return std::size_t The number of rows in the concentration matrix. + */ + std::size_t rows() const { return concentrationMatrix.rows(); } + + /** + * @brief Get the number of columns in the concentration matrix. + * + * @return std::size_t The number of columns in the concentration matrix. + */ + std::size_t cols() const { return concentrationMatrix.cols(); } + + /** + * @brief Returns the cell size in meter of the x-direction. + * + * This function returns the value of the delta column, which is used + * to represent the difference or change in the column value. + * + * @return T The cell size in meter of the x-direction. + */ + T deltaCol() const { return delta_col; } + + /** + * @brief Returns the cell size in meter of the y-direction. + * + * This function asserts that the grid is two-dimensional. If the grid is not + * two-dimensional, an assertion error is raised with the message "Grid is not + * two dimensional, there is no delta in y-direction!". + * + * @return The cell size in meter of the y-direction. + */ + T deltaRow() const { + tug_assert( + dim == 2, + "Grid is not two dimensional, there is no delta in y-direction!"); + + return delta_row; + } + + /** + * @brief Computes the domain size in the X direction. + * + * This function calculates the size of the domain in the X direction by + * multiplying the column spacing (delta_col) by the number of columns (cols). + * + * @return The size of the domain in the X direction. + */ + T domainX() const { return delta_col * cols(); } + + /** + * @brief Returns the size of the domain in the y-direction. + * + * This function calculates the size of the domain in the y-direction + * by multiplying the row spacing (delta_row) by the number of rows. + * It asserts that the grid is two-dimensional before performing the + * calculation. + * + * @return The size of the domain in the y-direction. + */ + T domainY() const { + tug_assert( + dim == 2, + "Grid is not two dimensional, there is no domain in y-direction!"); + + return delta_row * rows(); + } + + /** + * @brief Sets the domain length for a one-dimensional grid. + * + * This function sets the domain length for a one-dimensional grid and + * calculates the column width (delta_col) based on the given domain length + * and the number of columns. It asserts that the grid is one-dimensional and + * that the given domain length is positive. + * + * @param domain_length The length of the domain. Must be positive. + */ + void setDomain(T domain_length) { + tug_assert(dim == 1, "Grid is not one dimensional, use 2D domain setter!"); + tug_assert(domain_length > 0, "Given domain length is not positive!"); + + delta_col = domain_length / cols(); + } + + /** + * @brief Sets the domain size for a 2D grid simulation. + * + * This function sets the domain size in the x and y directions for a + * two-dimensional grid simulation. It asserts that the grid is indeed + * two-dimensional and that the provided domain sizes are positive. + * + * @tparam T The type of the domain size parameters. + * @param domain_row The size of the domain in the y-direction. + * @param domain_col The size of the domain in the x-direction. + */ + void setDomain(T domain_row, T domain_col) { + tug_assert(dim == 2, "Grid is not two dimensional, use 1D domain setter!"); + tug_assert(domain_col > 0, + "Given domain size in x-direction is not positive!"); + tug_assert(domain_row > 0, + "Given domain size in y-direction is not positive!"); + + delta_row = domain_row / rows(); + delta_col = domain_col / cols(); + } + + /** + * @brief Set the option to output the results to a CSV file. Off by default. + * + * + * @param csv_output Valid output option. The following options can be set + * here: + * - CSV_OUTPUT_OFF: do not produce csv output + * - CSV_OUTPUT_ON: produce csv output with last + * concentration matrix + * - CSV_OUTPUT_VERBOSE: produce csv output with all + * concentration matrices + * - CSV_OUTPUT_XTREME: produce csv output with all + * concentration matrices and simulation environment + */ + void setOutputCSV(CSV_OUTPUT csv_output) { this->csv_output = csv_output; } + + /** + * @brief Retrieves the CSV output. + * + * This function returns the CSV output associated with the simulation. + * + * @return CSV_OUTPUT The CSV output of the simulation. + */ + constexpr CSV_OUTPUT getOutputCSV() const { return this->csv_output; } + + /** + * @brief Set the options for outputting information to the console. Off by + * default. + * + * @param console_output Valid output option. The following options can be set + * here: + * - CONSOLE_OUTPUT_OFF: do not print any output to + * console + * - CONSOLE_OUTPUT_ON: print before and after + * concentrations to console + * - CONSOLE_OUTPUT_VERBOSE: print all concentration + * matrices to console + */ + void setOutputConsole(CONSOLE_OUTPUT console_output) { + this->console_output = console_output; + } + + /** + * @brief Retrieves the console output. + * + * This function returns the current state of the console output. + * + * @return CONSOLE_OUTPUT The current console output. + */ + constexpr CONSOLE_OUTPUT getOutputConsole() const { + return this->console_output; + } + + /** + * @brief Set the Time Measure option. Off by default. + * + * @param time_measure The following options are allowed: + * - TIME_MEASURE_OFF: Time of simulation is not printed + * to console + * - TIME_MEASURE_ON: Time of simulation run is printed to + * console + */ + void setTimeMeasure(TIME_MEASURE time_measure) { + this->time_measure = time_measure; + } + + /** + * @brief Retrieves the current time measurement. + * + * @return TIME_MEASURE The current time measurement. + */ + constexpr TIME_MEASURE getTimeMeasure() const { return this->time_measure; } + + /** + * @brief Set the desired iterations to be calculated. A value greater + * than zero must be specified here. Setting iterations is required. + * + * @param iterations Number of iterations to be simulated. + */ + void setIterations(int iterations) { + tug_assert(iterations > 0, + "Number of iterations must be greater than zero."); + + this->iterations = iterations; + } + + /** + * @brief Return the currently set iterations to be calculated. + * + * @return int Number of iterations. + */ + int getIterations() const { return this->iterations; } + + /** + * @brief Method starts the simulation process with the previously set + * parameters. + */ + virtual void run() = 0; + + virtual void setTimestep(T timestep) = 0; +}; +} // namespace tug diff --git a/include/tug/Core/FTCS.hpp b/include/tug/Core/FTCS.hpp deleted file mode 100644 index a3012b3..0000000 --- a/include/tug/Core/FTCS.hpp +++ /dev/null @@ -1,402 +0,0 @@ -/** - * @file FTCS.hpp - * @brief Implementation of heterogenous FTCS (forward time-centered space) - * solution of diffusion equation in 1D and 2D space. - * - */ - -#ifndef FTCS_H_ -#define FTCS_H_ - -#include "TugUtils.hpp" - -#include -#include -#include - -#ifdef _OPENMP -#include -#else -#define omp_get_thread_num() 0 -#endif - -namespace tug { - -// calculates horizontal change on one cell independent of boundary type -template -static inline T calcHorizontalChange(Grid &grid, int &row, int &col) { - - return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), - grid.getAlphaX()(row, col)) * - grid.getConcentrations()(row, col + 1) - - (calcAlphaIntercell(grid.getAlphaX()(row, col + 1), - grid.getAlphaX()(row, col)) + - calcAlphaIntercell(grid.getAlphaX()(row, col - 1), - grid.getAlphaX()(row, col))) * - grid.getConcentrations()(row, col) + - calcAlphaIntercell(grid.getAlphaX()(row, col - 1), - grid.getAlphaX()(row, col)) * - grid.getConcentrations()(row, col - 1); -} - -// calculates vertical change on one cell independent of boundary type -template -static inline T calcVerticalChange(Grid &grid, int &row, int &col) { - - return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), - grid.getAlphaY()(row, col)) * - grid.getConcentrations()(row + 1, col) - - (calcAlphaIntercell(grid.getAlphaY()(row + 1, col), - grid.getAlphaY()(row, col)) + - calcAlphaIntercell(grid.getAlphaY()(row - 1, col), - grid.getAlphaY()(row, col))) * - grid.getConcentrations()(row, col) + - calcAlphaIntercell(grid.getAlphaY()(row - 1, col), - grid.getAlphaY()(row, col)) * - grid.getConcentrations()(row - 1, col); -} - -// calculates horizontal change on one cell with a constant left boundary -template -static inline T calcHorizontalChangeLeftBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, int &col) { - - return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), - grid.getAlphaX()(row, col)) * - grid.getConcentrations()(row, col + 1) - - (calcAlphaIntercell(grid.getAlphaX()(row, col + 1), - grid.getAlphaX()(row, col)) + - 2 * grid.getAlphaX()(row, col)) * - grid.getConcentrations()(row, col) + - 2 * grid.getAlphaX()(row, col) * - bc.getBoundaryElementValue(BC_SIDE_LEFT, row); -} - -// calculates horizontal change on one cell with a closed left boundary -template -static inline T calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, - int &col) { - - return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), - grid.getAlphaX()(row, col)) * - (grid.getConcentrations()(row, col + 1) - - grid.getConcentrations()(row, col)); -} - -// checks boundary condition type for a cell on the left edge of grid -template -static inline T calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { - if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CONSTANT) { - return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col); - } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CLOSED) { - return calcHorizontalChangeLeftBoundaryClosed(grid, row, col); - } else { - throw_invalid_argument("Undefined Boundary Condition Type!"); - } -} - -// calculates horizontal change on one cell with a constant right boundary -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) - - (calcAlphaIntercell(grid.getAlphaX()(row, col - 1), - grid.getAlphaX()(row, col)) + - 2 * grid.getAlphaX()(row, col)) * - grid.getConcentrations()(row, col) + - calcAlphaIntercell(grid.getAlphaX()(row, col - 1), - grid.getAlphaX()(row, col)) * - grid.getConcentrations()(row, col - 1); -} - -// calculates horizontal change on one cell with a closed right boundary -template -static inline T calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, - int &col) { - - return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1), - grid.getAlphaX()(row, col)) * - (grid.getConcentrations()(row, col) - - grid.getConcentrations()(row, col - 1))); -} - -// checks boundary condition type for a cell on the right edge of grid -template -static inline T calcHorizontalChangeRightBoundary(Grid &grid, - Boundary &bc, int &row, - int &col) { - if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CONSTANT) { - return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col); - } else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CLOSED) { - return calcHorizontalChangeRightBoundaryClosed(grid, row, col); - } else { - throw_invalid_argument("Undefined Boundary Condition Type!"); - } -} - -// calculates vertical change on one cell with a constant top boundary -template -static inline T calcVerticalChangeTopBoundaryConstant(Grid &grid, - Boundary &bc, int &row, - int &col) { - - return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), - grid.getAlphaY()(row, col)) * - grid.getConcentrations()(row + 1, col) - - (calcAlphaIntercell(grid.getAlphaY()(row + 1, col), - grid.getAlphaY()(row, col)) + - 2 * grid.getAlphaY()(row, col)) * - grid.getConcentrations()(row, col) + - 2 * grid.getAlphaY()(row, col) * - bc.getBoundaryElementValue(BC_SIDE_TOP, col); -} - -// calculates vertical change on one cell with a closed top boundary -template -static inline T calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, - int &col) { - - return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), - grid.getAlphaY()(row, col)) * - (grid.getConcentrations()(row + 1, col) - - grid.getConcentrations()(row, col)); -} - -// checks boundary condition type for a cell on the top edge of grid -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) { - return calcVerticalChangeTopBoundaryClosed(grid, row, col); - } else { - throw_invalid_argument("Undefined Boundary Condition Type!"); - } -} - -// calculates vertical change on one cell with a constant bottom boundary -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) - - (calcAlphaIntercell(grid.getAlphaY()(row, col), - grid.getAlphaY()(row - 1, col)) + - 2 * grid.getAlphaY()(row, col)) * - grid.getConcentrations()(row, col) + - calcAlphaIntercell(grid.getAlphaY()(row, col), - grid.getAlphaY()(row - 1, col)) * - grid.getConcentrations()(row - 1, col); -} - -// calculates vertical change on one cell with a closed bottom boundary -template -static inline T calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, - int &col) { - - return -(calcAlphaIntercell(grid.getAlphaY()(row, col), - grid.getAlphaY()(row - 1, col)) * - (grid.getConcentrations()(row, col) - - grid.getConcentrations()(row - 1, col))); -} - -// checks boundary condition type for a cell on the bottom edge of grid -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) { - return calcVerticalChangeBottomBoundaryClosed(grid, row, col); - } else { - throw_invalid_argument("Undefined Boundary Condition Type!"); - } -} - -// FTCS solution for 1D grid -template -static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { - int colMax = grid.getCol(); - T deltaCol = grid.getDeltaCol(); - - // matrix for concentrations at time t+1 - RowMajMat concentrations_t1 = RowMajMat::Constant(1, colMax, 0); - - // only one row in 1D case -> row constant at index 0 - int row = 0; - - // inner cells - // independent of boundary condition type - for (int col = 1; col < colMax - 1; col++) { - concentrations_t1(row, col) = grid.getConcentrations()(row, col) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChange(grid, row, col)); - } - - // left boundary; hold column constant at index 0 - int col = 0; - concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChangeLeftBoundary(grid, bc, row, col)); - - // right boundary; hold column constant at max index - col = colMax - 1; - concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChangeRightBoundary(grid, bc, row, col)); - - // overwrite obsolete concentrations - grid.setConcentrations(concentrations_t1); -} - -// FTCS solution for 2D grid -template -static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, - int numThreads) { - int rowMax = grid.getRow(); - int colMax = grid.getCol(); - T deltaRow = grid.getDeltaRow(); - T deltaCol = grid.getDeltaCol(); - - // matrix for concentrations at time t+1 - RowMajMat concentrations_t1 = RowMajMat::Constant(rowMax, colMax, 0); - - // inner cells - // these are independent of the boundary condition type - // omp_set_num_threads(10); -#pragma omp parallel for num_threads(numThreads) - for (int row = 1; row < rowMax - 1; row++) { - for (int col = 1; col < colMax - 1; col++) { - concentrations_t1(row, col) = grid.getConcentrations()(row, col) + - timestep / (deltaRow * deltaRow) * - (calcVerticalChange(grid, row, col)) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChange(grid, row, col)); - } - } - - // boundary conditions - // left without corners / looping over rows - // hold column constant at index 0 - int col = 0; -#pragma omp parallel for num_threads(numThreads) - for (int row = 1; row < rowMax - 1; row++) { - concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + - timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col)); - } - - // right without corners / looping over rows - // hold column constant at max index - col = colMax - 1; -#pragma omp parallel for num_threads(numThreads) - for (int row = 1; row < rowMax - 1; row++) { - concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChangeRightBoundary(grid, bc, row, col)) + - timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col)); - } - - // top without corners / looping over columns - // hold row constant at index 0 - int row = 0; -#pragma omp parallel for num_threads(numThreads) - for (int col = 1; col < colMax - 1; col++) { - concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + - timestep / (deltaRow * deltaRow) * - (calcVerticalChangeTopBoundary(grid, bc, row, col)) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChange(grid, row, col)); - } - - // bottom without corners / looping over columns - // hold row constant at max index - row = rowMax - 1; -#pragma omp parallel for num_threads(numThreads) - for (int col = 1; col < colMax - 1; col++) { - concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + - timestep / (deltaRow * deltaRow) * - (calcVerticalChangeBottomBoundary(grid, bc, row, col)) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChange(grid, row, col)); - } - - // corner top left - // hold row and column constant at 0 - row = 0; - col = 0; - concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + - timestep / (deltaRow * deltaRow) * - (calcVerticalChangeTopBoundary(grid, bc, row, col)); - - // corner top right - // hold row constant at 0 and column constant at max index - row = 0; - col = colMax - 1; - concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChangeRightBoundary(grid, bc, row, col)) + - timestep / (deltaRow * deltaRow) * - (calcVerticalChangeTopBoundary(grid, bc, row, col)); - - // corner bottom left - // hold row constant at max index and column constant at 0 - row = rowMax - 1; - col = 0; - concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + - timestep / (deltaRow * deltaRow) * - (calcVerticalChangeBottomBoundary(grid, bc, row, col)); - - // corner bottom right - // hold row and column constant at max index - row = rowMax - 1; - col = colMax - 1; - concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + - timestep / (deltaCol * deltaCol) * - (calcHorizontalChangeRightBoundary(grid, bc, row, col)) + - timestep / (deltaRow * deltaRow) * - (calcVerticalChangeBottomBoundary(grid, bc, row, col)); - - // overwrite obsolete concentrations - grid.setConcentrations(concentrations_t1); - // } -} - -// entry point; differentiate between 1D and 2D grid -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) { - FTCS_2D(grid, bc, timestep, numThreads); - } else { - throw_invalid_argument( - "Error: Only 1- and 2-dimensional grids are defined!"); - } -} -} // namespace tug - -#endif // FTCS_H_ diff --git a/include/tug/Core/BTCS.hpp b/include/tug/Core/Numeric/BTCS.hpp similarity index 87% rename from include/tug/Core/BTCS.hpp rename to include/tug/Core/Numeric/BTCS.hpp index a1081af..5e561c5 100644 --- a/include/tug/Core/BTCS.hpp +++ b/include/tug/Core/Numeric/BTCS.hpp @@ -10,15 +10,16 @@ #ifndef BTCS_H_ #define BTCS_H_ -#include "Matrix.hpp" -#include "TugUtils.hpp" - #include #include -#include +#include +#include #include #include +#include +#include + #ifdef _OPENMP #include #else @@ -159,9 +160,9 @@ constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc, // creates a solution vector for next time step from the current state of // concentrations -template +template static Eigen::VectorX -createSolutionVector(const RowMajMat &concentrations, +createSolutionVector(const EigenType &concentrations, const RowMajMat &alphaX, const RowMajMat &alphaY, const std::vector> &bcLeft, const std::vector> &bcRight, @@ -351,25 +352,27 @@ static Eigen::VectorX ThomasAlgorithm(Eigen::SparseMatrix &A, // BTCS solution for 1D grid template -static void BTCS_1D(Grid &grid, Boundary &bc, T timestep, +static void BTCS_1D(SimulationInput &input, Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, Eigen::VectorX &b)) { - int length = grid.getLength(); - T sx = timestep / (grid.getDelta() * grid.getDelta()); + const std::size_t &length = input.colMax; + T sx = input.timestep / (input.deltaCol * input.deltaCol); Eigen::VectorX concentrations_t1(length); Eigen::SparseMatrix A; Eigen::VectorX b(length); - const auto &alpha = grid.getAlpha(); + const auto &alpha = input.alphaX; + + const auto &bc = input.boundaries; const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); const auto inner_bc = bc.getInnerBoundaryRow(0); - RowMajMat &concentrations = grid.getConcentrations(); + RowMajMatMap &concentrations = input.concentrations; int rowIndex = 0; A = createCoeffMatrix(alpha, bcLeft, bcRight, inner_bc, length, rowIndex, sx); // this is exactly same as in 2D @@ -396,29 +399,31 @@ static void BTCS_1D(Grid &grid, Boundary &bc, T timestep, // BTCS solution for 2D grid template -static void BTCS_2D(Grid &grid, Boundary &bc, T timestep, +static void BTCS_2D(SimulationInput &input, Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, Eigen::VectorX &b), int numThreads) { - int rowMax = grid.getRow(); - int colMax = grid.getCol(); - T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); - T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); + const std::size_t &rowMax = input.rowMax; + const std::size_t &colMax = input.colMax; + const T sx = input.timestep / (2 * input.deltaCol * input.deltaCol); + const T sy = input.timestep / (2 * input.deltaRow * input.deltaRow); RowMajMat concentrations_t1(rowMax, colMax); Eigen::SparseMatrix A; Eigen::VectorX b; - RowMajMat alphaX = grid.getAlphaX(); - RowMajMat alphaY = grid.getAlphaY(); + const RowMajMat &alphaX = input.alphaX; + const RowMajMat &alphaY = input.alphaY; + + const auto &bc = input.boundaries; 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); - RowMajMat &concentrations = grid.getConcentrations(); + RowMajMatMap &concentrations = input.concentrations; #pragma omp parallel for num_threads(numThreads) private(A, b) for (int i = 0; i < rowMax; i++) { @@ -432,44 +437,43 @@ static void BTCS_2D(Grid &grid, Boundary &bc, T timestep, } concentrations_t1.transposeInPlace(); - alphaX.transposeInPlace(); - alphaY.transposeInPlace(); + const RowMajMat alphaX_t = alphaX.transpose(); + const RowMajMat alphaY_t = alphaY.transpose(); #pragma omp parallel for num_threads(numThreads) private(A, b) for (int i = 0; i < colMax; i++) { auto inner_bc = bc.getInnerBoundaryCol(i); // swap alphas, boundary conditions and sx/sy for column-wise calculation - A = createCoeffMatrix(alphaY, bcTop, bcBottom, inner_bc, rowMax, i, sy); - b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, - bcLeft, bcRight, inner_bc, rowMax, i, sy, sx); + A = createCoeffMatrix(alphaY_t, bcTop, bcBottom, inner_bc, rowMax, i, sy); + b = createSolutionVector(concentrations_t1, alphaY_t, alphaX_t, bcTop, + bcBottom, bcLeft, bcRight, inner_bc, rowMax, i, sy, + sx); concentrations.col(i) = solverFunc(A, b); } } // entry point for EigenLU solver; differentiate between 1D and 2D grid -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) { - BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads); +template void BTCS_LU(SimulationInput &input, int numThreads) { + tug_assert(input.dim <= 2, + "Error: Only 1- and 2-dimensional grids are defined!"); + + if (input.dim == 1) { + BTCS_1D(input, EigenLUAlgorithm); } else { - throw_invalid_argument( - "Error: Only 1- and 2-dimensional grids are defined!"); + BTCS_2D(input.dim, EigenLUAlgorithm, numThreads); } } // entry point for Thomas algorithm solver; differentiate 1D and 2D grid -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) { - BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads); +template void BTCS_Thomas(SimulationInput &input, int numThreads) { + tug_assert(input.dim <= 2, + "Error: Only 1- and 2-dimensional grids are defined!"); + + if (input.dim == 1) { + BTCS_1D(input, ThomasAlgorithm); } else { - throw_invalid_argument( - "Error: Only 1- and 2-dimensional grids are defined!"); + BTCS_2D(input, ThomasAlgorithm, numThreads); } } } // namespace tug diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp new file mode 100644 index 0000000..ce1d5e8 --- /dev/null +++ b/include/tug/Core/Numeric/FTCS.hpp @@ -0,0 +1,240 @@ +/** + * @file FTCS.hpp + * @brief Implementation of heterogenous FTCS (forward time-centered space) + * solution of diffusion equation in 1D and 2D space. + * + */ + +#ifndef FTCS_H_ +#define FTCS_H_ + +#include "tug/Core/TugUtils.hpp" +#include +#include +#include +#include +#include + +#ifdef _OPENMP +#include +#else +#define omp_get_thread_num() 0 +#endif + +namespace tug { + +template +constexpr T calcChangeInner(T conc_c, T conc_left, T conc_right, T alpha_c, + T alpha_left, T alpha_right) { + const T alpha_center_left = calcAlphaIntercell(alpha_left, alpha_c); + const T alpha_center_right = calcAlphaIntercell(alpha_right, alpha_c); + + return alpha_center_left * conc_left - + (alpha_center_left + alpha_center_right) * conc_c + + alpha_center_right * conc_right; +} + +template +constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center, + T alpha_neighbor, const BoundaryElement &bc) { + const T alpha_center_neighbor = + calcAlphaIntercell(alpha_center, alpha_neighbor); + const T &conc_boundary = bc.getValue(); + + switch (bc.getType()) { + case BC_TYPE_CONSTANT: { + return 2 * alpha_center * conc_boundary - + (alpha_center_neighbor + 2 * alpha_center) * conc_c + + alpha_center_neighbor * conc_neighbor; + } + case BC_TYPE_CLOSED: { + return (alpha_center_neighbor * (conc_neighbor - conc_c)); + } + } + + tug_assert(false, "Undefined Boundary Condition Type!"); +} + +// FTCS solution for 1D grid +template static void FTCS_1D(SimulationInput &input) { + const std::size_t &colMax = input.colMax; + const T &deltaCol = input.deltaCol; + const T ×tep = input.timestep; + + RowMajMatMap &concentrations_grid = input.concentrations; + // matrix for concentrations at time t+1 + RowMajMat concentrations_t1 = concentrations_grid; + + const auto &alphaX = input.alphaX; + const auto &bc = input.boundaries; + + // only one row in 1D case -> row constant at index 0 + int row = 0; + + // inner cells + // independent of boundary condition type + for (int col = 1; col < colMax - 1; col++) { + const T &conc_c = concentrations_grid(row, col); + const T &conc_left = concentrations_grid(row, col - 1); + const T &conc_right = concentrations_grid(row, col + 1); + + const T &alpha_c = alphaX(row, col); + const T &alpha_left = alphaX(row, col - 1); + const T &alpha_right = alphaX(row, col + 1); + + concentrations_t1(row, col) = + concentrations_grid(row, col) + + timestep / (deltaCol * deltaCol) * + calcChangeInner(conc_c, conc_left, conc_right, alpha_c, alpha_left, + alpha_right); + } + + // left boundary; hold column constant at index 0 + { + int col = 0; + const T &conc_c = concentrations_grid(row, col); + const T &conc_right = concentrations_grid(row, col + 1); + const T &alpha_c = alphaX(row, col); + const T &alpha_right = alphaX(row, col + 1); + const BoundaryElement &bc_element = + input.boundaries.getBoundaryElement(BC_SIDE_LEFT, row); + + concentrations_t1(row, col) = + concentrations_grid(row, col) + + timestep / (deltaCol * deltaCol) * + calcChangeBoundary(conc_c, conc_right, alpha_c, alpha_right, + bc_element); + } + + // right boundary; hold column constant at max index + { + int col = colMax - 1; + const T &conc_c = concentrations_grid(row, col); + const T &conc_left = concentrations_grid(row, col - 1); + const T &alpha_c = alphaX(row, col); + const T &alpha_left = alphaX(row, col - 1); + const BoundaryElement &bc_element = + bc.getBoundaryElement(BC_SIDE_RIGHT, row); + + concentrations_t1(row, col) = + concentrations_grid(row, col) + + timestep / (deltaCol * deltaCol) * + calcChangeBoundary(conc_c, conc_left, alpha_c, alpha_left, + bc_element); + } + // overwrite obsolete concentrations + concentrations_grid = concentrations_t1; +} + +// FTCS solution for 2D grid +template +static void FTCS_2D(SimulationInput &input, int numThreads) { + const std::size_t &rowMax = input.rowMax; + const std::size_t &colMax = input.colMax; + const T &deltaRow = input.deltaRow; + const T &deltaCol = input.deltaCol; + const T ×tep = input.timestep; + + RowMajMatMap &concentrations_grid = input.concentrations; + + // matrix for concentrations at time t+1 + RowMajMat concentrations_t1 = concentrations_grid; + + const auto &alphaX = input.alphaX; + const auto &alphaY = input.alphaY; + const auto &bc = input.boundaries; + + const T sx = timestep / (deltaCol * deltaCol); + const T sy = timestep / (deltaRow * deltaRow); + +#pragma omp parallel for num_threads(numThreads) + for (std::size_t row_i = 0; row_i < rowMax; row_i++) { + for (std::size_t col_i = 0; col_i < colMax; col_i++) { + // horizontal change + T horizontal_change; + { + + const T &conc_c = concentrations_grid(row_i, col_i); + const T &alpha_c = alphaX(row_i, col_i); + + if (col_i == 0 || col_i == colMax - 1) { + // left or right boundary + const T &conc_neigbor = + concentrations_grid(row_i, col_i == 0 ? col_i + 1 : col_i - 1); + const T &alpha_neigbor = + alphaX(row_i, col_i == 0 ? col_i + 1 : col_i - 1); + + const BoundaryElement &bc_element = bc.getBoundaryElement( + col_i == 0 ? BC_SIDE_LEFT : BC_SIDE_RIGHT, row_i); + + horizontal_change = calcChangeBoundary(conc_c, conc_neigbor, alpha_c, + alpha_neigbor, bc_element); + } else { + // inner cell + const T &conc_left = concentrations_grid(row_i, col_i - 1); + const T &conc_right = concentrations_grid(row_i, col_i + 1); + + const T &alpha_left = alphaX(row_i, col_i - 1); + const T &alpha_right = alphaX(row_i, col_i + 1); + + horizontal_change = calcChangeInner(conc_c, conc_left, conc_right, + alpha_c, alpha_left, alpha_right); + } + } + + // vertical change + T vertical_change; + { + const T &conc_c = concentrations_grid(row_i, col_i); + const T &alpha_c = alphaY(row_i, col_i); + + if (row_i == 0 || row_i == rowMax - 1) { + // top or bottom boundary + const T &conc_neigbor = + concentrations_grid(row_i == 0 ? row_i + 1 : row_i - 1, col_i); + + const T &alpha_neigbor = + alphaY(row_i == 0 ? row_i + 1 : row_i - 1, col_i); + + const BoundaryElement &bc_element = bc.getBoundaryElement( + row_i == 0 ? BC_SIDE_TOP : BC_SIDE_BOTTOM, col_i); + + vertical_change = calcChangeBoundary(conc_c, conc_neigbor, alpha_c, + alpha_neigbor, bc_element); + } else { + // inner cell + const T &conc_bottom = concentrations_grid(row_i - 1, col_i); + const T &conc_top = concentrations_grid(row_i + 1, col_i); + + const T &alpha_bottom = alphaY(row_i - 1, col_i); + const T &alpha_top = alphaY(row_i + 1, col_i); + + vertical_change = calcChangeInner(conc_c, conc_bottom, conc_top, + alpha_c, alpha_bottom, alpha_top); + } + } + + concentrations_t1(row_i, col_i) = concentrations_grid(row_i, col_i) + + sx * horizontal_change + + sy * vertical_change; + } + } + + // overwrite obsolete concentrations + concentrations_grid = concentrations_t1; +} + +// entry point; differentiate between 1D and 2D grid +template void FTCS(SimulationInput &input, int &numThreads) { + tug_assert(input.dim <= 2, + "Error: Only 1- and 2-dimensional grids are defined!"); + + if (input.dim == 1) { + FTCS_1D(input); + } else { + FTCS_2D(input, numThreads); + } +} +} // namespace tug + +#endif // FTCS_H_ diff --git a/include/tug/Core/Numeric/SimulationInput.hpp b/include/tug/Core/Numeric/SimulationInput.hpp new file mode 100644 index 0000000..202619c --- /dev/null +++ b/include/tug/Core/Numeric/SimulationInput.hpp @@ -0,0 +1,21 @@ +#pragma once + +#include +#include + +namespace tug { + +template struct SimulationInput { + RowMajMatMap &concentrations; + const RowMajMat &alphaX; + const RowMajMat &alphaY; + const Boundary boundaries; + + const std::uint8_t dim; + const T timestep; + const std::size_t rowMax; + const std::size_t colMax; + const T deltaRow; + const T deltaCol; +}; +} // namespace tug \ No newline at end of file diff --git a/include/tug/Core/TugUtils.hpp b/include/tug/Core/TugUtils.hpp index 1fea096..2a0e9c2 100644 --- a/include/tug/Core/TugUtils.hpp +++ b/include/tug/Core/TugUtils.hpp @@ -1,9 +1,6 @@ -#ifndef TUGUTILS_H_ -#define TUGUTILS_H_ +#pragma once -#include -#include -#include +#include #define throw_invalid_argument(msg) \ throw std::invalid_argument(std::string(__FILE__) + ":" + \ @@ -24,6 +21,8 @@ duration.count(); \ }) +#define tug_assert(expr, msg) assert((expr) && msg) + // calculates arithmetic or harmonic mean of alpha between two cells template constexpr T calcAlphaIntercell(T alpha1, T alpha2, bool useHarmonic = true) { @@ -38,4 +37,3 @@ constexpr T calcAlphaIntercell(T alpha1, T alpha2, bool useHarmonic = true) { return 0.5 * (alpha1 + alpha2); } } -#endif // TUGUTILS_H_ diff --git a/include/tug/Simulation.hpp b/include/tug/Diffusion.hpp similarity index 54% rename from include/tug/Simulation.hpp rename to include/tug/Diffusion.hpp index 7f8632d..e353550 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Diffusion.hpp @@ -1,16 +1,14 @@ /** - * @file Simulation.hpp - * @brief API of Simulation class, that holds all information regarding a + * @file Diffusion.hpp + * @brief API of Diffusion class, that holds all information regarding a * specific simulation run like its timestep, number of iterations and output - * options. Simulation object also holds a predefined Grid and Boundary object. + * options. Diffusion object also holds a predefined Grid and Boundary object. * */ -#ifndef SIMULATION_H_ -#define SIMULATION_H_ +#pragma once -#include "Boundary.hpp" -#include "Grid.hpp" +#include "tug/Core/Matrix.hpp" #include #include #include @@ -19,9 +17,9 @@ #include #include -#include "Core/BTCS.hpp" -#include "Core/FTCS.hpp" -#include "Core/TugUtils.hpp" +#include +#include +#include #ifdef _OPENMP #include @@ -51,37 +49,6 @@ enum SOLVER { tridiagonal matrices */ }; -/** - * @brief Enum holding different options for .csv output. - * - */ -enum CSV_OUTPUT { - CSV_OUTPUT_OFF, /*!< do not produce csv output */ - CSV_OUTPUT_ON, /*!< produce csv output with last concentration matrix */ - CSV_OUTPUT_VERBOSE, /*!< produce csv output with all concentration matrices */ - CSV_OUTPUT_XTREME /*!< csv output like VERBOSE but additional boundary - conditions at beginning */ -}; - -/** - * @brief Enum holding different options for console output. - * - */ -enum CONSOLE_OUTPUT { - CONSOLE_OUTPUT_OFF, /*!< do not print any output to console */ - CONSOLE_OUTPUT_ON, /*!< print before and after concentrations to console */ - CONSOLE_OUTPUT_VERBOSE /*!< print all concentration matrices to console */ -}; - -/** - * @brief Enum holding different options for time measurement. - * - */ -enum TIME_MEASURE { - TIME_MEASURE_OFF, /*!< do not print any time measures */ - TIME_MEASURE_ON /*!< print time measure after last iteration */ -}; - /** * @brief The class forms the interface for performing the diffusion simulations * and contains all the methods for controlling the desired parameters, such as @@ -94,81 +61,94 @@ enum TIME_MEASURE { */ template -class Simulation { +class Diffusion : public BaseSimulationGrid { +private: + T timestep{-1}; + int innerIterations{1}; + int numThreads{omp_get_num_procs()}; + + RowMajMat alphaX; + RowMajMat alphaY; + + const std::vector approach_names = {"FTCS", "BTCS", "CRNI"}; + + static constexpr T DEFAULT_ALPHA = 1E-8; + + void init_alpha() { + this->alphaX = + RowMajMat::Constant(this->rows(), this->cols(), DEFAULT_ALPHA); + if (this->getDim() == 2) { + this->alphaY = + RowMajMat::Constant(this->rows(), this->cols(), DEFAULT_ALPHA); + } + } + public: /** - * @brief Set up a simulation environment. The timestep and number of - * iterations must be set. For the BTCS approach, the Thomas algorithm is used - * as the default linear equation solver as this is faster for tridiagonal - * matrices. CSV output, console output and time measure are off by - * default. Also, the number of cores is set to the maximum number of cores -1 - * by default. + * @brief Construct a new Diffusion object from a given Eigen matrix * - * @param grid Valid grid object - * @param bc Valid boundary condition object - * @param approach Approach to solving the problem. Either FTCS or BTCS. */ - Simulation(Grid &_grid, Boundary &_bc) : grid(_grid), bc(_bc){}; - - /** - * @brief Set the option to output the results to a CSV file. Off by default. - * - * - * @param csv_output Valid output option. The following options can be set - * here: - * - CSV_OUTPUT_OFF: do not produce csv output - * - CSV_OUTPUT_ON: produce csv output with last - * concentration matrix - * - CSV_OUTPUT_VERBOSE: produce csv output with all - * concentration matrices - * - CSV_OUTPUT_XTREME: produce csv output with all - * concentration matrices and simulation environment - */ - 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; + Diffusion(RowMajMat &origin) : BaseSimulationGrid(origin) { + init_alpha(); } /** - * @brief Set the options for outputting information to the console. Off by - * default. + * @brief Construct a new 2D Diffusion object from a given data pointer and + * the dimensions. * - * @param console_output Valid output option. The following options can be set - * here: - * - CONSOLE_OUTPUT_OFF: do not print any output to - * console - * - CONSOLE_OUTPUT_ON: print before and after - * concentrations to console - * - CONSOLE_OUTPUT_VERBOSE: print all concentration - * matrices to console */ - 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; + Diffusion(T *data, int rows, int cols) + : BaseSimulationGrid(data, rows, cols) { + init_alpha(); } /** - * @brief Set the Time Measure option. Off by default. + * @brief Construct a new 1D Diffusion object from a given data pointer and + * the length. * - * @param time_measure The following options are allowed: - * - TIME_MEASURE_OFF: Time of simulation is not printed - * to console - * - TIME_MEASURE_ON: Time of simulation run is printed to - * console */ - 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!"); - } + Diffusion(T *data, std::size_t length) : BaseSimulationGrid(data, length) { + init_alpha(); + } - this->time_measure = time_measure; + /** + * @brief Get the alphaX matrix. + * + * @return RowMajMat& Reference to the alphaX matrix. + */ + RowMajMat &getAlphaX() { return alphaX; } + + /** + * @brief Get the alphaY matrix. + * + * @return RowMajMat& Reference to the alphaY matrix. + */ + RowMajMat &getAlphaY() { + tug_assert( + this->getDim(), + "Grid is not two dimensional, there is no domain in y-direction!"); + + return alphaY; + } + + /** + * @brief Set the alphaX matrix. + * + * @param alphaX The new alphaX matrix. + */ + void setAlphaX(const RowMajMat &alphaX) { this->alphaX = alphaX; } + + /** + * @brief Set the alphaY matrix. + * + * @param alphaY The new alphaY matrix. + */ + void setAlphaY(const RowMajMat &alphaY) { + tug_assert( + this->getDim(), + "Grid is not two dimensional, there is no domain in y-direction!"); + + this->alphaY = alphaY; } /** @@ -177,33 +157,31 @@ public: * * @param timestep Valid timestep greater than zero. */ - void setTimestep(T timestep) { - if (timestep <= 0) { - throw_invalid_argument("Timestep has to be greater than zero."); - } + void setTimestep(T timestep) override { + tug_assert(timestep > 0, "Timestep has to be greater than zero."); if constexpr (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { T cfl; - if (grid.getDim() == 1) { + if (this->getDim() == 1) { - const T deltaSquare = grid.getDelta(); - const T maxAlpha = grid.getAlpha().maxCoeff(); + const T deltaSquare = this->deltaCol(); + const T maxAlpha = this->alphaX.maxCoeff(); // Courant-Friedrichs-Lewy condition cfl = deltaSquare / (4 * maxAlpha); - } else if (grid.getDim() == 2) { - const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + } else if (this->getDim() == 2) { + const T deltaColSquare = this->deltaCol() * this->deltaCol(); // will be 0 if 1D, else ... - const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + const T deltaRowSquare = this->deltaRow() * this->deltaRow(); const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); const T maxAlpha = - std::max(grid.getAlphaX().maxCoeff(), grid.getAlphaY().maxCoeff()); + std::max(this->alphaX.maxCoeff(), this->alphaY.maxCoeff()); cfl = minDeltaSquare / (4 * maxAlpha); } - const std::string dim = std::to_string(grid.getDim()) + "D"; + const std::string dim = std::to_string(this->getDim()) + "D"; const std::string &approachPrefix = this->approach_names[approach]; std::cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl @@ -244,20 +222,6 @@ public: */ T getTimestep() const { return this->timestep; } - /** - * @brief Set the desired iterations to be calculated. A value greater - * than zero must be specified here. Setting iterations is required. - * - * @param iterations Number of iterations to be simulated. - */ - 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 number of desired openMP Threads. * @@ -277,19 +241,12 @@ public: } } - /** - * @brief Return the currently set iterations to be calculated. - * - * @return int Number of iterations. - */ - int getIterations() const { return this->iterations; } - /** * @brief Outputs the current concentrations of the grid on the console. * */ - inline void printConcentrationsConsole() const { - std::cout << grid.getConcentrations() << std::endl; + void printConcentrationsConsole() const { + std::cout << this->getConcentrationMatrix() << std::endl; std::cout << std::endl; } @@ -309,9 +266,9 @@ public: // 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 row = std::to_string(this->rows()); + std::string col = std::to_string(this->cols()); + std::string numIterations = std::to_string(this->getIterations()); std::string filename = approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; @@ -330,7 +287,9 @@ public: // 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) { + if (this->getOutputCSV() == CSV_OUTPUT::XTREME) { + const auto &bc = this->getBoundaryConditions(); + Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "", " "); file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) @@ -365,7 +324,7 @@ public: } Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); - file << grid.getConcentrations().format(do_not_align) << std::endl; + file << this->getConcentrationMatrix().format(do_not_align) << std::endl; file << std::endl << std::endl; file.close(); } @@ -374,35 +333,42 @@ public: * @brief Method starts the simulation process with the previously set * parameters. */ - void run() { - if (this->timestep == -1) { - throw_invalid_argument("Timestep is not set!"); - } - if (this->iterations == -1) { - throw_invalid_argument("Number of iterations are not set!"); - } + void run() override { + tug_assert(this->getTimestep() > 0, "Timestep is not set!"); + tug_assert(this->getIterations() > 0, "Number of iterations are not set!"); std::string filename; - if (this->console_output > CONSOLE_OUTPUT_OFF) { + if (this->getOutputConsole() > CONSOLE_OUTPUT::OFF) { printConcentrationsConsole(); } - if (this->csv_output > CSV_OUTPUT_OFF) { + if (this->getOutputCSV() > CSV_OUTPUT::OFF) { filename = createCSVfile(); } auto begin = std::chrono::high_resolution_clock::now(); - if constexpr (approach == FTCS_APPROACH) { // FTCS case + SimulationInput sim_input = {.concentrations = + this->getConcentrationMatrix(), + .alphaX = this->getAlphaX(), + .alphaY = this->getAlphaY(), + .boundaries = this->getBoundaryConditions(), + .dim = this->getDim(), + .timestep = this->getTimestep(), + .rowMax = this->rows(), + .colMax = this->cols(), + .deltaRow = this->deltaRow(), + .deltaCol = this->deltaCol()}; - for (int i = 0; i < iterations * innerIterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + if constexpr (approach == FTCS_APPROACH) { // FTCS case + for (int i = 0; i < this->getIterations() * innerIterations; i++) { + if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) { printConcentrationsConsole(); } - if (csv_output >= CSV_OUTPUT_VERBOSE) { + if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) { printConcentrationsCSV(filename); } - FTCS(this->grid, this->bc, this->timestep, this->numThreads); + FTCS(sim_input, this->numThreads); // if (i % (iterations * innerIterations / 100) == 0) { // double percentage = (double)i / ((double)iterations * @@ -415,29 +381,28 @@ public: } else if constexpr (approach == BTCS_APPROACH) { // BTCS case if constexpr (solver == EIGEN_LU_SOLVER) { - for (int i = 0; i < iterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + for (int i = 0; i < this->getIterations(); i++) { + if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) { printConcentrationsConsole(); } - if (csv_output >= CSV_OUTPUT_VERBOSE) { + if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) { printConcentrationsCSV(filename); } - BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads); + BTCS_LU(sim_input, this->numThreads); } } else if constexpr (solver == THOMAS_ALGORITHM_SOLVER) { - for (int i = 0; i < iterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + for (int i = 0; i < this->getIterations(); i++) { + if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) { printConcentrationsConsole(); } - if (csv_output >= CSV_OUTPUT_VERBOSE) { + if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) { printConcentrationsCSV(filename); } - BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); + BTCS_Thomas(sim_input, this->numThreads); } } - } else if constexpr (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case @@ -449,22 +414,22 @@ public: RowMajMat concentrations; RowMajMat concentrationsFTCS; RowMajMat concentrationsResult; - for (int i = 0; i < iterations * innerIterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + for (int i = 0; i < this->getIterations() * innerIterations; i++) { + if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) { printConcentrationsConsole(); } - if (csv_output >= CSV_OUTPUT_VERBOSE) { + if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) { printConcentrationsCSV(filename); } - concentrations = grid.getConcentrations(); + concentrations = this->getConcentrationMatrix(); FTCS(this->grid, this->bc, this->timestep, this->numThreads); - concentrationsFTCS = grid.getConcentrations(); - grid.setConcentrations(concentrations); - BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); - concentrationsResult = - beta * concentrationsFTCS + (1 - beta) * grid.getConcentrations(); - grid.setConcentrations(concentrationsResult); + concentrationsFTCS = this->getConcentrationMatrix(); + this->getConcentrationMatrix() = concentrations; + BTCS_Thomas(sim_input, this->numThreads); + concentrationsResult = beta * concentrationsFTCS + + (1 - beta) * this->getConcentrationMatrix(); + this->getConcentrationMatrix() = concentrationsResult; } } @@ -472,33 +437,18 @@ public: auto milliseconds = std::chrono::duration_cast(end - begin); - if (this->console_output > CONSOLE_OUTPUT_OFF) { + if (this->getOutputConsole() > CONSOLE_OUTPUT::OFF) { printConcentrationsConsole(); } - if (this->csv_output > CSV_OUTPUT_OFF) { + if (this->getOutputCSV() > CSV_OUTPUT::OFF) { printConcentrationsCSV(filename); } - if (this->time_measure > TIME_MEASURE_OFF) { + if (this->getTimeMeasure() > TIME_MEASURE::OFF) { const std::string &approachString = this->approach_names[approach]; - const std::string dimString = std::to_string(grid.getDim()) + "D"; + const std::string dimString = std::to_string(this->getDim()) + "D"; std::cout << approachString << dimString << ":: run() finished in " << milliseconds.count() << "ms" << std::endl; } } - -private: - T timestep{-1}; - int iterations{-1}; - int innerIterations{1}; - int numThreads{omp_get_num_procs()}; - CSV_OUTPUT csv_output{CSV_OUTPUT_OFF}; - CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT_OFF}; - TIME_MEASURE time_measure{TIME_MEASURE_OFF}; - - Grid &grid; - Boundary &bc; - - const std::vector approach_names = {"FTCS", "BTCS", "CRNI"}; }; } // namespace tug -#endif // SIMULATION_H_ diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp deleted file mode 100644 index b8e471a..0000000 --- a/include/tug/Grid.hpp +++ /dev/null @@ -1,392 +0,0 @@ -#ifndef GRID_H_ -#define GRID_H_ - -/** - * @file Grid.hpp - * @brief API of Grid class, that holds a matrix with concenctrations and a - * respective matrix/matrices of alpha coefficients. - * - */ - -#include "Core/Matrix.hpp" -#include -#include -#include -#include -#include - -namespace tug { - -/** - * @brief Holds a matrix with concenctration and respective matrix/matrices of - * alpha coefficients. - * - * @tparam T Type to be used for matrices, e.g. double or float - */ -template class Grid { -public: - /** - * @brief Constructs a new 1D-Grid object of a given length, which holds a - * matrix with concentrations and a respective matrix of alpha coefficients. - * The domain length is per default the same as the length. The - * concentrations are all 20 by default and the alpha coefficients are 1. - * - * @param length Length of the 1D-Grid. Must be greater than 3. - */ - 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 = - static_cast(this->domainCol) / static_cast(this->col); // -> 1 - - this->concentrations = RowMajMat::Constant(1, col, MAT_INIT_VAL); - this->alphaX = RowMajMat::Constant(1, col, MAT_INIT_VAL); - } - - /** - * @brief Constructs a new 2D-Grid object of given dimensions, which holds a - * matrix with concentrations and the respective matrices of alpha coefficient - * for each direction. The domain in x- and y-direction is per default equal - * to the col length and row length, respectively. The concentrations are all - * 20 by default across the entire grid and the alpha coefficients 1 in both - * directions. - * - * @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) - : row(_row), col(_col), domainRow(_row), domainCol(_col) { - if (row <= 1 || col <= 1) { - throw std::invalid_argument( - "At least one dimension is 1. Use 1D grid for better results."); - } - - this->dim = 2; - this->deltaRow = - static_cast(this->domainRow) / static_cast(this->row); // -> 1 - this->deltaCol = - static_cast(this->domainCol) / static_cast(this->col); // -> 1 - - this->concentrations = RowMajMat::Constant(row, col, MAT_INIT_VAL); - this->alphaX = RowMajMat::Constant(row, col, MAT_INIT_VAL); - this->alphaY = RowMajMat::Constant(row, col, MAT_INIT_VAL); - } - - /** - * @brief Sets the concentrations matrix for a 1D or 2D-Grid. - * - * @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 RowMajMat &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 Sets the concentrations matrix for a 1D or 2D-Grid. - * - * @param concentrations A pointer to an array holding the concentrations. - * Array must have correct dimensions as defined in row and col. (Or length, - * in 1D case). There is no check for correct dimensions, so be careful! - */ - void setConcentrations(T *concentrations) { - tug::RowMajMatMap map(concentrations, this->row, this->col); - this->concentrations = map; - } - - /** - * @brief Gets the concentrations matrix for a Grid. - * - * @return An Eigen3 matrix holding the concentrations and having - * the same dimensions as the grid. - */ - auto &getConcentrations() { return this->concentrations; } - - /** - * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one - * dimensional. - * - * @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 RowMajMat &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 1D-Grid. Grid must be one - * dimensional. - * - * @param alpha A pointer to an array holding the alpha coefficients. Array - * must have correct dimensions as defined in length. There is no check for - * correct dimensions, so be careful! - */ - void setAlpha(T *alpha) { - if (dim != 1) { - throw std::invalid_argument( - "Grid is not one dimensional, you should probably " - "use 2D setter function!"); - } - RowMajMatMap map(alpha, 1, this->col); - this->alphaX = map; - } - - /** - * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two - * dimensional. - * - * @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 MatrixX holding the alpha coefficients in - * y-direction. Matrix must be of same size as the grid. - */ - void setAlpha(const RowMajMat &alphaX, const RowMajMat &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 Set the alpha coefficients of a 2D-Grid. Grid must be two - * dimensional. - * - * @param alphaX A pointer to an array holding the alpha coefficients in - * x-direction. Array must have correct dimensions as defined in row and col. - * There is no check for correct dimensions, so be careful! - * @param alphaY A pointer to an array holding the alpha coefficients in - * y-direction. Array must have correct dimensions as defined in row and col. - * There is no check for correct dimensions, so be careful! - */ - void setAlpha(T *alphaX, T *alphaY) { - if (dim != 2) { - throw std::invalid_argument( - "Grid is not two dimensional, you should probably " - "use 1D setter function!"); - } - RowMajMatMap mapX(alphaX, this->row, this->col); - RowMajMatMap mapY(alphaY, this->row, this->col); - this->alphaX = mapX; - this->alphaY = mapY; - } - - /** - * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one - * dimensional. - * - * @return A matrix with 1 row holding the alpha coefficients. - */ - const auto &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 A matrix holding the alpha coefficients in x-direction. - */ - const auto &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 A matrix holding the alpha coefficients in y-direction. - */ - const auto &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 Dimensions, either 1 or 2. - */ - int getDim() const { return this->dim; } - - /** - * @brief Gets length of 1D grid. Must be one dimensional grid. - * - * @return Length of 1D grid. - */ - 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 Number of rows. - */ - int getRow() const { return this->row; } - - /** - * @brief Gets the number of columns of the grid. - * - * @return Number of columns. - */ - 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) { - 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. - */ - 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 Delta value. - */ - 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 Delta value in x-direction. - */ - T getDeltaCol() const { return this->deltaCol; } - - /** - * @brief Gets the delta value in y-direction. Must be two dimensional grid. - * - * @return Delta value in y-direction. - */ - 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 - 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) - - RowMajMat concentrations; // Matrix holding grid concentrations - RowMajMat alphaX; // Matrix holding alpha coefficients in x-direction - RowMajMat alphaY; // Matrix holding alpha coefficients in y-direction - - static constexpr T MAT_INIT_VAL = 0; -}; - -using Grid64 = Grid; -using Grid32 = Grid; -} // namespace tug -#endif // GRID_H_ diff --git a/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp index 9ed7ccb..15f199f 100644 --- a/naaice/BTCS_2D_NAAICE.cpp +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -5,8 +5,7 @@ #include #include #include -#include -#include +#include #include #include "files.hpp" @@ -105,15 +104,14 @@ int main(int argc, char *argv[]) { // create a grid with a 5 x 10 field constexpr int row = 5; constexpr int col = 10; - 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); Eigen::MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col); - grid.setConcentrations(concentrations); + Grid64 grid(concentrations); + grid.setDomain(0.005, 0.01); const double sum_init = concentrations.sum(); // // (optional) set alphas of the grid, e.g.: @@ -141,7 +139,7 @@ int main(int argc, char *argv[]) { // // ************************ // set up a simulation environment - Simulation simulation = Simulation(grid, bc); // grid,boundary + Diffusion simulation(grid, bc); // grid,boundary // set the timestep of the simulation simulation.setTimestep(360); // timestep diff --git a/naaice/NAAICE_dble_vs_float.cpp b/naaice/NAAICE_dble_vs_float.cpp index ac3718c..2238375 100644 --- a/naaice/NAAICE_dble_vs_float.cpp +++ b/naaice/NAAICE_dble_vs_float.cpp @@ -7,11 +7,11 @@ #include #include #include -#include #include #include -#include "files.hpp" +#include +#include using namespace tug; @@ -114,15 +114,14 @@ template int doWork(int ngrid) { std::cout << name << " grid: " << ngrid << std::endl; - Grid grid(ngrid, ngrid); // Grid64 grid(ngrid, ngrid); // (optional) set the domain, e.g.: - grid.setDomain(0.1, 0.1); Eigen::MatrixX initConc64 = Eigen::MatrixX::Constant(ngrid, ngrid, 0); initConc64(50, 50) = 1E-6; - grid.setConcentrations(initConc64); + Grid grid(initConc64); + grid.setDomain(0.1, 0.1); const T sum_init64 = initConc64.sum(); @@ -142,8 +141,7 @@ template int doWork(int ngrid) { bc.setBoundarySideClosed(BC_SIDE_BOTTOM); // set up a simulation environment - Simulation Sim = - Simulation(grid, bc); // grid_64,boundary,simulation-approach + Diffusion Sim(grid, bc); // grid_64,boundary,simulation-approach // Sim64.setSolver(THOMAS_ALGORITHM_SOLVER); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b509370..2a8d780 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -10,12 +10,12 @@ FetchContent_MakeAvailable(googletest) add_executable(testTug -testSimulation.cpp -testGrid.cpp +setup.cpp +testDiffusion.cpp testFTCS.cpp testBoundary.cpp ) -target_link_libraries(testTug tug GTest::gtest_main) +target_link_libraries(testTug tug GTest::gtest) include(GoogleTest) gtest_discover_tests(testTug) @@ -25,4 +25,4 @@ get_filename_component(testSimulationCSV "FTCS_11_11_7000.csv" REALPATH) # set relative path in header file configure_file(testSimulation.hpp.in testSimulation.hpp) # include test directory with generated header file from above -target_include_directories(testTug PUBLIC "${CMAKE_CURRENT_BINARY_DIR}" "${PROJECT_SOURCE_DIR}/src") \ No newline at end of file +target_include_directories(testTug PUBLIC "${CMAKE_CURRENT_BINARY_DIR}" "${PROJECT_SOURCE_DIR}/src") diff --git a/test/TestUtils.hpp b/test/TestUtils.hpp index d9eeef9..b86b6d0 100644 --- a/test/TestUtils.hpp +++ b/test/TestUtils.hpp @@ -1,3 +1,4 @@ +#include "tug/Core/Matrix.hpp" #include #include #include @@ -8,7 +9,7 @@ #define TUG_TEST(x) TEST(Tug, x) -inline Eigen::MatrixXd CSV2Eigen(std::string file2Convert) { +inline tug::RowMajMat CSV2Eigen(std::string file2Convert) { std::vector matrixEntries; @@ -31,21 +32,20 @@ inline Eigen::MatrixXd CSV2Eigen(std::string file2Convert) { } } - return Eigen::Map< - Eigen::Matrix>( - matrixEntries.data(), matrixRowNumber, - matrixEntries.size() / matrixRowNumber); + return tug::RowMajMatMap(matrixEntries.data(), matrixRowNumber, + matrixEntries.size() / matrixRowNumber); } -inline bool checkSimilarity(Eigen::MatrixXd a, Eigen::MatrixXd b, +inline bool checkSimilarity(tug::RowMajMat &a, + tug::RowMajMatMap &b, double precision = 1e-5) { return a.isApprox(b, precision); } -inline bool checkSimilarityV2(Eigen::MatrixXd a, Eigen::MatrixXd b, - double maxDiff) { +inline bool checkSimilarityV2(tug::RowMajMat &a, + tug::RowMajMatMap &b, double maxDiff) { - Eigen::MatrixXd diff = a - b; + tug::RowMajMat diff = a - b; double maxCoeff = diff.maxCoeff(); return abs(maxCoeff) < maxDiff; } diff --git a/test/setup.cpp b/test/setup.cpp index 0a3f254..b885e52 100644 --- a/test/setup.cpp +++ b/test/setup.cpp @@ -1,2 +1,7 @@ -#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN -#include +#include + +int main(int argc, char **argv) { + ::testing::InitGoogleTest(&argc, argv); + GTEST_FLAG_SET(death_test_style, "threadsafe"); + return RUN_ALL_TESTS(); +} \ No newline at end of file diff --git a/test/testBoundary.cpp b/test/testBoundary.cpp index 0b2601b..38cc14a 100644 --- a/test/testBoundary.cpp +++ b/test/testBoundary.cpp @@ -1,3 +1,5 @@ +#include "gtest/gtest.h" +#include #include #include #include @@ -15,7 +17,7 @@ BOUNDARY_TEST(Element) { EXPECT_NO_THROW(BoundaryElement()); EXPECT_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED); EXPECT_DOUBLE_EQ(boundaryElementClosed.getValue(), -1); - EXPECT_ANY_THROW(boundaryElementClosed.setValue(0.2)); + EXPECT_THROW(boundaryElementClosed.setValue(0.2), std::invalid_argument); BoundaryElement boundaryElementConstant = BoundaryElement(0.1); EXPECT_NO_THROW(BoundaryElement(0.1)); @@ -26,10 +28,8 @@ BOUNDARY_TEST(Element) { } BOUNDARY_TEST(Class) { - Grid grid1D = Grid64(10); - Grid grid2D = Grid64(10, 12); - Boundary boundary1D = Boundary(grid1D); - Boundary boundary2D = Boundary(grid2D); + Boundary boundary1D(10); + Boundary boundary2D(10, 12); vector> boundary1DVector(1, BoundaryElement(1.0)); constexpr double inner_condition_value = -5; @@ -43,31 +43,34 @@ BOUNDARY_TEST(Class) { col_ibc[0] = innerBoundary; { - EXPECT_NO_THROW(Boundary boundary(grid1D)); EXPECT_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1); EXPECT_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1); EXPECT_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CLOSED); - EXPECT_ANY_THROW(boundary1D.getBoundarySide(BC_SIDE_TOP)); - EXPECT_ANY_THROW(boundary1D.getBoundarySide(BC_SIDE_BOTTOM)); + EXPECT_DEATH(boundary1D.getBoundarySide(BC_SIDE_TOP), + ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*"); + EXPECT_DEATH(boundary1D.getBoundarySide(BC_SIDE_BOTTOM), + ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*"); EXPECT_NO_THROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT)); - EXPECT_ANY_THROW(boundary1D.setBoundarySideClosed(BC_SIDE_TOP)); + EXPECT_DEATH(boundary1D.setBoundarySideClosed(BC_SIDE_TOP), + ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*"); EXPECT_NO_THROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); EXPECT_DOUBLE_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); - EXPECT_ANY_THROW(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2)); + EXPECT_DEATH(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2), + ".*Index is selected either too large or too small.*"); EXPECT_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); EXPECT_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); EXPECT_NO_THROW(boundary1D.setInnerBoundary(0, inner_condition_value)); - EXPECT_ANY_THROW(boundary1D.setInnerBoundary(0, 0, inner_condition_value)); + EXPECT_DEATH(boundary1D.setInnerBoundary(0, 0, inner_condition_value), + ".*only available for 2D grids.*"); EXPECT_EQ(boundary1D.getInnerBoundary(0), innerBoundary); EXPECT_FALSE(boundary1D.getInnerBoundary(1).first); } { - EXPECT_NO_THROW(Boundary boundary(grid1D)); EXPECT_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10); EXPECT_EQ(boundary2D.getBoundarySide(BC_SIDE_RIGHT).size(), 10); EXPECT_EQ(boundary2D.getBoundarySide(BC_SIDE_TOP).size(), 12); @@ -80,13 +83,15 @@ BOUNDARY_TEST(Class) { EXPECT_NO_THROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP)); EXPECT_NO_THROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); EXPECT_DOUBLE_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); - EXPECT_ANY_THROW(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12)); + EXPECT_DEATH(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12), + ".*too large or too small.*"); EXPECT_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); EXPECT_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); - EXPECT_ANY_THROW(boundary2D.setInnerBoundary(0, inner_condition_value)); + EXPECT_DEATH(boundary2D.setInnerBoundary(0, inner_condition_value), + ".* 1D .*"); EXPECT_NO_THROW(boundary2D.setInnerBoundary(0, 1, inner_condition_value)); EXPECT_EQ(boundary2D.getInnerBoundary(0, 1), innerBoundary); EXPECT_FALSE(boundary2D.getInnerBoundary(0, 2).first); diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp new file mode 100644 index 0000000..15fa18b --- /dev/null +++ b/test/testDiffusion.cpp @@ -0,0 +1,216 @@ +#include "TestUtils.hpp" +#include "tug/Core/Matrix.hpp" +#include "gtest/gtest.h" +#include +#include +#include + +#include +#include + +// include the configured header file +#include + +#define DIFFUSION_TEST(x) TEST(Diffusion, x) + +using namespace Eigen; +using namespace std; +using namespace tug; + +constexpr int row = 11; +constexpr int col = 11; + +template +Diffusion setupSimulation(RowMajMat &concentrations, + double timestep, int iterations) { + int domain_row = 10; + int domain_col = 10; + + // Grid + // RowMajMat concentrations = MatrixXd::Constant(row, col, 0); + concentrations(5, 5) = 1; + + Diffusion diffusiongrid(concentrations); + + diffusiongrid.getConcentrationMatrix() = concentrations; + diffusiongrid.setDomain(domain_row, domain_col); + + diffusiongrid.setTimestep(timestep); + diffusiongrid.setIterations(iterations); + diffusiongrid.setDomain(domain_row, domain_col); + + MatrixXd alpha = MatrixXd::Constant(row, col, 1); + for (int i = 0; i < 5; i++) { + for (int j = 0; j < 6; j++) { + alpha(i, j) = 0.01; + } + } + for (int i = 0; i < 5; i++) { + for (int j = 6; j < 11; j++) { + alpha(i, j) = 0.001; + } + } + for (int i = 5; i < 11; i++) { + for (int j = 6; j < 11; j++) { + alpha(i, j) = 0.1; + } + } + diffusiongrid.setAlphaX(alpha); + diffusiongrid.setAlphaY(alpha); + + return diffusiongrid; +} + +constexpr double timestep = 0.001; +constexpr double iterations = 7000; + +DIFFUSION_TEST(EqualityFTCS) { + // set string from the header file + string test_path = testSimulationCSVDir; + RowMajMat reference = CSV2Eigen(test_path); + cout << "FTCS Test: " << endl; + + RowMajMat concentrations = MatrixXd::Constant(row, col, 0); + + Diffusion sim = + setupSimulation(concentrations, timestep, iterations); + + // Boundary bc = Boundary(grid); + + // Simulation + + // Diffusion sim(grid, bc); + // sim.setOutputConsole(CONSOLE_OUTPUT_ON); + // sim.setTimestep(timestep); + // sim.setIterations(iterations); + sim.run(); + + cout << endl; + EXPECT_TRUE(checkSimilarity(reference, sim.getConcentrationMatrix(), 0.1)); +} + +DIFFUSION_TEST(EqualityBTCS) { + // set string from the header file + string test_path = testSimulationCSVDir; + RowMajMat reference = CSV2Eigen(test_path); + cout << "BTCS Test: " << endl; + + RowMajMat concentrations = MatrixXd::Constant(row, col, 0); + + Diffusion sim = + setupSimulation(concentrations, timestep, + iterations); // Boundary + + // Boundary bc = Boundary(grid); + + // Simulation + // Diffusion sim(grid, bc); + // sim.setOutputConsole(CONSOLE_OUTPUT_ON); + // sim.setTimestep(timestep); + // sim.setIterations(iterations); + sim.run(); + + cout << endl; + EXPECT_TRUE(checkSimilarityV2(reference, sim.getConcentrationMatrix(), 0.01)); +} + +DIFFUSION_TEST(InitializeEnvironment) { + int rc = 12; + RowMajMat concentrations(rc, rc); + // Grid64 grid(concentrations); + // Boundary boundary(grid); + + EXPECT_NO_FATAL_FAILURE(Diffusion sim(concentrations)); +} + +// DIFFUSION_TEST(SimulationEnvironment) { +// int rc = 12; +// Eigen::MatrixXd concentrations(rc, rc); +// Grid64 grid(concentrations); +// grid.initAlpha(); +// Boundary boundary(grid); +// Diffusion sim(grid, boundary); + +// EXPECT_EQ(sim.getIterations(), 1); + +// EXPECT_NO_THROW(sim.setIterations(2000)); +// EXPECT_EQ(sim.getIterations(), 2000); +// EXPECT_THROW(sim.setIterations(-300), std::invalid_argument); + +// EXPECT_NO_THROW(sim.setTimestep(0.1)); +// EXPECT_DOUBLE_EQ(sim.getTimestep(), 0.1); +// EXPECT_DEATH(sim.setTimestep(-0.3), ".* greater than zero.*"); +// } + +DIFFUSION_TEST(ClosedBoundaries) { + constexpr std::uint32_t nrows = 5; + constexpr std::uint32_t ncols = 5; + + RowMajMat concentrations = + RowMajMat::Constant(nrows, ncols, 1.0); + RowMajMat alphax = RowMajMat::Constant(nrows, ncols, 1E-5); + RowMajMat alphay = RowMajMat::Constant(nrows, ncols, 1E-5); + + Diffusion sim(concentrations); + sim.getAlphaX() = alphax; + sim.getAlphaY() = alphay; + + // tug::Grid64 grid(concentrations); + + // grid.setAlpha(alphax, alphay); + + // tug::Boundary bc(grid); + auto &bc = sim.getBoundaryConditions(); + bc.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1.0); + bc.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1.0); + bc.setBoundarySideConstant(tug::BC_SIDE_TOP, 1.0); + bc.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1.0); + + // tug::Diffusion sim(grid, bc); + sim.setTimestep(1); + sim.setIterations(1); + + RowMajMat input_values(concentrations); + sim.run(); + + EXPECT_TRUE( + checkSimilarityV2(input_values, sim.getConcentrationMatrix(), 1E-12)); +} + +DIFFUSION_TEST(ConstantInnerCell) { + constexpr std::uint32_t nrows = 5; + constexpr std::uint32_t ncols = 5; + + RowMajMat concentrations = + RowMajMat::Constant(nrows, ncols, 1.0); + RowMajMat alphax = RowMajMat::Constant(nrows, ncols, 1E-5); + RowMajMat alphay = RowMajMat::Constant(nrows, ncols, 1E-5); + + Diffusion sim(concentrations); + sim.getAlphaX() = alphax; + sim.getAlphaY() = alphay; + + // tug::Grid64 grid(concentrations); + // grid.setAlpha(alphax, alphay); + + // tug::Boundary bc(grid); + auto &bc = sim.getBoundaryConditions(); + // inner + bc.setInnerBoundary(2, 2, 0); + + // tug::Diffusion sim(grid, bc); + sim.setTimestep(1); + sim.setIterations(1); + + MatrixXd input_values(concentrations); + sim.run(); + + const auto &concentrations_result = sim.getConcentrationMatrix(); + + EXPECT_DOUBLE_EQ(concentrations_result(2, 2), 0); + EXPECT_LT(concentrations_result.sum(), input_values.sum()); + + EXPECT_FALSE((concentrations_result.array() > 1.0).any()); + + EXPECT_FALSE((concentrations_result.array() < 0.0).any()); +} diff --git a/test/testGrid.cpp b/test/testGrid.cpp deleted file mode 100644 index 4539389..0000000 --- a/test/testGrid.cpp +++ /dev/null @@ -1,262 +0,0 @@ -#include -#include -#include - -#include - -using namespace Eigen; -using namespace std; -using namespace tug; - -#define GRID_TEST(x) TEST(Grid, x) - -GRID_TEST(InvalidConstructorParams) { - EXPECT_ANY_THROW(Grid64(2)); - EXPECT_ANY_THROW(Grid64(1, 4)); - EXPECT_ANY_THROW(Grid64(4, 1)); -} - -GRID_TEST(Grid64OneDimensional) { - int l = 12; - Grid64 grid(l); - - { - EXPECT_EQ(grid.getDim(), 1); - EXPECT_EQ(grid.getLength(), l); - EXPECT_EQ(grid.getCol(), l); - EXPECT_EQ(grid.getRow(), 1); - - EXPECT_EQ(grid.getConcentrations().rows(), 1); - EXPECT_EQ(grid.getConcentrations().cols(), l); - EXPECT_EQ(grid.getAlpha().rows(), 1); - EXPECT_EQ(grid.getAlpha().cols(), l); - EXPECT_EQ(grid.getDeltaCol(), 1); - - EXPECT_ANY_THROW(grid.getAlphaX()); - EXPECT_ANY_THROW(grid.getAlphaY()); - EXPECT_ANY_THROW(grid.getDeltaRow()); - } - - { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(1, l, 3); - EXPECT_NO_THROW(grid.setConcentrations(concentrations)); - - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - } - - { - // correct alpha matrix - MatrixXd alpha = MatrixXd::Constant(1, l, 3); - EXPECT_NO_THROW(grid.setAlpha(alpha)); - - EXPECT_ANY_THROW(grid.setAlpha(alpha, alpha)); - - grid.setAlpha(alpha); - EXPECT_EQ(grid.getAlpha(), alpha); - EXPECT_ANY_THROW(grid.getAlphaX()); - EXPECT_ANY_THROW(grid.getAlphaY()); - - // false alpha matrix - MatrixXd wAlpha = MatrixXd::Constant(3, l, 2); - EXPECT_ANY_THROW(grid.setAlpha(wAlpha)); - } - - { - int d = 8; - // set 1D domain - EXPECT_NO_THROW(grid.setDomain(d)); - - // set 2D domain - EXPECT_ANY_THROW(grid.setDomain(d, d)); - - grid.setDomain(d); - EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(d) / double(l)); - EXPECT_ANY_THROW(grid.getDeltaRow()); - - // set too small domain - d = 0; - EXPECT_ANY_THROW(grid.setDomain(d)); - d = -2; - EXPECT_ANY_THROW(grid.setDomain(d)); - } -} - -GRID_TEST(Grid64Quadratic) { - int rc = 12; - Grid64 grid(rc, rc); - - { - EXPECT_EQ(grid.getDim(), 2); - EXPECT_ANY_THROW(grid.getLength()); - EXPECT_EQ(grid.getCol(), rc); - EXPECT_EQ(grid.getRow(), rc); - - EXPECT_EQ(grid.getConcentrations().rows(), rc); - EXPECT_EQ(grid.getConcentrations().cols(), rc); - EXPECT_ANY_THROW(grid.getAlpha()); - - EXPECT_EQ(grid.getAlphaX().rows(), rc); - EXPECT_EQ(grid.getAlphaX().cols(), rc); - EXPECT_EQ(grid.getAlphaY().rows(), rc); - EXPECT_EQ(grid.getAlphaY().cols(), rc); - EXPECT_EQ(grid.getDeltaRow(), 1); - EXPECT_EQ(grid.getDeltaCol(), 1); - } - - { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2); - EXPECT_NO_THROW(grid.setConcentrations(concentrations)); - - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(rc, rc + 3, 1); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(rc + 3, rc, 4); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(rc + 2, rc + 4, 6); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - } - - { - // correct alpha matrices - MatrixXd alphax = MatrixXd::Constant(rc, rc, 2); - MatrixXd alphay = MatrixXd::Constant(rc, rc, 4); - EXPECT_NO_THROW(grid.setAlpha(alphax, alphay)); - - EXPECT_ANY_THROW(grid.setAlpha(alphax)); - - grid.setAlpha(alphax, alphay); - EXPECT_EQ(grid.getAlphaX(), alphax); - EXPECT_EQ(grid.getAlphaY(), alphay); - EXPECT_ANY_THROW(grid.getAlpha()); - - // false alpha matrices - alphax = MatrixXd::Constant(rc + 3, rc + 1, 3); - EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay)); - alphay = MatrixXd::Constant(rc + 2, rc + 1, 3); - EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay)); - } - - { - int dr = 8; - int dc = 9; - - // set 1D domain - EXPECT_ANY_THROW(grid.setDomain(dr)); - - // set 2D domain - EXPECT_NO_THROW(grid.setDomain(dr, dc)); - - grid.setDomain(dr, dc); - EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(dc) / double(rc)); - EXPECT_DOUBLE_EQ(grid.getDeltaRow(), double(dr) / double(rc)); - - // set too small domain - dr = 0; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); - dr = 8; - dc = 0; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); - dr = -2; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); - } -} - -GRID_TEST(Grid64NonQuadratic) { - int r = 12; - int c = 15; - Grid64 grid(r, c); - - { - EXPECT_EQ(grid.getDim(), 2); - EXPECT_ANY_THROW(grid.getLength()); - EXPECT_EQ(grid.getCol(), c); - EXPECT_EQ(grid.getRow(), r); - - EXPECT_EQ(grid.getConcentrations().rows(), r); - EXPECT_EQ(grid.getConcentrations().cols(), c); - EXPECT_ANY_THROW(grid.getAlpha()); - - EXPECT_EQ(grid.getAlphaX().rows(), r); - EXPECT_EQ(grid.getAlphaX().cols(), c); - EXPECT_EQ(grid.getAlphaY().rows(), r); - EXPECT_EQ(grid.getAlphaY().cols(), c); - EXPECT_EQ(grid.getDeltaRow(), 1); - EXPECT_EQ(grid.getDeltaCol(), 1); - } - - { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(r, c, 2); - EXPECT_NO_THROW(grid.setConcentrations(concentrations)); - - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(r, c + 3, 6); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(r + 3, c, 3); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(r + 2, c + 4, 2); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - } - - { - // correct alpha matrices - MatrixXd alphax = MatrixXd::Constant(r, c, 2); - MatrixXd alphay = MatrixXd::Constant(r, c, 4); - EXPECT_NO_THROW(grid.setAlpha(alphax, alphay)); - - EXPECT_ANY_THROW(grid.setAlpha(alphax)); - - grid.setAlpha(alphax, alphay); - EXPECT_EQ(grid.getAlphaX(), alphax); - EXPECT_EQ(grid.getAlphaY(), alphay); - EXPECT_ANY_THROW(grid.getAlpha()); - - // false alpha matrices - alphax = MatrixXd::Constant(r + 3, c + 1, 3); - EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay)); - alphay = MatrixXd::Constant(r + 2, c + 1, 5); - EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay)); - } - - { - int dr = 8; - int dc = 9; - - // set 1D domain - EXPECT_ANY_THROW(grid.setDomain(dr)); - - // set 2D domain - EXPECT_NO_THROW(grid.setDomain(dr, dc)); - - grid.setDomain(dr, dc); - EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(dc) / double(c)); - EXPECT_DOUBLE_EQ(grid.getDeltaRow(), double(dr) / double(r)); - - // set too small domain - dr = 0; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); - dr = 8; - dc = -1; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); - dr = -2; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); - } - - { - std::vector concentrations(r * c); - - for (int i = 0; i < r * c; i++) { - concentrations[i] = i; - } - - grid.setConcentrations(concentrations.data()); - - EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 0), 0); - EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 1), 1); - EXPECT_DOUBLE_EQ(grid.getConcentrations()(1, 0), c); - } -} diff --git a/test/testSimulation.cpp b/test/testSimulation.cpp deleted file mode 100644 index 6411412..0000000 --- a/test/testSimulation.cpp +++ /dev/null @@ -1,182 +0,0 @@ -#include "TestUtils.hpp" -#include -#include -#include - -#include -#include - -// include the configured header file -#include - -#define DIFFUSION_TEST(x) TEST(Diffusion, x) - -using namespace Eigen; -using namespace std; -using namespace tug; - -Grid64 setupSimulation(double timestep, int iterations) { - int row = 11; - int col = 11; - int domain_row = 10; - int domain_col = 10; - - // Grid - Grid grid = Grid64(row, col); - grid.setDomain(domain_row, domain_col); - - MatrixXd concentrations = MatrixXd::Constant(row, col, 0); - concentrations(5, 5) = 1; - grid.setConcentrations(concentrations); - - MatrixXd alpha = MatrixXd::Constant(row, col, 1); - for (int i = 0; i < 5; i++) { - for (int j = 0; j < 6; j++) { - alpha(i, j) = 0.01; - } - } - for (int i = 0; i < 5; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.001; - } - } - for (int i = 5; i < 11; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.1; - } - } - grid.setAlpha(alpha, alpha); - - return grid; -} - -constexpr double timestep = 0.001; -constexpr double iterations = 7000; - -DIFFUSION_TEST(EqualityFTCS) { - // set string from the header file - string test_path = testSimulationCSVDir; - MatrixXd reference = CSV2Eigen(test_path); - cout << "FTCS Test: " << endl; - - Grid grid = setupSimulation(timestep, iterations); // Boundary - Boundary bc = Boundary(grid); - - // Simulation - - Simulation sim = Simulation(grid, bc); - // sim.setOutputConsole(CONSOLE_OUTPUT_ON); - sim.setTimestep(timestep); - sim.setIterations(iterations); - sim.run(); - - cout << endl; - EXPECT_TRUE(checkSimilarity(reference, grid.getConcentrations(), 0.1)); -} - -DIFFUSION_TEST(EqualityBTCS) { - // set string from the header file - string test_path = testSimulationCSVDir; - MatrixXd reference = CSV2Eigen(test_path); - cout << "BTCS Test: " << endl; - - Grid grid = setupSimulation(timestep, iterations); // Boundary - Boundary bc = Boundary(grid); - - // Simulation - Simulation sim = Simulation(grid, bc); - // sim.setOutputConsole(CONSOLE_OUTPUT_ON); - sim.setTimestep(timestep); - sim.setIterations(iterations); - sim.run(); - - cout << endl; - EXPECT_TRUE(checkSimilarityV2(reference, grid.getConcentrations(), 0.01)); -} - -DIFFUSION_TEST(InitializeEnvironment) { - int rc = 12; - Grid64 grid(rc, rc); - Boundary boundary(grid); - - EXPECT_NO_THROW(Simulation sim(grid, boundary)); -} - -DIFFUSION_TEST(SimulationEnvironment) { - int rc = 12; - Grid64 grid(rc, rc); - Boundary boundary(grid); - Simulation sim(grid, boundary); - - EXPECT_EQ(sim.getIterations(), -1); - - EXPECT_NO_THROW(sim.setIterations(2000)); - EXPECT_EQ(sim.getIterations(), 2000); - EXPECT_THROW(sim.setIterations(-300), std::invalid_argument); - - EXPECT_NO_THROW(sim.setTimestep(0.1)); - EXPECT_DOUBLE_EQ(sim.getTimestep(), 0.1); - EXPECT_THROW(sim.setTimestep(-0.3), std::invalid_argument); -} - -DIFFUSION_TEST(ClosedBoundaries) { - - constexpr std::uint32_t nrows = 5; - constexpr std::uint32_t ncols = 5; - - tug::Grid64 grid(nrows, ncols); - - auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0); - auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); - auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); - - grid.setConcentrations(concentrations); - grid.setAlpha(alphax, alphay); - - tug::Boundary bc(grid); - bc.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1.0); - bc.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1.0); - bc.setBoundarySideConstant(tug::BC_SIDE_TOP, 1.0); - bc.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1.0); - - tug::Simulation sim(grid, bc); - sim.setTimestep(1); - sim.setIterations(1); - - MatrixXd input_values(concentrations); - sim.run(); - - EXPECT_TRUE(checkSimilarityV2(input_values, grid.getConcentrations(), 1E-12)); -} - -DIFFUSION_TEST(ConstantInnerCell) { - constexpr std::uint32_t nrows = 5; - constexpr std::uint32_t ncols = 5; - - tug::Grid64 grid(nrows, ncols); - - auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0); - auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); - auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); - - grid.setConcentrations(concentrations); - grid.setAlpha(alphax, alphay); - - tug::Boundary bc(grid); - // inner - bc.setInnerBoundary(2, 2, 0); - - tug::Simulation sim(grid, bc); - sim.setTimestep(1); - sim.setIterations(1); - - MatrixXd input_values(concentrations); - sim.run(); - - EXPECT_DOUBLE_EQ(grid.getConcentrations()(2, 2), 0); - EXPECT_LT(grid.getConcentrations().sum(), input_values.sum()); - - EXPECT_FALSE((grid.getConcentrations().array() > 1.0).any()); - - EXPECT_FALSE((grid.getConcentrations().array() < 0.0).any()); -} \ No newline at end of file