From 16fbefa9edeb440c5cd14d6f5c6ed50a27145185 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 12 Jun 2024 12:19:46 +0200 Subject: [PATCH 1/4] BREAKING CHANGE: Rename Simulation to Diffusion chore: Cleanup of applications --- docs_sphinx/{Simulation.rst => Diffusion.rst} | 4 +- docs_sphinx/user.rst | 2 +- examples/BTCS_1D_proto_example.cpp | 51 ------ examples/BTCS_2D_proto_example.cpp | 5 +- examples/CMakeLists.txt | 17 +- examples/CRNI_2D_proto_example.cpp | 29 ---- examples/FTCS_1D_proto_example.cpp | 51 ------ examples/FTCS_2D_proto_closed_mdl.cpp | 6 +- examples/FTCS_2D_proto_example.cpp | 92 ---------- examples/FTCS_2D_proto_example_mdl.cpp | 6 +- examples/profiling_openmp.cpp | 70 -------- examples/profiling_speedup.cpp | 70 -------- examples/reference-FTCS_2D_closed.cpp | 53 ------ include/tug/Core/BaseSimulation.hpp | 135 +++++++++++++++ include/tug/Core/{ => Numeric}/BTCS.hpp | 4 +- include/tug/Core/{ => Numeric}/FTCS.hpp | 4 +- include/tug/{Simulation.hpp => Diffusion.hpp} | 157 +++--------------- naaice/BTCS_2D_NAAICE.cpp | 5 +- naaice/NAAICE_dble_vs_float.cpp | 7 +- test/CMakeLists.txt | 2 +- .../{testSimulation.cpp => testDiffusion.cpp} | 17 +- 21 files changed, 184 insertions(+), 603 deletions(-) rename docs_sphinx/{Simulation.rst => Diffusion.rst} (79%) delete mode 100644 examples/BTCS_1D_proto_example.cpp delete mode 100644 examples/CRNI_2D_proto_example.cpp delete mode 100644 examples/FTCS_1D_proto_example.cpp delete mode 100644 examples/FTCS_2D_proto_example.cpp delete mode 100644 examples/profiling_openmp.cpp delete mode 100644 examples/profiling_speedup.cpp delete mode 100644 examples/reference-FTCS_2D_closed.cpp create mode 100644 include/tug/Core/BaseSimulation.hpp rename include/tug/Core/{ => Numeric}/BTCS.hpp (99%) rename include/tug/Core/{ => Numeric}/FTCS.hpp (99%) rename include/tug/{Simulation.hpp => Diffusion.hpp} (74%) rename test/{testSimulation.cpp => testDiffusion.cpp} (92%) 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/user.rst b/docs_sphinx/user.rst index ddd4941..cb9c2c8 100644 --- a/docs_sphinx/user.rst +++ b/docs_sphinx/user.rst @@ -5,4 +5,4 @@ User API 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..9cefd10 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; @@ -61,8 +61,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..a5bd8f4 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; @@ -69,8 +69,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..3f750eb 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; @@ -60,8 +60,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/Core/BaseSimulation.hpp b/include/tug/Core/BaseSimulation.hpp new file mode 100644 index 0000000..5c6ca11 --- /dev/null +++ b/include/tug/Core/BaseSimulation.hpp @@ -0,0 +1,135 @@ +#pragma once + +#include + +namespace tug { + +/** + * @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 */ +}; + +class BaseSimulation { +protected: + CSV_OUTPUT csv_output{CSV_OUTPUT_OFF}; + CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT_OFF}; + TIME_MEASURE time_measure{TIME_MEASURE_OFF}; + + int iterations{1}; + +public: + /** + * @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; + } + + /** + * @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) { + if (console_output < CONSOLE_OUTPUT_OFF && + console_output > CONSOLE_OUTPUT_VERBOSE) { + throw std::invalid_argument("Invalid console output option given!"); + } + + this->console_output = console_output; + } + + /** + * @brief Set the Time Measure option. Off by default. + * + * @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!"); + } + + this->time_measure = 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) { + if (iterations <= 0) { + throw std::invalid_argument( + "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; +}; +} // namespace tug \ No newline at end of file diff --git a/include/tug/Core/BTCS.hpp b/include/tug/Core/Numeric/BTCS.hpp similarity index 99% rename from include/tug/Core/BTCS.hpp rename to include/tug/Core/Numeric/BTCS.hpp index a1081af..3aca4a8 100644 --- a/include/tug/Core/BTCS.hpp +++ b/include/tug/Core/Numeric/BTCS.hpp @@ -10,8 +10,8 @@ #ifndef BTCS_H_ #define BTCS_H_ -#include "Matrix.hpp" -#include "TugUtils.hpp" +#include "../Matrix.hpp" +#include "../TugUtils.hpp" #include #include diff --git a/include/tug/Core/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp similarity index 99% rename from include/tug/Core/FTCS.hpp rename to include/tug/Core/Numeric/FTCS.hpp index a3012b3..41eeb80 100644 --- a/include/tug/Core/FTCS.hpp +++ b/include/tug/Core/Numeric/FTCS.hpp @@ -8,10 +8,8 @@ #ifndef FTCS_H_ #define FTCS_H_ -#include "TugUtils.hpp" +#include "../TugUtils.hpp" -#include -#include #include #ifdef _OPENMP diff --git a/include/tug/Simulation.hpp b/include/tug/Diffusion.hpp similarity index 74% rename from include/tug/Simulation.hpp rename to include/tug/Diffusion.hpp index 7f8632d..4d2f0f0 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Diffusion.hpp @@ -1,13 +1,12 @@ /** - * @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" @@ -19,9 +18,10 @@ #include #include -#include "Core/BTCS.hpp" -#include "Core/FTCS.hpp" +#include "Core/Numeric/BTCS.hpp" +#include "Core/Numeric/FTCS.hpp" #include "Core/TugUtils.hpp" +#include "tug/Core/BaseSimulation.hpp" #ifdef _OPENMP #include @@ -51,37 +51,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,7 +63,17 @@ enum TIME_MEASURE { */ template -class Simulation { +class Diffusion : public BaseSimulation { +private: + T timestep{-1}; + int innerIterations{1}; + int numThreads{omp_get_num_procs()}; + + Grid &grid; + Boundary &bc; + + const std::vector approach_names = {"FTCS", "BTCS", "CRNI"}; + public: /** * @brief Set up a simulation environment. The timestep and number of @@ -108,68 +87,7 @@ public: * @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; - } - - /** - * @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) { - if (console_output < CONSOLE_OUTPUT_OFF && - console_output > CONSOLE_OUTPUT_VERBOSE) { - throw std::invalid_argument("Invalid console output option given!"); - } - - this->console_output = console_output; - } - - /** - * @brief Set the Time Measure option. Off by default. - * - * @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!"); - } - - this->time_measure = time_measure; - } + Diffusion(Grid &_grid, Boundary &_bc) : grid(_grid), bc(_bc){}; /** * @brief Setting the time step for each iteration step. Time step must be @@ -244,20 +162,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,13 +181,6 @@ 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. * @@ -393,7 +290,6 @@ public: auto begin = std::chrono::high_resolution_clock::now(); if constexpr (approach == FTCS_APPROACH) { // FTCS case - for (int i = 0; i < iterations * innerIterations; i++) { if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { printConcentrationsConsole(); @@ -485,20 +381,5 @@ public: << 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/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp index 9ed7ccb..f344893 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" @@ -141,7 +140,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..c432f0f 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; @@ -142,8 +142,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 49b3d73..6a374c9 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,6 +1,6 @@ find_package(doctest REQUIRED) -add_executable(testTug setup.cpp testSimulation.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp) +add_executable(testTug setup.cpp testDiffusion.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp) target_link_libraries(testTug doctest::doctest tug) # get relative path of the CSV file diff --git a/test/testSimulation.cpp b/test/testDiffusion.cpp similarity index 92% rename from test/testSimulation.cpp rename to test/testDiffusion.cpp index 7e6e30f..0aca414 100644 --- a/test/testSimulation.cpp +++ b/test/testDiffusion.cpp @@ -1,9 +1,8 @@ #include "TestUtils.hpp" -#include +#include #include #include -#include #include // include the configured header file @@ -62,7 +61,7 @@ TEST_CASE("equality to reference matrix with FTCS") { // Simulation - Simulation sim = Simulation(grid, bc); + Diffusion sim(grid, bc); // sim.setOutputConsole(CONSOLE_OUTPUT_ON); sim.setTimestep(timestep); sim.setIterations(iterations); @@ -82,7 +81,7 @@ TEST_CASE("equality to reference matrix with BTCS") { Boundary bc = Boundary(grid); // Simulation - Simulation sim = Simulation(grid, bc); + Diffusion sim(grid, bc); // sim.setOutputConsole(CONSOLE_OUTPUT_ON); sim.setTimestep(timestep); sim.setIterations(iterations); @@ -97,16 +96,16 @@ TEST_CASE("Initialize environment") { Grid64 grid(rc, rc); Boundary boundary(grid); - CHECK_NOTHROW(Simulation sim(grid, boundary)); + CHECK_NOTHROW(Diffusion sim(grid, boundary)); } TEST_CASE("Simulation environment") { int rc = 12; Grid64 grid(rc, rc); Boundary boundary(grid); - Simulation sim(grid, boundary); + Diffusion sim(grid, boundary); - SUBCASE("default paremeters") { CHECK_EQ(sim.getIterations(), -1); } + SUBCASE("default paremeters") { CHECK_EQ(sim.getIterations(), 1); } SUBCASE("set iterations") { CHECK_NOTHROW(sim.setIterations(2000)); @@ -141,7 +140,7 @@ TEST_CASE("Closed Boundaries - no change expected") { bc.setBoundarySideConstant(tug::BC_SIDE_TOP, 1.0); bc.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1.0); - tug::Simulation sim(grid, bc); + tug::Diffusion sim(grid, bc); sim.setTimestep(1); sim.setIterations(1); @@ -169,7 +168,7 @@ TEST_CASE("Constant inner cell - 'absorbing' concentration") { // inner bc.setInnerBoundary(2, 2, 0); - tug::Simulation sim(grid, bc); + tug::Diffusion sim(grid, bc); sim.setTimestep(1); sim.setIterations(1); From 9aafb3e18411ff915595d99476932235705209f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 12 Jun 2024 13:55:20 +0200 Subject: [PATCH 2/4] refactor: Use assert instead of custom throw for invalid argument in TugUtils.hpp --- include/tug/Core/TugUtils.hpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) 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_ From 8c18e82b155bf8e511cd838146d0df145dfd712e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Thu, 13 Jun 2024 09:16:51 +0200 Subject: [PATCH 3/4] refactor: Update testDiffusion.cpp and Diffusion.hpp Refactor testDiffusion.cpp and Diffusion.hpp to improve code readability and maintainability. Remove unnecessary exception throwing and replace with assert statements for invalid arguments. --- .gitlab-ci.yml | 2 +- examples/BTCS_2D_proto_example.cpp | 3 +- examples/FTCS_2D_proto_closed_mdl.cpp | 3 +- examples/FTCS_2D_proto_example_mdl.cpp | 3 +- include/tug/Boundary.hpp | 126 ++++------ include/tug/Core/Numeric/BTCS.hpp | 6 +- include/tug/Core/Numeric/FTCS.hpp | 35 +-- include/tug/Diffusion.hpp | 16 +- include/tug/Grid.hpp | 324 +++++++++++-------------- naaice/BTCS_2D_NAAICE.cpp | 5 +- naaice/NAAICE_dble_vs_float.cpp | 5 +- test/testBoundary.cpp | 17 +- test/testDiffusion.cpp | 27 +-- test/testGrid.cpp | 200 +++++---------- 14 files changed, 300 insertions(+), 472 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 4f3100c..c2307fe 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -14,7 +14,7 @@ build_release: expire_in: 100s script: - mkdir build && cd build - - cmake -DCMAKE_BUILD_TYPE=Release -DTUG_ENABLE_TESTING=ON .. + - cmake -DCMAKE_BUILD_TYPE=Debug -DTUG_ENABLE_TESTING=ON .. - make -j$(nproc) - cp ../test/FTCS_11_11_7000.csv test/ diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index 9cefd10..c09d93d 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -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 diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index a5bd8f4..69b31c3 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -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 diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp index 3f750eb..41ddb91 100644 --- a/examples/FTCS_2D_proto_example_mdl.cpp +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -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 diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index ebcde98..cdfe3be 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -8,6 +8,7 @@ #define BOUNDARY_H_ #include "Grid.hpp" +#include "tug/Core/TugUtils.hpp" #include #include @@ -161,13 +162,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 +185,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 +209,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 +229,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 +245,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 +287,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 +303,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 +321,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 +368,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 +382,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 +396,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 +416,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 +437,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 +463,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/Numeric/BTCS.hpp b/include/tug/Core/Numeric/BTCS.hpp index 3aca4a8..460119d 100644 --- a/include/tug/Core/Numeric/BTCS.hpp +++ b/include/tug/Core/Numeric/BTCS.hpp @@ -354,15 +354,15 @@ template static void BTCS_1D(Grid &grid, Boundary &bc, T timestep, Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, Eigen::VectorX &b)) { - int length = grid.getLength(); - T sx = timestep / (grid.getDelta() * grid.getDelta()); + int length = grid.getCol(); + T sx = timestep / (grid.getDeltaCol() * grid.getDeltaCol()); Eigen::VectorX concentrations_t1(length); Eigen::SparseMatrix A; Eigen::VectorX b(length); - const auto &alpha = grid.getAlpha(); + const auto &alpha = grid.getAlphaX(); const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp index 41eeb80..9e899aa 100644 --- a/include/tug/Core/Numeric/FTCS.hpp +++ b/include/tug/Core/Numeric/FTCS.hpp @@ -9,7 +9,9 @@ #define FTCS_H_ #include "../TugUtils.hpp" +#include "tug/Core/Matrix.hpp" +#include #include #ifdef _OPENMP @@ -225,6 +227,7 @@ static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { int colMax = grid.getCol(); T deltaCol = grid.getDeltaCol(); + RowMajMat &concentrations_grid = grid.getConcentrations(); // matrix for concentrations at time t+1 RowMajMat concentrations_t1 = RowMajMat::Constant(1, colMax, 0); @@ -234,7 +237,7 @@ static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { // inner cells // independent of boundary condition type for (int col = 1; col < colMax - 1; col++) { - concentrations_t1(row, col) = grid.getConcentrations()(row, col) + + concentrations_t1(row, col) = concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChange(grid, row, col)); } @@ -242,19 +245,20 @@ static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { // left boundary; hold column constant at index 0 int col = 0; concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(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) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeRightBoundary(grid, bc, row, col)); // overwrite obsolete concentrations - grid.setConcentrations(concentrations_t1); + std::memcpy(concentrations_grid.data(), concentrations_t1.data(), + colMax * sizeof(T)); } // FTCS solution for 2D grid @@ -266,6 +270,8 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, T deltaRow = grid.getDeltaRow(); T deltaCol = grid.getDeltaCol(); + RowMajMat &concentrations_grid = grid.getConcentrations(); + // matrix for concentrations at time t+1 RowMajMat concentrations_t1 = RowMajMat::Constant(rowMax, colMax, 0); @@ -275,7 +281,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, #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) + + concentrations_t1(row, col) = concentrations_grid(row, col) + timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col)) + timestep / (deltaCol * deltaCol) * @@ -290,7 +296,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, #pragma omp parallel for num_threads(numThreads) for (int row = 1; row < rowMax - 1; row++) { concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col)); @@ -302,7 +308,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, #pragma omp parallel for num_threads(numThreads) for (int row = 1; row < rowMax - 1; row++) { concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeRightBoundary(grid, bc, row, col)) + timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col)); @@ -314,7 +320,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, #pragma omp parallel for num_threads(numThreads) for (int col = 1; col < colMax - 1; col++) { concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaRow * deltaRow) * (calcVerticalChangeTopBoundary(grid, bc, row, col)) + timestep / (deltaCol * deltaCol) * @@ -327,7 +333,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, #pragma omp parallel for num_threads(numThreads) for (int col = 1; col < colMax - 1; col++) { concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaRow * deltaRow) * (calcVerticalChangeBottomBoundary(grid, bc, row, col)) + timestep / (deltaCol * deltaCol) * @@ -339,7 +345,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, row = 0; col = 0; concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + timestep / (deltaRow * deltaRow) * @@ -350,7 +356,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, row = 0; col = colMax - 1; concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeRightBoundary(grid, bc, row, col)) + timestep / (deltaRow * deltaRow) * @@ -361,7 +367,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, row = rowMax - 1; col = 0; concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + timestep / (deltaRow * deltaRow) * @@ -372,14 +378,15 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, row = rowMax - 1; col = colMax - 1; concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(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); + std::memcpy(concentrations_grid.data(), concentrations_t1.data(), + rowMax * colMax * sizeof(T)); // } } diff --git a/include/tug/Diffusion.hpp b/include/tug/Diffusion.hpp index 4d2f0f0..675605c 100644 --- a/include/tug/Diffusion.hpp +++ b/include/tug/Diffusion.hpp @@ -96,17 +96,15 @@ 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."); - } + 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) { - const T deltaSquare = grid.getDelta(); - const T maxAlpha = grid.getAlpha().maxCoeff(); + const T deltaSquare = grid.getDeltaCol(); + const T maxAlpha = grid.getAlphaX().maxCoeff(); // Courant-Friedrichs-Lewy condition cfl = deltaSquare / (4 * maxAlpha); @@ -272,12 +270,8 @@ public: * 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!"); - } + tug_assert(this->timestep > 0, "Timestep is not set!"); + tug_assert(this->iterations > 0, "Number of iterations are not set!"); std::string filename; if (this->console_output > CONSOLE_OUTPUT_OFF) { diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index b8e471a..f7ee8ac 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -1,5 +1,4 @@ -#ifndef GRID_H_ -#define GRID_H_ +#pragma once /** * @file Grid.hpp @@ -9,11 +8,12 @@ */ #include "Core/Matrix.hpp" +#include "tug/Core/TugUtils.hpp" #include #include #include #include -#include +#include namespace tug { @@ -26,83 +26,97 @@ namespace tug { 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. + * @brief Construct a new Grid object. * - * @param length Length of the 1D-Grid. Must be greater than 3. + * Constructs a new Grid object with given concentrations, defined by an + * Eigen::Matrix. The domain length is per default the same as the length. The + * alpha coefficients are set to 1. The dimensions of the grid are determined + * by the given matrix, which can also be an Eigen::Vector for a 1D-Grid. + * + * @param concentrations An Eigen3 MatrixX holding the concentrations. */ - Grid(int length) : col(length), domainCol(length) { - if (length <= 3) { - throw std::invalid_argument( - "Given grid length too small. Must be greater than 3."); + Grid(const RowMajMat &concentrations) { + if (concentrations.rows() == 1) { + this->dim = 1; + this->domainCol = static_cast(concentrations.cols()); + this->deltaCol = static_cast(this->domainCol) / + static_cast(concentrations.cols()); // -> 1 + + this->concentrations = concentrations; + return; } - this->dim = 1; - this->deltaCol = - static_cast(this->domainCol) / static_cast(this->col); // -> 1 + if (concentrations.cols() == 1) { + this->dim = 1; + this->domainCol = static_cast(concentrations.rows()); + this->deltaCol = static_cast(this->domainCol) / + static_cast(concentrations.rows()); // -> 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->concentrations = concentrations.transpose(); + return; } 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->domainRow = static_cast(concentrations.rows()); + this->domainCol = static_cast(concentrations.cols()); + this->deltaRow = static_cast(this->domainRow) / + static_cast(concentrations.rows()); // -> 1 + this->deltaCol = static_cast(this->domainCol) / + static_cast(concentrations.cols()); // -> 1 this->concentrations = concentrations; + // this->alphaX = RowMajMat::Constant(concentrations.rows(), + // concentrations.cols(), + // MAT_INIT_VAL); + // this->alphaY = RowMajMat::Constant(concentrations.rows(), + // concentrations.cols(), + // MAT_INIT_VAL); } /** - * @brief Sets the concentrations matrix for a 1D or 2D-Grid. + * @brief Construct a new Grid object. * - * @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! + * Constructs a new 1D Grid object with given concentrations, defined by a + * pointer to consecutive memory and the length of the array. The domain + * length is per default the same as the count of grid cells (length of + * array). The memory region is mapped internally, changes will affect the + * original array and the memory shall not be freed. There is no check for + * correct memory size! + * + * @param concentrations Pointer to consecutive memory holding concentrations. + * @param length Length of the array/the 1D grid. */ - void setConcentrations(T *concentrations) { - tug::RowMajMatMap map(concentrations, this->row, this->col); - this->concentrations = map; + Grid(T *concentrations, std::size_t length) : dim(1) { + this->domainCol = static_cast(length); // -> 1 + this->deltaCol = + static_cast(this->domainCol) / static_cast(length); // -> 1 + + this->concentrations = RowMajMatMap(concentrations, 1, length); + } + + /** + * @brief Construct a new Grid object. + * + * Constructs a new 2D Grid object with given concentrations, defined by a + * pointer to consecutive memory and the number of rows and columns. The + * domain size is per default the same as the number of rows and columns. The + * memory region is mapped internally, changes will affect the original array + * and the memory shall not be freed. There is no check for correct memory + * size! + * + * @param concentrations Pointer to consecutive memory holding concentrations. + * @param row Number of rows. + * @param col Number of columns. + */ + Grid(T *concentrations, std::size_t row, std::size_t col) : dim(2) { + this->domainRow = static_cast(row); // -> 1 + this->domainCol = static_cast(col); // -> 1 + this->deltaCol = + static_cast(this->domainCol) / static_cast(col); // -> 1 + this->deltaRow = + static_cast(this->domainRow) / static_cast(row); // -> 1 + + this->concentrations = RowMajMatMap(concentrations, row, col); } /** @@ -113,6 +127,17 @@ public: */ auto &getConcentrations() { return this->concentrations; } + void initAlpha() { + this->alphaX = RowMajMat::Constant( + this->concentrations.rows(), this->concentrations.cols(), MAT_INIT_VAL); + if (dim > 1) { + + this->alphaY = + RowMajMat::Constant(this->concentrations.rows(), + this->concentrations.cols(), MAT_INIT_VAL); + } + } + /** * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one * dimensional. @@ -121,15 +146,12 @@ public: * 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!"); - } + tug_assert(dim == 1, + "Grid is not one dimensional, use 2D setter function!"); + + tug_assert( + alpha.rows() == 1 || alpha.cols() == this->concentrations.cols(), + "Given matrix of alpha coefficients mismatch with Grid dimensions!"); this->alphaX = alpha; } @@ -143,11 +165,9 @@ public: * 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!"); - } + tug_assert(dim == 1, + "Grid is not one dimensional, use 2D setter function!"); + RowMajMatMap map(alpha, 1, this->col); this->alphaX = map; } @@ -162,21 +182,21 @@ public: * 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!"); - } + tug_assert(dim == 2, + "Grid is not two dimensional, use 1D setter function!"); + + tug_assert(alphaX.rows() == this->concentrations.rows(), + "Alpha in x-direction " + "has wrong number of rows!"); + tug_assert(alphaX.cols() == this->concentrations.cols(), + "Alpha in x-direction has wrong number of columns!"); + + tug_assert(alphaY.rows() == this->concentrations.rows(), + "Alpha in y-direction " + "has wrong number of rows!"); + + tug_assert(alphaY.cols() == this->concentrations.cols(), + "Alpha in y-direction has wrong number of columns!"); this->alphaX = alphaX; this->alphaY = alphaY; @@ -194,33 +214,14 @@ public: * 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!"); - } + tug_assert(dim == 2, "Grid is not two dimensional, there is no alphaY!"); + 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. @@ -228,12 +229,7 @@ public: * @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()!"); - } - + tug_assert(this->alphaX.size() > 0, "AlphaX is empty!"); return this->alphaX; } @@ -244,11 +240,8 @@ public: * @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()!"); - } + tug_assert(dim == 2, "Grid is not two dimensional, there is no alphaY!"); + tug_assert(this->alphaY.size() > 0, "AlphaY is empty!"); return this->alphaY; } @@ -260,34 +253,19 @@ public: */ 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; } + int getRow() const { return this->concentrations.rows(); } /** * @brief Gets the number of columns of the grid. * * @return Number of columns. */ - int getCol() const { return this->col; } + int getCol() const { return this->concentrations.cols(); } /** * @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional. @@ -295,17 +273,12 @@ public: * @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!"); - } + tug_assert(dim == 1, "Grid is not one dimensional, use 2D domain setter!"); + tug_assert(domainLength > 0, "Given domain length is not positive!"); this->domainCol = domainLength; - this->deltaCol = double(this->domainCol) / double(this->col); + this->deltaCol = + double(this->domainCol) / double(this->concentrations.cols()); } /** @@ -317,35 +290,18 @@ public: * 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!"); - } + tug_assert(dim == 2, "Grid is not two dimensional, use 1D domain setter!"); + tug_assert(domainCol > 0, + "Given domain size in x-direction is not positive!"); + tug_assert(domainRow > 0, + "Given domain size in y-direction 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; + this->deltaRow = + double(this->domainRow) / double(this->concentrations.rows()); + this->deltaCol = + double(this->domainCol) / double(this->concentrations.cols()); } /** @@ -361,23 +317,18 @@ public: * @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!"); - } + tug_assert(dim == 2, "Grid is not two dimensional, 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) + T deltaRow; // delta in y-direction (between rows) RowMajMat concentrations; // Matrix holding grid concentrations RowMajMat alphaX; // Matrix holding alpha coefficients in x-direction @@ -389,4 +340,3 @@ private: 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 f344893..15f199f 100644 --- a/naaice/BTCS_2D_NAAICE.cpp +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -104,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.: diff --git a/naaice/NAAICE_dble_vs_float.cpp b/naaice/NAAICE_dble_vs_float.cpp index c432f0f..2238375 100644 --- a/naaice/NAAICE_dble_vs_float.cpp +++ b/naaice/NAAICE_dble_vs_float.cpp @@ -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(); diff --git a/test/testBoundary.cpp b/test/testBoundary.cpp index 471c985..29d00ae 100644 --- a/test/testBoundary.cpp +++ b/test/testBoundary.cpp @@ -1,3 +1,5 @@ +#include +#include #include #include #include @@ -31,8 +33,10 @@ TEST_CASE("BoundaryElement") { } TEST_CASE("Boundary Class") { - Grid grid1D = Grid64(10); - Grid grid2D = Grid64(10, 12); + Eigen::VectorXd conc(10); + Grid grid1D = Grid64(conc); + Eigen::MatrixXd conc2D(10, 12); + Grid grid2D = Grid64(conc2D); Boundary boundary1D = Boundary(grid1D); Boundary boundary2D = Boundary(grid2D); vector> boundary1DVector(1, BoundaryElement(1.0)); @@ -53,26 +57,21 @@ TEST_CASE("Boundary Class") { CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1); CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CLOSED); - CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_TOP)); - CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_BOTTOM)); CHECK_NOTHROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT)); - CHECK_THROWS(boundary1D.setBoundarySideClosed(BC_SIDE_TOP)); CHECK_NOTHROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); CHECK_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); - CHECK_THROWS(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2)); CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); CHECK_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); CHECK_NOTHROW(boundary1D.setInnerBoundary(0, inner_condition_value)); - CHECK_THROWS(boundary1D.setInnerBoundary(0, 0, inner_condition_value)); CHECK_EQ(boundary1D.getInnerBoundary(0), innerBoundary); CHECK_EQ(boundary1D.getInnerBoundary(1).first, false); } SUBCASE("Boundaries 2D case") { - CHECK_NOTHROW(Boundary boundary(grid1D)); + CHECK_NOTHROW(Boundary boundary(grid2D)); CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10); CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_RIGHT).size(), 10); CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_TOP).size(), 12); @@ -85,13 +84,11 @@ TEST_CASE("Boundary Class") { CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP)); CHECK_NOTHROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); CHECK_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); - CHECK_THROWS(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12)); CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); - CHECK_THROWS(boundary2D.setInnerBoundary(0, inner_condition_value)); CHECK_NOTHROW(boundary2D.setInnerBoundary(0, 1, inner_condition_value)); CHECK_EQ(boundary2D.getInnerBoundary(0, 1), innerBoundary); CHECK_EQ(boundary2D.getInnerBoundary(0, 2).first, false); diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index 0aca414..48632a0 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -19,12 +19,11 @@ Grid64 setupSimulation(double timestep, int iterations) { 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); + + Grid grid = Grid64(concentrations); + grid.setDomain(domain_row, domain_col); MatrixXd alpha = MatrixXd::Constant(row, col, 1); for (int i = 0; i < 5; i++) { @@ -93,7 +92,8 @@ TEST_CASE("equality to reference matrix with BTCS") { TEST_CASE("Initialize environment") { int rc = 12; - Grid64 grid(rc, rc); + Eigen::MatrixXd concentrations(rc, rc); + Grid64 grid(concentrations); Boundary boundary(grid); CHECK_NOTHROW(Diffusion sim(grid, boundary)); @@ -101,7 +101,9 @@ TEST_CASE("Initialize environment") { TEST_CASE("Simulation environment") { int rc = 12; - Grid64 grid(rc, rc); + Eigen::MatrixXd concentrations(rc, rc); + Grid64 grid(concentrations); + grid.initAlpha(); Boundary boundary(grid); Diffusion sim(grid, boundary); @@ -110,13 +112,11 @@ TEST_CASE("Simulation environment") { SUBCASE("set iterations") { CHECK_NOTHROW(sim.setIterations(2000)); CHECK_EQ(sim.getIterations(), 2000); - CHECK_THROWS(sim.setIterations(-300)); } SUBCASE("set timestep") { CHECK_NOTHROW(sim.setTimestep(0.1)); CHECK_EQ(sim.getTimestep(), 0.1); - CHECK_THROWS(sim.setTimestep(-0.3)); } } @@ -125,13 +125,12 @@ TEST_CASE("Closed Boundaries - no change expected") { 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); + tug::Grid64 grid(concentrations); + grid.setAlpha(alphax, alphay); tug::Boundary bc(grid); @@ -155,13 +154,11 @@ TEST_CASE("Constant inner cell - 'absorbing' concentration") { 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); + tug::Grid64 grid(concentrations); grid.setAlpha(alphax, alphay); tug::Boundary bc(grid); @@ -183,4 +180,4 @@ TEST_CASE("Constant inner cell - 'absorbing' concentration") { const bool less_zero = (grid.getConcentrations().array() < 0.0).any(); CHECK(less_zero == false); -} \ No newline at end of file +} diff --git a/test/testGrid.cpp b/test/testGrid.cpp index 9eee265..c580665 100644 --- a/test/testGrid.cpp +++ b/test/testGrid.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -7,50 +8,23 @@ using namespace Eigen; using namespace std; using namespace tug; -TEST_CASE("1D Grid, too small length") { - int l = 2; - CHECK_THROWS(Grid64(l)); -} - -TEST_CASE("2D Grid64, too small side") { - int r = 1; - int c = 4; - CHECK_THROWS(Grid64(r, c)); - - r = 4; - c = 1; - CHECK_THROWS(Grid64(r, c)); -} - TEST_CASE("1D Grid64") { int l = 12; - Grid64 grid(l); + Eigen::VectorXd conc(l); + Grid64 grid(conc); + grid.initAlpha(); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 1); - CHECK_EQ(grid.getLength(), l); + CHECK_EQ(grid.getCol(), l); CHECK_EQ(grid.getCol(), l); CHECK_EQ(grid.getRow(), 1); CHECK_EQ(grid.getConcentrations().rows(), 1); CHECK_EQ(grid.getConcentrations().cols(), l); - CHECK_EQ(grid.getAlpha().rows(), 1); - CHECK_EQ(grid.getAlpha().cols(), l); + CHECK_EQ(grid.getAlphaX().rows(), 1); + CHECK_EQ(grid.getAlphaX().cols(), l); CHECK_EQ(grid.getDeltaCol(), 1); - - CHECK_THROWS(grid.getAlphaX()); - CHECK_THROWS(grid.getAlphaY()); - CHECK_THROWS(grid.getDeltaRow()); - } - - SUBCASE("setting concentrations") { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(1, l, 3); - CHECK_NOTHROW(grid.setConcentrations(concentrations)); - - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); } SUBCASE("setting alpha") { @@ -58,16 +32,8 @@ TEST_CASE("1D Grid64") { MatrixXd alpha = MatrixXd::Constant(1, l, 3); CHECK_NOTHROW(grid.setAlpha(alpha)); - CHECK_THROWS(grid.setAlpha(alpha, alpha)); - grid.setAlpha(alpha); - CHECK_EQ(grid.getAlpha(), alpha); - CHECK_THROWS(grid.getAlphaX()); - CHECK_THROWS(grid.getAlphaY()); - - // false alpha matrix - MatrixXd wAlpha = MatrixXd::Constant(3, l, 2); - CHECK_THROWS(grid.setAlpha(wAlpha)); + CHECK_EQ(grid.getAlphaX(), alpha); } SUBCASE("setting domain") { @@ -75,34 +41,24 @@ TEST_CASE("1D Grid64") { // set 1D domain CHECK_NOTHROW(grid.setDomain(d)); - // set 2D domain - CHECK_THROWS(grid.setDomain(d, d)); - grid.setDomain(d); CHECK_EQ(grid.getDeltaCol(), double(d) / double(l)); - CHECK_THROWS(grid.getDeltaRow()); - - // set too small domain - d = 0; - CHECK_THROWS(grid.setDomain(d)); - d = -2; - CHECK_THROWS(grid.setDomain(d)); } } TEST_CASE("2D Grid64 quadratic") { int rc = 12; - Grid64 grid(rc, rc); + Eigen::MatrixXd conc(rc, rc); + Grid64 grid(conc); + grid.initAlpha(); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 2); - CHECK_THROWS(grid.getLength()); CHECK_EQ(grid.getCol(), rc); CHECK_EQ(grid.getRow(), rc); CHECK_EQ(grid.getConcentrations().rows(), rc); CHECK_EQ(grid.getConcentrations().cols(), rc); - CHECK_THROWS(grid.getAlpha()); CHECK_EQ(grid.getAlphaX().rows(), rc); CHECK_EQ(grid.getAlphaX().cols(), rc); @@ -112,79 +68,44 @@ TEST_CASE("2D Grid64 quadratic") { CHECK_EQ(grid.getDeltaCol(), 1); } - SUBCASE("setting concentrations") { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2); - CHECK_NOTHROW(grid.setConcentrations(concentrations)); - - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(rc, rc + 3, 1); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(rc + 3, rc, 4); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(rc + 2, rc + 4, 6); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - } - SUBCASE("setting alphas") { // correct alpha matrices MatrixXd alphax = MatrixXd::Constant(rc, rc, 2); MatrixXd alphay = MatrixXd::Constant(rc, rc, 4); CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); - CHECK_THROWS(grid.setAlpha(alphax)); - grid.setAlpha(alphax, alphay); CHECK_EQ(grid.getAlphaX(), alphax); CHECK_EQ(grid.getAlphaY(), alphay); - CHECK_THROWS(grid.getAlpha()); - - // false alpha matrices - alphax = MatrixXd::Constant(rc + 3, rc + 1, 3); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); - alphay = MatrixXd::Constant(rc + 2, rc + 1, 3); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); } SUBCASE("setting domain") { int dr = 8; int dc = 9; - // set 1D domain - CHECK_THROWS(grid.setDomain(dr)); - // set 2D domain CHECK_NOTHROW(grid.setDomain(dr, dc)); grid.setDomain(dr, dc); CHECK_EQ(grid.getDeltaCol(), double(dc) / double(rc)); CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc)); - - // set too small domain - dr = 0; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = 8; - dc = 0; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = -2; - CHECK_THROWS(grid.setDomain(dr, dc)); } } TEST_CASE("2D Grid64 non-quadratic") { int r = 12; int c = 15; - Grid64 grid(r, c); + Eigen::MatrixXd conc(r, c); + Grid64 grid(conc); + grid.initAlpha(); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 2); - CHECK_THROWS(grid.getLength()); CHECK_EQ(grid.getCol(), c); CHECK_EQ(grid.getRow(), r); CHECK_EQ(grid.getConcentrations().rows(), r); CHECK_EQ(grid.getConcentrations().cols(), c); - CHECK_THROWS(grid.getAlpha()); CHECK_EQ(grid.getAlphaX().rows(), r); CHECK_EQ(grid.getAlphaX().cols(), c); @@ -194,18 +115,55 @@ TEST_CASE("2D Grid64 non-quadratic") { CHECK_EQ(grid.getDeltaCol(), 1); } - SUBCASE("setting concentrations") { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(r, c, 2); - CHECK_NOTHROW(grid.setConcentrations(concentrations)); + SUBCASE("setting alphas") { + // correct alpha matrices + MatrixXd alphax = MatrixXd::Constant(r, c, 2); + MatrixXd alphay = MatrixXd::Constant(r, c, 4); + CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(r, c + 3, 6); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(r + 3, c, 3); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(r + 2, c + 4, 2); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); + grid.setAlpha(alphax, alphay); + CHECK_EQ(grid.getAlphaX(), alphax); + CHECK_EQ(grid.getAlphaY(), alphay); + } + + SUBCASE("setting domain") { + int dr = 8; + int dc = 9; + + // set 2D domain + CHECK_NOTHROW(grid.setDomain(dr, dc)); + + grid.setDomain(dr, dc); + CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c)); + CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r)); + } +} + +TEST_CASE("2D Grid64 non-quadratic from pointer") { + int r = 4; + int c = 5; + std::vector concentrations(r * c); + + for (int i = 0; i < r * c; i++) { + concentrations[i] = i; + } + Grid64 grid(concentrations.data(), r, c); + grid.initAlpha(); + + SUBCASE("correct construction") { + CHECK_EQ(grid.getDim(), 2); + CHECK_EQ(grid.getCol(), c); + CHECK_EQ(grid.getRow(), r); + + CHECK_EQ(grid.getConcentrations().rows(), r); + CHECK_EQ(grid.getConcentrations().cols(), c); + + CHECK_EQ(grid.getAlphaX().rows(), r); + CHECK_EQ(grid.getAlphaX().cols(), c); + CHECK_EQ(grid.getAlphaY().rows(), r); + CHECK_EQ(grid.getAlphaY().cols(), c); + CHECK_EQ(grid.getDeltaRow(), 1); + CHECK_EQ(grid.getDeltaCol(), 1); } SUBCASE("setting alphas") { @@ -214,55 +172,27 @@ TEST_CASE("2D Grid64 non-quadratic") { MatrixXd alphay = MatrixXd::Constant(r, c, 4); CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); - CHECK_THROWS(grid.setAlpha(alphax)); - grid.setAlpha(alphax, alphay); CHECK_EQ(grid.getAlphaX(), alphax); CHECK_EQ(grid.getAlphaY(), alphay); - CHECK_THROWS(grid.getAlpha()); - - // false alpha matrices - alphax = MatrixXd::Constant(r + 3, c + 1, 3); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); - alphay = MatrixXd::Constant(r + 2, c + 1, 5); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); } SUBCASE("setting domain") { int dr = 8; int dc = 9; - // set 1D domain - CHECK_THROWS(grid.setDomain(dr)); - // set 2D domain CHECK_NOTHROW(grid.setDomain(dr, dc)); grid.setDomain(dr, dc); CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c)); CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r)); - - // set too small domain - dr = 0; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = 8; - dc = -1; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = -2; - CHECK_THROWS(grid.setDomain(dr, dc)); } - SUBCASE("set concentration from pointer") { - std::vector concentrations(r * c); - - for (int i = 0; i < r * c; i++) { - concentrations[i] = i; - } - - grid.setConcentrations(concentrations.data()); - + SUBCASE("correct values") { CHECK_EQ(grid.getConcentrations()(0, 0), 0); CHECK_EQ(grid.getConcentrations()(0, 1), 1); CHECK_EQ(grid.getConcentrations()(1, 0), c); + CHECK_EQ(grid.getConcentrations()(2, 1), 2 * c + 1); } -} +} \ No newline at end of file From 0e840c51b308c88cd2f365ab767499818b4beb3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Mon, 17 Jun 2024 08:47:11 +0200 Subject: [PATCH 4/4] [wip] save state --- examples/BTCS_2D_proto_example.cpp | 2 +- examples/FTCS_2D_proto_closed_mdl.cpp | 2 +- examples/FTCS_2D_proto_example_mdl.cpp | 2 +- include/tug/Boundary.hpp | 8 +- include/tug/Core/Matrix.hpp | 45 +++- include/tug/Core/Numeric/BTCS.hpp | 10 +- include/tug/Core/Numeric/FTCS.hpp | 51 ++-- include/tug/Diffusion.hpp | 8 +- include/tug/Grid.hpp | 342 ------------------------- include/tug/UniformAlpha.hpp | 30 +++ include/tug/UniformGrid.hpp | 117 +++++++++ naaice/BTCS_2D_NAAICE.cpp | 2 +- naaice/NAAICE_dble_vs_float.cpp | 5 +- test/CMakeLists.txt | 3 +- test/testBoundary.cpp | 4 +- test/testDiffusion.cpp | 16 +- test/testGrid.cpp | 274 ++++++++++---------- 17 files changed, 387 insertions(+), 534 deletions(-) delete mode 100644 include/tug/Grid.hpp create mode 100644 include/tug/UniformAlpha.hpp create mode 100644 include/tug/UniformGrid.hpp diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index c09d93d..33e7cbe 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -23,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; - Grid64 grid(concentrations); + UniformGrid64 grid(concentrations); // (optional) set alphas of the grid, e.g.: // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index 69b31c3..83a7466 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -43,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; - Grid64 grid(concentrations); + UniformGrid64 grid(concentrations); // (optional) set alphas of the grid, e.g.: MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp index 41ddb91..e16ee5b 100644 --- a/examples/FTCS_2D_proto_example_mdl.cpp +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -34,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; - Grid64 grid(concentrations); + UniformGrid64 grid(concentrations); // (optional) set alphas of the grid, e.g.: MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index cdfe3be..885f40a 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -7,7 +7,7 @@ #ifndef BOUNDARY_H_ #define BOUNDARY_H_ -#include "Grid.hpp" +#include "UniformGrid.hpp" #include "tug/Core/TugUtils.hpp" #include @@ -115,7 +115,7 @@ public: * * @param length Length of the grid */ - Boundary(std::uint32_t length) : Boundary(Grid(length)){}; + Boundary(std::uint32_t length) : Boundary(UnfiormGrid(length)){}; /** * @brief Creates a boundary object for a 2D grid @@ -124,7 +124,7 @@ public: * @param cols Number of columns of the grid */ Boundary(std::uint32_t rows, std::uint32_t cols) - : Boundary(Grid(rows, cols)){}; + : Boundary(UnfiormGrid(rows, cols)){}; /** * @brief Creates a boundary object based on the passed grid object and @@ -133,7 +133,7 @@ public: * @param grid Grid object on the basis of which the simulation takes place * and from which the dimensions (in 2D case) are taken. */ - Boundary(const Grid &grid) + Boundary(const UnfiormGrid &grid) : dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) { if (this->dim == 1) { this->boundaries = std::vector>>( diff --git a/include/tug/Core/Matrix.hpp b/include/tug/Core/Matrix.hpp index 9b1686e..3c5aa82 100644 --- a/include/tug/Core/Matrix.hpp +++ b/include/tug/Core/Matrix.hpp @@ -1,6 +1,8 @@ #pragma once #include +#include +#include namespace tug { /** @@ -13,9 +15,48 @@ namespace tug { * * @tparam T The type of the matrix elements. */ +// template +// using RowMajMat = +// Eigen::Matrix; + template -using RowMajMat = - Eigen::Matrix; +class RowMajMat + : public Eigen::Matrix { +protected: + std::uint8_t dim; + +public: + RowMajMat(Eigen::Index rows, Eigen::Index cols, T initial_value) : dim(2) { + *this = Eigen::Matrix::Constant(rows, cols, initial_value); + }; + + RowMajMat(Eigen::Index n_cells, T initial_value) : dim(1) { + *this = Eigen::Matrix::Constant(1, n_cells, initial_value); + }; + + /** + * @brief Gets the number of rows of the grid. + * + * @return Number of rows. + */ + Eigen::Index getRow() const { return this->rows(); } + + /** + * @brief Gets the number of columns of the grid. + * + * @return Number of columns. + */ + Eigen::Index getCol() const { return this->cols(); } + + /** + * @brief Gets the dimensions of the grid. + * + * @return Dimensions, either 1 or 2. + */ + int getDim() const { return this->dim; } +}; template using RowMajMatMap = Eigen::Map>; } // namespace tug \ No newline at end of file diff --git a/include/tug/Core/Numeric/BTCS.hpp b/include/tug/Core/Numeric/BTCS.hpp index 460119d..8692e53 100644 --- a/include/tug/Core/Numeric/BTCS.hpp +++ b/include/tug/Core/Numeric/BTCS.hpp @@ -351,7 +351,7 @@ 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(UnfiormGrid &grid, Boundary &bc, T timestep, Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, Eigen::VectorX &b)) { int length = grid.getCol(); @@ -396,7 +396,7 @@ 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(UnfiormGrid &grid, Boundary &bc, T timestep, Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, Eigen::VectorX &b), int numThreads) { @@ -449,7 +449,8 @@ static void BTCS_2D(Grid &grid, Boundary &bc, T timestep, // entry point for EigenLU solver; differentiate between 1D and 2D grid template -void BTCS_LU(Grid &grid, Boundary &bc, T timestep, int numThreads) { +void BTCS_LU(UnfiormGrid &grid, Boundary &bc, T timestep, + int numThreads) { if (grid.getDim() == 1) { BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); } else if (grid.getDim() == 2) { @@ -462,7 +463,8 @@ void BTCS_LU(Grid &grid, Boundary &bc, T timestep, int numThreads) { // entry point for Thomas algorithm solver; differentiate 1D and 2D grid template -void BTCS_Thomas(Grid &grid, Boundary &bc, T timestep, int numThreads) { +void BTCS_Thomas(UnfiormGrid &grid, Boundary &bc, T timestep, + int numThreads) { if (grid.getDim() == 1) { BTCS_1D(grid, bc, timestep, ThomasAlgorithm); } else if (grid.getDim() == 2) { diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp index 9e899aa..9bbea7c 100644 --- a/include/tug/Core/Numeric/FTCS.hpp +++ b/include/tug/Core/Numeric/FTCS.hpp @@ -24,7 +24,7 @@ namespace tug { // calculates horizontal change on one cell independent of boundary type template -static inline T calcHorizontalChange(Grid &grid, int &row, int &col) { +static inline T calcHorizontalChange(UnfiormGrid &grid, int &row, int &col) { return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), grid.getAlphaX()(row, col)) * @@ -41,7 +41,7 @@ static inline T calcHorizontalChange(Grid &grid, int &row, int &col) { // calculates vertical change on one cell independent of boundary type template -static inline T calcVerticalChange(Grid &grid, int &row, int &col) { +static inline T calcVerticalChange(UnfiormGrid &grid, int &row, int &col) { return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), grid.getAlphaY()(row, col)) * @@ -58,7 +58,7 @@ static inline T calcVerticalChange(Grid &grid, int &row, int &col) { // calculates horizontal change on one cell with a constant left boundary template -static inline T calcHorizontalChangeLeftBoundaryConstant(Grid &grid, +static inline T calcHorizontalChangeLeftBoundaryConstant(UnfiormGrid &grid, Boundary &bc, int &row, int &col) { @@ -75,8 +75,8 @@ static inline T calcHorizontalChangeLeftBoundaryConstant(Grid &grid, // calculates horizontal change on one cell with a closed left boundary template -static inline T calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, - int &col) { +static inline T calcHorizontalChangeLeftBoundaryClosed(UnfiormGrid &grid, + int &row, int &col) { return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), grid.getAlphaX()(row, col)) * @@ -86,8 +86,9 @@ static inline T calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, // 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) { +static inline T calcHorizontalChangeLeftBoundary(UnfiormGrid &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) { @@ -99,7 +100,7 @@ static inline T calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, // calculates horizontal change on one cell with a constant right boundary template -static inline T calcHorizontalChangeRightBoundaryConstant(Grid &grid, +static inline T calcHorizontalChangeRightBoundaryConstant(UnfiormGrid &grid, Boundary &bc, int &row, int &col) { @@ -116,8 +117,8 @@ static inline T calcHorizontalChangeRightBoundaryConstant(Grid &grid, // calculates horizontal change on one cell with a closed right boundary template -static inline T calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, - int &col) { +static inline T calcHorizontalChangeRightBoundaryClosed(UnfiormGrid &grid, + int &row, int &col) { return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1), grid.getAlphaX()(row, col)) * @@ -127,7 +128,7 @@ static inline T calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, // checks boundary condition type for a cell on the right edge of grid template -static inline T calcHorizontalChangeRightBoundary(Grid &grid, +static inline T calcHorizontalChangeRightBoundary(UnfiormGrid &grid, Boundary &bc, int &row, int &col) { if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CONSTANT) { @@ -141,7 +142,7 @@ static inline T calcHorizontalChangeRightBoundary(Grid &grid, // calculates vertical change on one cell with a constant top boundary template -static inline T calcVerticalChangeTopBoundaryConstant(Grid &grid, +static inline T calcVerticalChangeTopBoundaryConstant(UnfiormGrid &grid, Boundary &bc, int &row, int &col) { @@ -158,8 +159,8 @@ static inline T calcVerticalChangeTopBoundaryConstant(Grid &grid, // calculates vertical change on one cell with a closed top boundary template -static inline T calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, - int &col) { +static inline T calcVerticalChangeTopBoundaryClosed(UnfiormGrid &grid, + int &row, int &col) { return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), grid.getAlphaY()(row, col)) * @@ -169,8 +170,9 @@ static inline T calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, // 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) { +static inline T calcVerticalChangeTopBoundary(UnfiormGrid &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) { @@ -182,7 +184,7 @@ static inline T calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, // calculates vertical change on one cell with a constant bottom boundary template -static inline T calcVerticalChangeBottomBoundaryConstant(Grid &grid, +static inline T calcVerticalChangeBottomBoundaryConstant(UnfiormGrid &grid, Boundary &bc, int &row, int &col) { @@ -199,8 +201,8 @@ static inline T calcVerticalChangeBottomBoundaryConstant(Grid &grid, // calculates vertical change on one cell with a closed bottom boundary template -static inline T calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, - int &col) { +static inline T calcVerticalChangeBottomBoundaryClosed(UnfiormGrid &grid, + int &row, int &col) { return -(calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row - 1, col)) * @@ -210,8 +212,9 @@ static inline T calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, // 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) { +static inline T calcVerticalChangeBottomBoundary(UnfiormGrid &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) { @@ -223,7 +226,7 @@ static inline T calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, // FTCS solution for 1D grid template -static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { +static void FTCS_1D(UnfiormGrid &grid, Boundary &bc, T timestep) { int colMax = grid.getCol(); T deltaCol = grid.getDeltaCol(); @@ -263,7 +266,7 @@ static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { // FTCS solution for 2D grid template -static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, +static void FTCS_2D(UnfiormGrid &grid, Boundary &bc, T timestep, int numThreads) { int rowMax = grid.getRow(); int colMax = grid.getCol(); @@ -392,7 +395,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, // entry point; differentiate between 1D and 2D grid template -void FTCS(Grid &grid, Boundary &bc, T timestep, int &numThreads) { +void FTCS(UnfiormGrid &grid, Boundary &bc, T timestep, int &numThreads) { if (grid.getDim() == 1) { FTCS_1D(grid, bc, timestep); } else if (grid.getDim() == 2) { diff --git a/include/tug/Diffusion.hpp b/include/tug/Diffusion.hpp index 675605c..8094623 100644 --- a/include/tug/Diffusion.hpp +++ b/include/tug/Diffusion.hpp @@ -9,7 +9,7 @@ #pragma once #include "Boundary.hpp" -#include "Grid.hpp" +#include "UniformGrid.hpp" #include #include #include @@ -69,7 +69,7 @@ private: int innerIterations{1}; int numThreads{omp_get_num_procs()}; - Grid &grid; + UnfiormGrid &grid; Boundary &bc; const std::vector approach_names = {"FTCS", "BTCS", "CRNI"}; @@ -87,7 +87,7 @@ public: * @param bc Valid boundary condition object * @param approach Approach to solving the problem. Either FTCS or BTCS. */ - Diffusion(Grid &_grid, Boundary &_bc) : grid(_grid), bc(_bc){}; + Diffusion(UnfiormGrid &_grid, Boundary &_bc) : grid(_grid), bc(_bc){}; /** * @brief Setting the time step for each iteration step. Time step must be @@ -96,7 +96,7 @@ public: * @param timestep Valid timestep greater than zero. */ void setTimestep(T timestep) { - tug_assert(timestep > 0, "Timestep has to be greater than zero."); + tug_assert(timestep > 0, "Timestep has to be greater than zero."); if constexpr (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp deleted file mode 100644 index f7ee8ac..0000000 --- a/include/tug/Grid.hpp +++ /dev/null @@ -1,342 +0,0 @@ -#pragma once - -/** - * @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 "tug/Core/TugUtils.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 Construct a new Grid object. - * - * Constructs a new Grid object with given concentrations, defined by an - * Eigen::Matrix. The domain length is per default the same as the length. The - * alpha coefficients are set to 1. The dimensions of the grid are determined - * by the given matrix, which can also be an Eigen::Vector for a 1D-Grid. - * - * @param concentrations An Eigen3 MatrixX holding the concentrations. - */ - Grid(const RowMajMat &concentrations) { - if (concentrations.rows() == 1) { - this->dim = 1; - this->domainCol = static_cast(concentrations.cols()); - this->deltaCol = static_cast(this->domainCol) / - static_cast(concentrations.cols()); // -> 1 - - this->concentrations = concentrations; - return; - } - - if (concentrations.cols() == 1) { - this->dim = 1; - this->domainCol = static_cast(concentrations.rows()); - this->deltaCol = static_cast(this->domainCol) / - static_cast(concentrations.rows()); // -> 1 - - this->concentrations = concentrations.transpose(); - return; - } - - this->dim = 2; - this->domainRow = static_cast(concentrations.rows()); - this->domainCol = static_cast(concentrations.cols()); - this->deltaRow = static_cast(this->domainRow) / - static_cast(concentrations.rows()); // -> 1 - this->deltaCol = static_cast(this->domainCol) / - static_cast(concentrations.cols()); // -> 1 - - this->concentrations = concentrations; - // this->alphaX = RowMajMat::Constant(concentrations.rows(), - // concentrations.cols(), - // MAT_INIT_VAL); - // this->alphaY = RowMajMat::Constant(concentrations.rows(), - // concentrations.cols(), - // MAT_INIT_VAL); - } - - /** - * @brief Construct a new Grid object. - * - * Constructs a new 1D Grid object with given concentrations, defined by a - * pointer to consecutive memory and the length of the array. The domain - * length is per default the same as the count of grid cells (length of - * array). The memory region is mapped internally, changes will affect the - * original array and the memory shall not be freed. There is no check for - * correct memory size! - * - * @param concentrations Pointer to consecutive memory holding concentrations. - * @param length Length of the array/the 1D grid. - */ - Grid(T *concentrations, std::size_t length) : dim(1) { - this->domainCol = static_cast(length); // -> 1 - this->deltaCol = - static_cast(this->domainCol) / static_cast(length); // -> 1 - - this->concentrations = RowMajMatMap(concentrations, 1, length); - } - - /** - * @brief Construct a new Grid object. - * - * Constructs a new 2D Grid object with given concentrations, defined by a - * pointer to consecutive memory and the number of rows and columns. The - * domain size is per default the same as the number of rows and columns. The - * memory region is mapped internally, changes will affect the original array - * and the memory shall not be freed. There is no check for correct memory - * size! - * - * @param concentrations Pointer to consecutive memory holding concentrations. - * @param row Number of rows. - * @param col Number of columns. - */ - Grid(T *concentrations, std::size_t row, std::size_t col) : dim(2) { - this->domainRow = static_cast(row); // -> 1 - this->domainCol = static_cast(col); // -> 1 - this->deltaCol = - static_cast(this->domainCol) / static_cast(col); // -> 1 - this->deltaRow = - static_cast(this->domainRow) / static_cast(row); // -> 1 - - this->concentrations = RowMajMatMap(concentrations, row, col); - } - - /** - * @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; } - - void initAlpha() { - this->alphaX = RowMajMat::Constant( - this->concentrations.rows(), this->concentrations.cols(), MAT_INIT_VAL); - if (dim > 1) { - - this->alphaY = - RowMajMat::Constant(this->concentrations.rows(), - this->concentrations.cols(), MAT_INIT_VAL); - } - } - - /** - * @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) { - tug_assert(dim == 1, - "Grid is not one dimensional, use 2D setter function!"); - - tug_assert( - alpha.rows() == 1 || alpha.cols() == this->concentrations.cols(), - "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) { - tug_assert(dim == 1, - "Grid is not one dimensional, 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) { - tug_assert(dim == 2, - "Grid is not two dimensional, use 1D setter function!"); - - tug_assert(alphaX.rows() == this->concentrations.rows(), - "Alpha in x-direction " - "has wrong number of rows!"); - tug_assert(alphaX.cols() == this->concentrations.cols(), - "Alpha in x-direction has wrong number of columns!"); - - tug_assert(alphaY.rows() == this->concentrations.rows(), - "Alpha in y-direction " - "has wrong number of rows!"); - - tug_assert(alphaY.cols() == this->concentrations.cols(), - "Alpha in y-direction has wrong number of columns!"); - - 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) { - tug_assert(dim == 2, "Grid is not two dimensional, there is no alphaY!"); - - 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 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 { - tug_assert(this->alphaX.size() > 0, "AlphaX is empty!"); - 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 { - tug_assert(dim == 2, "Grid is not two dimensional, there is no alphaY!"); - tug_assert(this->alphaY.size() > 0, "AlphaY is empty!"); - - return this->alphaY; - } - - /** - * @brief Gets the dimensions of the grid. - * - * @return Dimensions, either 1 or 2. - */ - int getDim() const { return this->dim; } - - /** - * @brief Gets the number of rows of the grid. - * - * @return Number of rows. - */ - int getRow() const { return this->concentrations.rows(); } - - /** - * @brief Gets the number of columns of the grid. - * - * @return Number of columns. - */ - int getCol() const { return this->concentrations.cols(); } - - /** - * @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) { - tug_assert(dim == 1, "Grid is not one dimensional, use 2D domain setter!"); - tug_assert(domainLength > 0, "Given domain length is not positive!"); - - this->domainCol = domainLength; - this->deltaCol = - double(this->domainCol) / double(this->concentrations.cols()); - } - - /** - * @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) { - tug_assert(dim == 2, "Grid is not two dimensional, use 1D domain setter!"); - tug_assert(domainCol > 0, - "Given domain size in x-direction is not positive!"); - tug_assert(domainRow > 0, - "Given domain size in y-direction is not positive!"); - - this->domainRow = domainRow; - this->domainCol = domainCol; - this->deltaRow = - double(this->domainRow) / double(this->concentrations.rows()); - this->deltaCol = - double(this->domainCol) / double(this->concentrations.cols()); - } - - /** - * @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 { - tug_assert(dim == 2, "Grid is not two dimensional, there is no delta in " - "y-direction!"); - - return this->deltaRow; - } - -private: - 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; // 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 diff --git a/include/tug/UniformAlpha.hpp b/include/tug/UniformAlpha.hpp new file mode 100644 index 0000000..11b7dd2 --- /dev/null +++ b/include/tug/UniformAlpha.hpp @@ -0,0 +1,30 @@ +#pragma once + +/** + * @file Grid.hpp + * @brief API of Grid class, that holds a matrix with concenctrations and a + * respective matrix/matrices of alpha coefficients. + * + */ + +#include "tug/UniformGrid.hpp" +#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 UniformAlpha : public RowMajMat { +public: + UniformAlpha(const UniformGrid &grid, T init_value = 0) + : RowMajMat::Constant(grid.rows(), grid.cols(), init_value) {} + + UniformAlpha(const UniformGrid &grid, T *alpha) + : tug::RowMajMatMap(alpha, grid.rows(), grid.cols()) {} +}; + +} // namespace tug diff --git a/include/tug/UniformGrid.hpp b/include/tug/UniformGrid.hpp new file mode 100644 index 0000000..da32588 --- /dev/null +++ b/include/tug/UniformGrid.hpp @@ -0,0 +1,117 @@ +#pragma once + +/** + * @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 "tug/Core/TugUtils.hpp" +#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 UniformGrid : public RowMajMat { +public: + UniformGrid(Eigen::Index n_cells, T spatial_size, T initial_value) + : RowMajMat(n_cells, initial_value), domainCol(spatial_size) { + this->deltaCol = static_cast(this->domainCol) / static_cast(n_cells); + } + + UniformGrid(Eigen::Index row, Eigen::Index col, T spatial_row_size, + T spatial_col_size, T initial_value) + : RowMajMat(row, col, initial_value), domainRow(spatial_row_size), + domainCol(spatial_col_size) { + this->deltaRow = static_cast(this->domainRow) / static_cast(row); + this->deltaCol = static_cast(this->domainCol) / static_cast(col); + } + + UniformGrid(T *concentrations, Eigen::Index n_cells, T spatial_size) + : RowMajMatMap(concentrations, 1, n_cells), domainCol(spatial_size) { + this->deltaCol = static_cast(this->domainCol) / static_cast(n_cells); + } + + UniformGrid(T *concentrations, Eigen::Index row, Eigen::Index col, + T spatial_row_size, T spatial_col_size) + : RowMajMatMap(concentrations, row, col), domainRow(spatial_row_size), + domainCol(spatial_col_size) { + this->deltaRow = static_cast(this->domainRow) / static_cast(row); + this->deltaCol = static_cast(this->domainCol) / static_cast(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) { + tug_assert(this->dim == 1, + "Grid is not one dimensional, use 2D domain setter!"); + tug_assert(domainLength > 0, "Given domain length is not positive!"); + + this->domainCol = domainLength; + this->deltaCol = double(this->domainCol) / double(this->cols()); + } + + /** + * @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) { + tug_assert(this->dim == 2, + "Grid is not two dimensional, use 1D domain setter!"); + tug_assert(domainCol > 0, + "Given domain size in x-direction is not positive!"); + tug_assert(domainRow > 0, + "Given domain size in y-direction is not positive!"); + + this->domainRow = domainRow; + this->domainCol = domainCol; + this->deltaRow = + double(this->domainRow) / double(this->concentrations.rows()); + this->deltaCol = + double(this->domainCol) / double(this->concentrations.cols()); + } + + /** + * @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 { + tug_assert(this->dim == 2, + "Grid is not two dimensional, there is no delta in " + "y-direction!"); + + return this->deltaRow; + } + +private: + T domainCol; // number of domain columns + T domainRow{0}; // number of domain rows + T deltaCol; // delta in x-direction (between columns) + T deltaRow; // delta in y-direction (between rows) +}; + +using UniformGrid64 = UniformGrid; +using UniformGrid32 = UniformGrid; +} // namespace tug diff --git a/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp index 15f199f..0f4285e 100644 --- a/naaice/BTCS_2D_NAAICE.cpp +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -109,7 +109,7 @@ int main(int argc, char *argv[]) { const auto init_values_vec = CSVToVector(INPUT_CONC_FILE); Eigen::MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col); - Grid64 grid(concentrations); + UniformGrid64 grid(concentrations); grid.setDomain(0.005, 0.01); const double sum_init = concentrations.sum(); diff --git a/naaice/NAAICE_dble_vs_float.cpp b/naaice/NAAICE_dble_vs_float.cpp index 2238375..e279386 100644 --- a/naaice/NAAICE_dble_vs_float.cpp +++ b/naaice/NAAICE_dble_vs_float.cpp @@ -1,3 +1,4 @@ +#include "tug/UniformGrid.hpp" #include #include #include @@ -120,8 +121,8 @@ template int doWork(int ngrid) { Eigen::MatrixX initConc64 = Eigen::MatrixX::Constant(ngrid, ngrid, 0); initConc64(50, 50) = 1E-6; - Grid grid(initConc64); - grid.setDomain(0.1, 0.1); + + UniformGrid grid(initConc64.data(), ngrid, ngrid); const T sum_init64 = initConc64.sum(); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6a374c9..62e7272 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,6 +1,7 @@ find_package(doctest REQUIRED) -add_executable(testTug setup.cpp testDiffusion.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp) +# add_executable(testTug setup.cpp testDiffusion.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp) +add_executable(testTug setup.cpp testGrid.cpp) target_link_libraries(testTug doctest::doctest tug) # get relative path of the CSV file diff --git a/test/testBoundary.cpp b/test/testBoundary.cpp index 29d00ae..1c2f017 100644 --- a/test/testBoundary.cpp +++ b/test/testBoundary.cpp @@ -34,9 +34,9 @@ TEST_CASE("BoundaryElement") { TEST_CASE("Boundary Class") { Eigen::VectorXd conc(10); - Grid grid1D = Grid64(conc); + UnfiormGrid grid1D = UniformGrid64(conc); Eigen::MatrixXd conc2D(10, 12); - Grid grid2D = Grid64(conc2D); + UnfiormGrid grid2D = UniformGrid64(conc2D); Boundary boundary1D = Boundary(grid1D); Boundary boundary2D = Boundary(grid2D); vector> boundary1DVector(1, BoundaryElement(1.0)); diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index 48632a0..39e23fc 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -12,7 +12,7 @@ using namespace Eigen; using namespace std; using namespace tug; -Grid64 setupSimulation(double timestep, int iterations) { +UniformGrid64 setupSimulation(double timestep, int iterations) { int row = 11; int col = 11; int domain_row = 10; @@ -22,7 +22,7 @@ Grid64 setupSimulation(double timestep, int iterations) { MatrixXd concentrations = MatrixXd::Constant(row, col, 0); concentrations(5, 5) = 1; - Grid grid = Grid64(concentrations); + UnfiormGrid grid = UniformGrid64(concentrations); grid.setDomain(domain_row, domain_col); MatrixXd alpha = MatrixXd::Constant(row, col, 1); @@ -55,7 +55,7 @@ TEST_CASE("equality to reference matrix with FTCS") { MatrixXd reference = CSV2Eigen(test_path); cout << "FTCS Test: " << endl; - Grid grid = setupSimulation(timestep, iterations); // Boundary + UnfiormGrid grid = setupSimulation(timestep, iterations); // Boundary Boundary bc = Boundary(grid); // Simulation @@ -76,7 +76,7 @@ TEST_CASE("equality to reference matrix with BTCS") { MatrixXd reference = CSV2Eigen(test_path); cout << "BTCS Test: " << endl; - Grid grid = setupSimulation(timestep, iterations); // Boundary + UnfiormGrid grid = setupSimulation(timestep, iterations); // Boundary Boundary bc = Boundary(grid); // Simulation @@ -93,7 +93,7 @@ TEST_CASE("equality to reference matrix with BTCS") { TEST_CASE("Initialize environment") { int rc = 12; Eigen::MatrixXd concentrations(rc, rc); - Grid64 grid(concentrations); + UniformGrid64 grid(concentrations); Boundary boundary(grid); CHECK_NOTHROW(Diffusion sim(grid, boundary)); @@ -102,7 +102,7 @@ TEST_CASE("Initialize environment") { TEST_CASE("Simulation environment") { int rc = 12; Eigen::MatrixXd concentrations(rc, rc); - Grid64 grid(concentrations); + UniformGrid64 grid(concentrations); grid.initAlpha(); Boundary boundary(grid); Diffusion sim(grid, boundary); @@ -129,7 +129,7 @@ TEST_CASE("Closed Boundaries - no change expected") { auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); - tug::Grid64 grid(concentrations); + tug::UniformGrid64 grid(concentrations); grid.setAlpha(alphax, alphay); @@ -158,7 +158,7 @@ TEST_CASE("Constant inner cell - 'absorbing' concentration") { auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); - tug::Grid64 grid(concentrations); + tug::UniformGrid64 grid(concentrations); grid.setAlpha(alphax, alphay); tug::Boundary bc(grid); diff --git a/test/testGrid.cpp b/test/testGrid.cpp index c580665..c169dc4 100644 --- a/test/testGrid.cpp +++ b/test/testGrid.cpp @@ -1,7 +1,8 @@ +#include "tug/UniformGrid.hpp" #include #include #include -#include +#include #include using namespace Eigen; @@ -11,30 +12,29 @@ using namespace tug; TEST_CASE("1D Grid64") { int l = 12; Eigen::VectorXd conc(l); - Grid64 grid(conc); - grid.initAlpha(); + UniformGrid64 grid(l, l, 1); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 1); CHECK_EQ(grid.getCol(), l); CHECK_EQ(grid.getCol(), l); CHECK_EQ(grid.getRow(), 1); - - CHECK_EQ(grid.getConcentrations().rows(), 1); - CHECK_EQ(grid.getConcentrations().cols(), l); - CHECK_EQ(grid.getAlphaX().rows(), 1); - CHECK_EQ(grid.getAlphaX().cols(), l); CHECK_EQ(grid.getDeltaCol(), 1); + + // CHECK_EQ(grid.getConcentrations().rows(), 1); + // CHECK_EQ(grid.getConcentrations().cols(), l); + // CHECK_EQ(grid.getAlphaX().rows(), 1); + // CHECK_EQ(grid.getAlphaX().cols(), l); } - SUBCASE("setting alpha") { - // correct alpha matrix - MatrixXd alpha = MatrixXd::Constant(1, l, 3); - CHECK_NOTHROW(grid.setAlpha(alpha)); + // SUBCASE("setting alpha") { + // // correct alpha matrix + // MatrixXd alpha = MatrixXd::Constant(1, l, 3); + // CHECK_NOTHROW(grid.setAlpha(alpha)); - grid.setAlpha(alpha); - CHECK_EQ(grid.getAlphaX(), alpha); - } + // grid.setAlpha(alpha); + // CHECK_EQ(grid.getAlphaX(), alpha); + // } SUBCASE("setting domain") { int d = 8; @@ -46,153 +46,153 @@ TEST_CASE("1D Grid64") { } } -TEST_CASE("2D Grid64 quadratic") { - int rc = 12; - Eigen::MatrixXd conc(rc, rc); - Grid64 grid(conc); - grid.initAlpha(); +// TEST_CASE("2D Grid64 quadratic") { +// int rc = 12; +// Eigen::MatrixXd conc(rc, rc); +// Grid64 grid(conc); +// grid.initAlpha(); - SUBCASE("correct construction") { - CHECK_EQ(grid.getDim(), 2); - CHECK_EQ(grid.getCol(), rc); - CHECK_EQ(grid.getRow(), rc); +// SUBCASE("correct construction") { +// CHECK_EQ(grid.getDim(), 2); +// CHECK_EQ(grid.getCol(), rc); +// CHECK_EQ(grid.getRow(), rc); - CHECK_EQ(grid.getConcentrations().rows(), rc); - CHECK_EQ(grid.getConcentrations().cols(), rc); +// CHECK_EQ(grid.getConcentrations().rows(), rc); +// CHECK_EQ(grid.getConcentrations().cols(), rc); - CHECK_EQ(grid.getAlphaX().rows(), rc); - CHECK_EQ(grid.getAlphaX().cols(), rc); - CHECK_EQ(grid.getAlphaY().rows(), rc); - CHECK_EQ(grid.getAlphaY().cols(), rc); - CHECK_EQ(grid.getDeltaRow(), 1); - CHECK_EQ(grid.getDeltaCol(), 1); - } +// CHECK_EQ(grid.getAlphaX().rows(), rc); +// CHECK_EQ(grid.getAlphaX().cols(), rc); +// CHECK_EQ(grid.getAlphaY().rows(), rc); +// CHECK_EQ(grid.getAlphaY().cols(), rc); +// CHECK_EQ(grid.getDeltaRow(), 1); +// CHECK_EQ(grid.getDeltaCol(), 1); +// } - SUBCASE("setting alphas") { - // correct alpha matrices - MatrixXd alphax = MatrixXd::Constant(rc, rc, 2); - MatrixXd alphay = MatrixXd::Constant(rc, rc, 4); - CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); +// SUBCASE("setting alphas") { +// // correct alpha matrices +// MatrixXd alphax = MatrixXd::Constant(rc, rc, 2); +// MatrixXd alphay = MatrixXd::Constant(rc, rc, 4); +// CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); - grid.setAlpha(alphax, alphay); - CHECK_EQ(grid.getAlphaX(), alphax); - CHECK_EQ(grid.getAlphaY(), alphay); - } +// grid.setAlpha(alphax, alphay); +// CHECK_EQ(grid.getAlphaX(), alphax); +// CHECK_EQ(grid.getAlphaY(), alphay); +// } - SUBCASE("setting domain") { - int dr = 8; - int dc = 9; +// SUBCASE("setting domain") { +// int dr = 8; +// int dc = 9; - // set 2D domain - CHECK_NOTHROW(grid.setDomain(dr, dc)); +// // set 2D domain +// CHECK_NOTHROW(grid.setDomain(dr, dc)); - grid.setDomain(dr, dc); - CHECK_EQ(grid.getDeltaCol(), double(dc) / double(rc)); - CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc)); - } -} +// grid.setDomain(dr, dc); +// CHECK_EQ(grid.getDeltaCol(), double(dc) / double(rc)); +// CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc)); +// } +// } -TEST_CASE("2D Grid64 non-quadratic") { - int r = 12; - int c = 15; - Eigen::MatrixXd conc(r, c); - Grid64 grid(conc); - grid.initAlpha(); +// TEST_CASE("2D Grid64 non-quadratic") { +// int r = 12; +// int c = 15; +// Eigen::MatrixXd conc(r, c); +// Grid64 grid(conc); +// grid.initAlpha(); - SUBCASE("correct construction") { - CHECK_EQ(grid.getDim(), 2); - CHECK_EQ(grid.getCol(), c); - CHECK_EQ(grid.getRow(), r); +// SUBCASE("correct construction") { +// CHECK_EQ(grid.getDim(), 2); +// CHECK_EQ(grid.getCol(), c); +// CHECK_EQ(grid.getRow(), r); - CHECK_EQ(grid.getConcentrations().rows(), r); - CHECK_EQ(grid.getConcentrations().cols(), c); +// CHECK_EQ(grid.getConcentrations().rows(), r); +// CHECK_EQ(grid.getConcentrations().cols(), c); - CHECK_EQ(grid.getAlphaX().rows(), r); - CHECK_EQ(grid.getAlphaX().cols(), c); - CHECK_EQ(grid.getAlphaY().rows(), r); - CHECK_EQ(grid.getAlphaY().cols(), c); - CHECK_EQ(grid.getDeltaRow(), 1); - CHECK_EQ(grid.getDeltaCol(), 1); - } +// CHECK_EQ(grid.getAlphaX().rows(), r); +// CHECK_EQ(grid.getAlphaX().cols(), c); +// CHECK_EQ(grid.getAlphaY().rows(), r); +// CHECK_EQ(grid.getAlphaY().cols(), c); +// CHECK_EQ(grid.getDeltaRow(), 1); +// CHECK_EQ(grid.getDeltaCol(), 1); +// } - SUBCASE("setting alphas") { - // correct alpha matrices - MatrixXd alphax = MatrixXd::Constant(r, c, 2); - MatrixXd alphay = MatrixXd::Constant(r, c, 4); - CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); +// SUBCASE("setting alphas") { +// // correct alpha matrices +// MatrixXd alphax = MatrixXd::Constant(r, c, 2); +// MatrixXd alphay = MatrixXd::Constant(r, c, 4); +// CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); - grid.setAlpha(alphax, alphay); - CHECK_EQ(grid.getAlphaX(), alphax); - CHECK_EQ(grid.getAlphaY(), alphay); - } +// grid.setAlpha(alphax, alphay); +// CHECK_EQ(grid.getAlphaX(), alphax); +// CHECK_EQ(grid.getAlphaY(), alphay); +// } - SUBCASE("setting domain") { - int dr = 8; - int dc = 9; +// SUBCASE("setting domain") { +// int dr = 8; +// int dc = 9; - // set 2D domain - CHECK_NOTHROW(grid.setDomain(dr, dc)); +// // set 2D domain +// CHECK_NOTHROW(grid.setDomain(dr, dc)); - grid.setDomain(dr, dc); - CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c)); - CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r)); - } -} +// grid.setDomain(dr, dc); +// CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c)); +// CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r)); +// } +// } -TEST_CASE("2D Grid64 non-quadratic from pointer") { - int r = 4; - int c = 5; - std::vector concentrations(r * c); +// TEST_CASE("2D Grid64 non-quadratic from pointer") { +// int r = 4; +// int c = 5; +// std::vector concentrations(r * c); - for (int i = 0; i < r * c; i++) { - concentrations[i] = i; - } - Grid64 grid(concentrations.data(), r, c); - grid.initAlpha(); +// for (int i = 0; i < r * c; i++) { +// concentrations[i] = i; +// } +// Grid64 grid(concentrations.data(), r, c); +// grid.initAlpha(); - SUBCASE("correct construction") { - CHECK_EQ(grid.getDim(), 2); - CHECK_EQ(grid.getCol(), c); - CHECK_EQ(grid.getRow(), r); +// SUBCASE("correct construction") { +// CHECK_EQ(grid.getDim(), 2); +// CHECK_EQ(grid.getCol(), c); +// CHECK_EQ(grid.getRow(), r); - CHECK_EQ(grid.getConcentrations().rows(), r); - CHECK_EQ(grid.getConcentrations().cols(), c); +// CHECK_EQ(grid.getConcentrations().rows(), r); +// CHECK_EQ(grid.getConcentrations().cols(), c); - CHECK_EQ(grid.getAlphaX().rows(), r); - CHECK_EQ(grid.getAlphaX().cols(), c); - CHECK_EQ(grid.getAlphaY().rows(), r); - CHECK_EQ(grid.getAlphaY().cols(), c); - CHECK_EQ(grid.getDeltaRow(), 1); - CHECK_EQ(grid.getDeltaCol(), 1); - } +// CHECK_EQ(grid.getAlphaX().rows(), r); +// CHECK_EQ(grid.getAlphaX().cols(), c); +// CHECK_EQ(grid.getAlphaY().rows(), r); +// CHECK_EQ(grid.getAlphaY().cols(), c); +// CHECK_EQ(grid.getDeltaRow(), 1); +// CHECK_EQ(grid.getDeltaCol(), 1); +// } - SUBCASE("setting alphas") { - // correct alpha matrices - MatrixXd alphax = MatrixXd::Constant(r, c, 2); - MatrixXd alphay = MatrixXd::Constant(r, c, 4); - CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); +// SUBCASE("setting alphas") { +// // correct alpha matrices +// MatrixXd alphax = MatrixXd::Constant(r, c, 2); +// MatrixXd alphay = MatrixXd::Constant(r, c, 4); +// CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); - grid.setAlpha(alphax, alphay); - CHECK_EQ(grid.getAlphaX(), alphax); - CHECK_EQ(grid.getAlphaY(), alphay); - } +// grid.setAlpha(alphax, alphay); +// CHECK_EQ(grid.getAlphaX(), alphax); +// CHECK_EQ(grid.getAlphaY(), alphay); +// } - SUBCASE("setting domain") { - int dr = 8; - int dc = 9; +// SUBCASE("setting domain") { +// int dr = 8; +// int dc = 9; - // set 2D domain - CHECK_NOTHROW(grid.setDomain(dr, dc)); +// // set 2D domain +// CHECK_NOTHROW(grid.setDomain(dr, dc)); - grid.setDomain(dr, dc); - CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c)); - CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r)); - } +// grid.setDomain(dr, dc); +// CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c)); +// CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r)); +// } - SUBCASE("correct values") { - CHECK_EQ(grid.getConcentrations()(0, 0), 0); - CHECK_EQ(grid.getConcentrations()(0, 1), 1); - CHECK_EQ(grid.getConcentrations()(1, 0), c); - CHECK_EQ(grid.getConcentrations()(2, 1), 2 * c + 1); - } -} \ No newline at end of file +// SUBCASE("correct values") { +// CHECK_EQ(grid.getConcentrations()(0, 0), 0); +// CHECK_EQ(grid.getConcentrations()(0, 1), 1); +// CHECK_EQ(grid.getConcentrations()(1, 0), c); +// CHECK_EQ(grid.getConcentrations()(2, 1), 2 * c + 1); +// } +// } \ No newline at end of file