diff --git a/CMakeLists.txt b/CMakeLists.txt index 11e56ec..78e45c9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,50 +1,102 @@ -#debian stable (currently bullseye) +# debian stable (currently bullseye) cmake_minimum_required(VERSION 3.18) -project(tug CXX) +project( + tug + VERSION 0.4 + LANGUAGES CXX) set(CMAKE_CXX_STANDARD 17) find_package(Eigen3 REQUIRED NO_MODULE) find_package(OpenMP) -# find_package(easy_profiler) -# option(EASY_OPTION_LOG "Verbose easy_profiler" 1) -## SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -mfma") -option(TUG_USE_OPENMP "Compile with OpenMP support" ON) +include(GNUInstallDirs) -set(CMAKE_CXX_FLAGS_GENERICOPT "-O3 -march=native" CACHE STRING - "Flags used by the C++ compiler during opt builds." - FORCE) +# find_package(easy_profiler) option(EASY_OPTION_LOG "Verbose easy_profiler" 1) -set(CMAKE_BUILD_TYPE "${CMAKE_BUILD_TYPE}" CACHE STRING - "Choose the type of build, options are: None Debug Release RelWithDebInfo MinSizeRel GenericOpt." - FORCE) +# SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -mfma") +option(TUG_USE_OPENMP "Compile tug with OpenMP support" ON) -option(TUG_USE_UNSAFE_MATH_OPT - "Use compiler options to break IEEE compliances by +set(CMAKE_CXX_FLAGS_GENERICOPT + "-O3 -march=native" + CACHE STRING "Flags used by the C++ compiler during opt builds." FORCE) + +set(CMAKE_BUILD_TYPE + "${CMAKE_BUILD_TYPE}" + CACHE + STRING + "Choose the type of build, options are: None Debug Release RelWithDebInfo MinSizeRel GenericOpt." + FORCE) + +option( + TUG_USE_UNSAFE_MATH_OPT "Use compiler options to break IEEE compliances by oenabling reordering of instructions when adding/multiplying of floating - points." - OFF) + points." OFF) -if(TUG_USE_UNSAFE_MATH_OPT) - add_compile_options(-ffast-math) + +option(TUG_ENABLE_TESTING "Run tests after succesfull compilation" OFF) + +option(TUG_HANNESPHILIPP_EXAMPLES "Compile example applications" OFF) + +option(TUG_NAAICE_EXAMPLE "Enables NAAICE examples with export of CSV files" + OFF) + +add_library(tug INTERFACE) +target_include_directories( + tug INTERFACE $ + $) + +target_link_libraries(tug INTERFACE Eigen3::Eigen) + +target_compile_features(tug INTERFACE cxx_std_17) + +if(TUG_USE_OPENMP AND OpenMP_CXX_FOUND) + target_link_libraries(tug INTERFACE OpenMP::OpenMP_CXX) endif() -option(TUG_ENABLE_TESTING - "Run tests after succesfull compilation" - OFF) +if(TUG_USE_UNSAFE_MATH_OPT) + target_compile_options(tug INTERFACE -ffast-math) +endif() -option(TUG_ENABLE_EXAMPLES - "Compile example applications" - OFF) +install( + TARGETS tug + EXPORT ${PROJECT_NAME}_Targets + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) -add_subdirectory(src) +include(CMakePackageConfigHelpers) +write_basic_package_version_file( + "tugConfigVersion.cmake" + VERSION ${PROJECT_VERSION} + COMPATIBILITY SameMajorVersion) + +configure_package_config_file( + "${PROJECT_SOURCE_DIR}/cmake/${PROJECT_NAME}Config.cmake.in" + "${PROJECT_BINARY_DIR}/${PROJECT_NAME}Config.cmake" + INSTALL_DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/${PROJECT_NAME}/cmake) + +install( + EXPORT ${PROJECT_NAME}_Targets + FILE ${PROJECT_NAME}Targets.cmake + NAMESPACE ${PROJECT_NAME}:: + DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/${PROJECT_NAME}/cmake) + +install(FILES "${PROJECT_BINARY_DIR}/${PROJECT_NAME}Config.cmake" + "${PROJECT_BINARY_DIR}/${PROJECT_NAME}ConfigVersion.cmake" + DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/${PROJECT_NAME}/cmake) + +install(DIRECTORY ${PROJECT_SOURCE_DIR}/include/tug DESTINATION include) if(TUG_ENABLE_TESTING) add_subdirectory(test) endif() -if(TUG_ENABLE_EXAMPLES) +if(TUG_HANNESPHILIPP_EXAMPLES) add_subdirectory(examples) endif() + +if(TUG_NAAICE_EXAMPLE) + add_subdirectory(naaice) +endif() diff --git a/cmake/tugConfig.cmake.in b/cmake/tugConfig.cmake.in new file mode 100644 index 0000000..9c15f36 --- /dev/null +++ b/cmake/tugConfig.cmake.in @@ -0,0 +1,4 @@ +@PACKAGE_INIT@ + +include("${CMAKE_CURRENT_LIST_DIR}/@PROJECT_NAME@Targets.cmake") +check_required_components("@PROJECT_NAME@") diff --git a/docs_sphinx/Boundary.rst b/docs_sphinx/Boundary.rst index d22e133..7f3ce75 100644 --- a/docs_sphinx/Boundary.rst +++ b/docs_sphinx/Boundary.rst @@ -1,4 +1,8 @@ Boundary ======== -.. doxygenclass:: Boundary \ No newline at end of file +.. doxygenenum:: tug::BC_TYPE +.. doxygenenum:: tug::BC_SIDE + +.. doxygenclass:: tug::Boundary +.. doxygenclass:: tug::BoundaryElement diff --git a/docs_sphinx/Grid.rst b/docs_sphinx/Grid.rst index c680ab4..e4dc497 100644 --- a/docs_sphinx/Grid.rst +++ b/docs_sphinx/Grid.rst @@ -1,4 +1,4 @@ Grid ==== -.. doxygenclass:: Grid \ No newline at end of file +.. doxygenclass:: tug::Grid diff --git a/docs_sphinx/Simulation.rst b/docs_sphinx/Simulation.rst index a31ec9b..a461073 100644 --- a/docs_sphinx/Simulation.rst +++ b/docs_sphinx/Simulation.rst @@ -1,4 +1,10 @@ Simulation ========== -.. doxygenclass:: Simulation \ No newline at end of file +.. doxygenenum:: tug::APPROACH +.. doxygenenum:: tug::SOLVER +.. doxygenenum:: tug::CSV_OUTPUT +.. doxygenenum:: tug::CONSOLE_OUTPUT +.. doxygenenum:: tug::TIME_MEASURE + +.. doxygenclass:: tug::Simulation diff --git a/docs_sphinx/developper.rst b/docs_sphinx/developper.rst deleted file mode 100644 index 6b97cc5..0000000 --- a/docs_sphinx/developper.rst +++ /dev/null @@ -1,2 +0,0 @@ -Developper API -============== \ No newline at end of file diff --git a/docs_sphinx/examples.rst b/docs_sphinx/examples.rst index e9736f2..9de5a8f 100644 --- a/docs_sphinx/examples.rst +++ b/docs_sphinx/examples.rst @@ -11,18 +11,21 @@ Two dimensional grid with constant boundaries and FTCS method ------------------------------------------------------------- **Initialization of the grid** -For example, the initalization of a grid with 20 by 20 cells and a domain size (physical extent of the grid) of -also 20 by 20 length units can be done as follows. The setting of the domain is optional here and is set to the -same size as the number of cells in the standard case. As seen in the code, the cells of the grid are set to an -initial value of 0 and only in the upper left corner (0,0) the starting concentration is set to the value 20. +For example, the initalization of a grid with 20 by 20 cells using double values +and a domain size (physical extent of the grid) of also 20 by 20 length units +can be done as follows. The setting of the domain is optional here and is set to +the same size as the number of cells in the standard case. As seen in the code, +the cells of the grid are set to an initial value of 0 and only in the upper +left corner (0,0) the starting concentration is set to the value 20. .. code-block:: cpp int row = 20 int col = 20; - Grid grid = Grid(row,col); + Grid grid(row,col); grid.setDomain(row, col); MatrixXd concentrations = MatrixXd::Constant(row,col,0); + // or MatrixX concentrations = MatrixX::Constant(row,col,0); concentrations(0,0) = 20; grid.setConcentrations(concentrations); @@ -40,14 +43,17 @@ of the grid are set as constant edges with a concentration of 0. bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); **Setting of the simulation parameters and simulation start** -In the last block, a simulation class is created and the objects of the grid and the boundary conditions are passed. The solution -method is also specified (either FCTS or BTCS). Furthermore, the desired time step and the number of iterations are set. The penultimate -parameter specifies the output of the simulated results in a CSV file. In the present case, the result of each iteration step is written -one below the other into the corresponding CSV file. + +In the last block, a simulation class is created and the objects of the grid and +the boundary conditions are passed. The solution method is also specified +(either FCTS or BTCS). Furthermore, the desired time step and the number of +iterations are set. The penultimate parameter specifies the output of the +simulated results in a CSV file. In the present case, the result of each +iteration step is written one below the other into the corresponding CSV file. .. code-block:: cpp - Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); + Simulation simulation(grid, bc); simulation.setTimestep(0.1); simulation.setIterations(1000); simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); @@ -59,4 +65,4 @@ one below the other into the corresponding CSV file. Setting special boundary conditions on individual cells -------------------------------------------------------- \ No newline at end of file +------------------------------------------------------- diff --git a/docs_sphinx/images/class_diagram.svg b/docs_sphinx/images/class_diagram.svg index 3392a01..d33462e 100644 --- a/docs_sphinx/images/class_diagram.svg +++ b/docs_sphinx/images/class_diagram.svg @@ -1,4 +1,4 @@ -
Grid
Grid
- dim:int
- col:int
- row:int
- domain_col: int
- domain_row: int
- delta_col: double
- delta_row: double
- dim:int...
+ Grid(col: int)
+ Grid(row: int, col: int)
+ setConcentrations(concentrations: MatrixXd)
+ getConcentrations(): MatrixXd
+ setAlpha(alpha: MatrixXd)
+ setAlpha(alpha_x: MatrixXd, alpha_y: MatrixXd)
+ getAlphaX(): MatrixXd
+ getAlphaY(): MatrixXd
+ getDim(): int
+ getRow(): int
+ getCol(): int
+ setDomain(domain_col: int)
+ setDomain(domain_row: int, domain_col: int)
+ getDeltaCol(): double
+ getDeltaRow(): double
+ Grid(col: int)...
Boundary
Boundary
- grid: Grid
- type: BC_TYPE
- left: VectorXd
- right: VectorXd
- top: VectorXd
- bottom: VectorXd
- grid: Grid...
+ Boundary(grid: Grid, type: BC_TYPE)
+ getBoundaryConditionType(): BC_TYPE
+ setBoundaryConditionValue(side: BC_SIDE, value: VectorXd)
+ getBoundaryConditionValue(): VectorXd
+ Boundary(grid: Grid, type: BC_TYPE)...
Simulation
Simulation
- timestep: double
- iterations: int
- csv_output: CSV_OUTPUT
- grid: GRID
- bc: Boundary
- approach: APPROACH
- timestep: double...
+ Simulation(grid: Grid, bc: Boundary, approach: APPROACH
+ setOutputCSV(csv_output: CSV_OUTPUT)
+ setTimestep(timestep: double)
+ getTimestep(): double
+ setIterations(iterations: int)
+ getIterations(): int
+ printConcentrationsConsole()
+ printConcentrationsCSV(ident: string, appendMode: bool)
+ run()
+ Simulation(grid: Grid, bc: Boundary, approach: APPROACH...
FTCS
FTCS
BTCS
BTCS
functionalities
functionalities
TUG API
TUG API
BoundaryElement
BoundaryElement
- type: BC_TYPE
- value: double
- type: BC_TYPE...
+ BoundaryElement()
+ BoundaryElement(value: double)
+ setType(type: BC_TYPE)
+ setValue(value: double)
+ getType(): BC_TYPE
+ getValue(): double
+ BoundaryElement()...Text is not SVG - cannot display \ No newline at end of file +
Grid<T>
Grid<T>
- dim:int
- col:int
- row:int
- domain_col: int
- domain_row: int
- delta_col: T
- delta_row: T
- dim:int...
+ Grid(col: int)
+ Grid(row: int, col: int)
+ setConcentrations(concentrations: MatrixX<T>)
+ getConcentrations(): MatrixX<T>
+ setAlpha(alpha: MatrixX<T>)
+ setAlpha(alpha_x: MatrixX<T>, alpha_y:   MatrixX<T>)
+ getAlphaX(): MatrixX<T>
+ getAlphaY(): MatrixX<T>
+ getDim(): int
+ getRow(): int
+ getCol(): int
+ setDomain(domain_col: int)
+ setDomain(domain_row: int, domain_col: int)
+ getDeltaCol(): T
+ getDeltaRow(): T
+ Grid(col: int)...
Boundary<T>
Boundary<T>
- grid: Grid
- type: BC_TYPE
- left: std::vector<T>
- right: std::vector<T>
- top: std::vector<T>
- bottom: std::vector<T>
- grid: Grid...
+ Boundary(grid: Grid, type: BC_TYPE)
+ getBoundaryConditionType(): BC_TYPE
+ setBoundaryConditionValue(side: BC_SIDE, 
value: std::vector<T>)
+ getBoundaryConditionValue(): std::vector<T>
+ Boundary(grid: Grid, type: BC_TYPE)...
Simulation<T, approach, solver>
Simulation<T, approach, solver>
- timestep: double
- iterations: int
- csv_output: CSV_OUTPUT
- grid: GRID
- bc: Boundary
- timestep: double...
+ Simulation(grid: Grid, bc: Boundary, approach: APPROACH
+ setOutputCSV(csv_output: CSV_OUTPUT)
+ setTimestep(timestep: T)
+ getTimestep(): T
+ setIterations(iterations: int)
+ getIterations(): int
+ printConcentrationsConsole()
+ printConcentrationsCSV(ident: string, appendMode: bool)
+ run()
+ Simulation(grid: Grid, bc: Boundary, approach: APPROACH...
FTCS
FTCS
BTCS
BTCS
functionalities
functionalities
TUG API
TUG API
BoundaryElement<T>
BoundaryElement<T>
- type: BC_TYPE
- value: T
- type: BC_TYPE...
+ BoundaryElement()
+ BoundaryElement(value: T)
+ setType(type: BC_TYPE)
+ setValue(value: T)
+ getType(): BC_TYPE
+ getValue(): T
+ BoundaryElement()...Text is not SVG - cannot display \ No newline at end of file diff --git a/docs_sphinx/theory.rst b/docs_sphinx/theory.rst index 27c8146..4d3468f 100644 --- a/docs_sphinx/theory.rst +++ b/docs_sphinx/theory.rst @@ -12,4 +12,4 @@ Numerical Solver **Backward Time-Centered Space (BTCS) Method** -**Forward Time-Centered Space (BTCS) Method** \ No newline at end of file +**Forward Time-Centered Space (FTCS) Method** diff --git a/docs_sphinx/user.rst b/docs_sphinx/user.rst index 8d73236..ddd4941 100644 --- a/docs_sphinx/user.rst +++ b/docs_sphinx/user.rst @@ -3,7 +3,6 @@ User API .. toctree:: - :maxdepth: 2 Grid Boundary - Simulation \ No newline at end of file + Simulation diff --git a/examples/BTCS_1D_proto_example.cpp b/examples/BTCS_1D_proto_example.cpp index 1a53738..af9888e 100644 --- a/examples/BTCS_1D_proto_example.cpp +++ b/examples/BTCS_1D_proto_example.cpp @@ -2,6 +2,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { // ************** @@ -10,7 +11,7 @@ int main(int argc, char *argv[]) { // create a linear grid with 20 cells int cells = 20; - Grid grid = Grid(cells); + Grid64 grid(cells); MatrixXd concentrations = MatrixXd::Constant(1, 20, 0); concentrations(0, 0) = 2000; @@ -31,8 +32,7 @@ int main(int argc, char *argv[]) { // ************************ // set up a simulation environment - Simulation simulation = - Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + Simulation simulation = Simulation(grid, bc); // grid,boundary // set the timestep of the simulation simulation.setTimestep(0.1); // timestep diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index 9a103cb..98acdb9 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -2,6 +2,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { // EASY_PROFILER_ENABLE; @@ -13,7 +14,7 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int row = 40; int col = 50; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); @@ -61,7 +62,7 @@ int main(int argc, char *argv[]) { // set up a simulation environment Simulation simulation = - Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + Simulation(grid, bc); // grid,boundary // set the timestep of the simulation simulation.setTimestep(0.1); // timestep diff --git a/examples/CRNI_2D_proto_example.cpp b/examples/CRNI_2D_proto_example.cpp index 39262e7..fc3f046 100644 --- a/examples/CRNI_2D_proto_example.cpp +++ b/examples/CRNI_2D_proto_example.cpp @@ -2,11 +2,12 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { int row = 20; int col = 20; - Grid grid(row, col); + Grid64 grid(row, col); MatrixXd concentrations = MatrixXd::Constant(row, col, 0); concentrations(10, 10) = 2000; @@ -18,7 +19,8 @@ int main(int argc, char *argv[]) { bc.setBoundarySideClosed(BC_SIDE_TOP); bc.setBoundarySideClosed(BC_SIDE_BOTTOM); - Simulation simulation = Simulation(grid, bc, CRANK_NICOLSON_APPROACH); + Simulation simulation = + Simulation(grid, bc); simulation.setTimestep(0.1); simulation.setIterations(50); simulation.setOutputCSV(CSV_OUTPUT_XTREME); diff --git a/examples/FTCS_1D_proto_example.cpp b/examples/FTCS_1D_proto_example.cpp index e97f0be..b3dc905 100644 --- a/examples/FTCS_1D_proto_example.cpp +++ b/examples/FTCS_1D_proto_example.cpp @@ -3,6 +3,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { // ************** @@ -11,7 +12,7 @@ int main(int argc, char *argv[]) { // create a linear grid with 20 cells int cells = 20; - Grid grid = Grid(cells); + Grid64 grid(cells); MatrixXd concentrations = MatrixXd::Constant(1, 20, 20); grid.setConcentrations(concentrations); @@ -31,7 +32,7 @@ int main(int argc, char *argv[]) { // set up a simulation environment Simulation simulation = - Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + Simulation(grid, bc); // grid,boundary,simulation-approach // (optional) set the timestep of the simulation simulation.setTimestep(0.1); // timestep diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index 845d36f..c6ed02f 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -12,6 +12,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { @@ -30,7 +31,7 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int n2 = row / 2 - 1; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); @@ -69,7 +70,7 @@ int main(int argc, char *argv[]) { // set up a simulation environment Simulation simulation = - Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + 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 index 5f5e520..9b8987a 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -10,6 +10,7 @@ #include using namespace Eigen; +using namespace tug; // #include // #define EASY_PROFILER_ENABLE ::profiler::setEnabled(true); @@ -23,7 +24,7 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int row = 20; int col = 20; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); @@ -67,7 +68,7 @@ int main(int argc, char *argv[]) { // set up a simulation environment Simulation simulation = - Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + Simulation(grid, bc); // grid,boundary,simulation-approach // set the timestep of the simulation simulation.setTimestep(0.1); // timestep diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp index 2c63c0c..fdf40d1 100644 --- a/examples/FTCS_2D_proto_example_mdl.cpp +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -10,6 +10,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { @@ -21,7 +22,7 @@ int main(int argc, char *argv[]) { int row = 64; int col = 64; int n2 = row / 2 - 1; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); @@ -60,7 +61,7 @@ int main(int argc, char *argv[]) { // set up a simulation environment Simulation simulation = - Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + 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 index 5dc13ef..8ad62e0 100644 --- a/examples/profiling_openmp.cpp +++ b/examples/profiling_openmp.cpp @@ -7,6 +7,7 @@ using namespace Eigen; using namespace std; +using namespace tug; int main(int argc, char *argv[]) { @@ -29,7 +30,7 @@ int main(int argc, char *argv[]) { // myfile << "Iterations: " << iterations[j] << endl; for (int k = 0; k < repetition; k++) { cout << "Wiederholung: " << k << endl; - Grid grid = Grid(n[i], n[i]); + Grid64 grid(n[i], n[i]); grid.setDomain(1, 1); MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); @@ -39,8 +40,7 @@ int main(int argc, char *argv[]) { Boundary bc = Boundary(grid); - Simulation sim = Simulation(grid, bc, BTCS_APPROACH); - sim.setSolver(THOMAS_ALGORITHM_SOLVER); + Simulation sim = Simulation(grid, bc); if (argc == 2) { int numThreads = atoi(argv[1]); diff --git a/examples/profiling_speedup.cpp b/examples/profiling_speedup.cpp index 5dc13ef..8ad62e0 100644 --- a/examples/profiling_speedup.cpp +++ b/examples/profiling_speedup.cpp @@ -7,6 +7,7 @@ using namespace Eigen; using namespace std; +using namespace tug; int main(int argc, char *argv[]) { @@ -29,7 +30,7 @@ int main(int argc, char *argv[]) { // myfile << "Iterations: " << iterations[j] << endl; for (int k = 0; k < repetition; k++) { cout << "Wiederholung: " << k << endl; - Grid grid = Grid(n[i], n[i]); + Grid64 grid(n[i], n[i]); grid.setDomain(1, 1); MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); @@ -39,8 +40,7 @@ int main(int argc, char *argv[]) { Boundary bc = Boundary(grid); - Simulation sim = Simulation(grid, bc, BTCS_APPROACH); - sim.setSolver(THOMAS_ALGORITHM_SOLVER); + Simulation sim = Simulation(grid, bc); if (argc == 2) { int numThreads = atoi(argv[1]); diff --git a/examples/reference-FTCS_2D_closed.cpp b/examples/reference-FTCS_2D_closed.cpp index 340a6fd..f08c689 100644 --- a/examples/reference-FTCS_2D_closed.cpp +++ b/examples/reference-FTCS_2D_closed.cpp @@ -4,6 +4,7 @@ using namespace std; using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { int row = 50; @@ -12,7 +13,7 @@ int main(int argc, char *argv[]) { int domain_col = 10; // Grid - Grid grid = Grid(row, col); + Grid64 grid(row, col); grid.setDomain(domain_row, domain_col); MatrixXd concentrations = MatrixXd::Constant(row, col, 0); @@ -41,7 +42,7 @@ int main(int argc, char *argv[]) { Boundary bc = Boundary(grid); // Simulation - Simulation sim = Simulation(grid, bc, FTCS_APPROACH); + Simulation sim = Simulation(grid, bc); sim.setTimestep(0.001); sim.setIterations(10000); sim.setOutputCSV(CSV_OUTPUT_OFF); diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index b803e75..3a7f9b8 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -9,6 +9,10 @@ #include "Grid.hpp" #include +#include +#include + +namespace tug { /** * @brief Enum defining the two implemented boundary conditions. @@ -26,8 +30,10 @@ enum BC_SIDE { BC_SIDE_LEFT, BC_SIDE_RIGHT, BC_SIDE_TOP, BC_SIDE_BOTTOM }; * This class defines the boundary conditions of individual boundary elements. * These can be flexibly used and combined later in other classes. * The class serves as an auxiliary class for structuring the Boundary class. + * + * @tparam T Data type of the boundary condition element */ -class BoundaryElement { +template class BoundaryElement { public: /** * @brief Construct a new Boundary Element object for the closed case. @@ -35,7 +41,7 @@ public: * BC_TYPE_CLOSED, where the value takes -1 and does not hold any * physical meaning. */ - BoundaryElement(); + BoundaryElement(){}; /** * @brief Construct a new Boundary Element object for the constant case. @@ -45,7 +51,7 @@ public: * @param value Value of the constant concentration to be assumed at the * corresponding boundary element. */ - BoundaryElement(double value); + BoundaryElement(T _value) : value(_value), type(BC_TYPE_CONSTANT) {} /** * @brief Allows changing the boundary type of a corresponding @@ -54,7 +60,7 @@ public: * @param type Type of boundary condition. Either BC_TYPE_CONSTANT or BC_TYPE_CLOSED. */ - void setType(BC_TYPE type); + void setType(BC_TYPE type) { this->type = type; }; /** * @brief Sets the value of a boundary condition for the constant case. @@ -62,34 +68,46 @@ public: * @param value Concentration to be considered constant for the * corresponding boundary element. */ - void setValue(double value); + void setValue(double value) { + if (value < 0) { + throw std::invalid_argument("No negative concentration allowed."); + } + if (type == BC_TYPE_CLOSED) { + throw std::invalid_argument( + "No constant boundary concentrations can be set for closed " + "boundaries. Please change type first."); + } + this->value = value; + } /** * @brief Return the type of the boundary condition, i.e. whether the * boundary is considered closed or constant. * - * @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or + * @return Type of boundary condition, either BC_TYPE_CLOSED or BC_TYPE_CONSTANT. */ - BC_TYPE getType(); + BC_TYPE getType() const { return this->type; } /** * @brief Return the concentration value for the constant boundary condition. * - * @return double Value of the concentration. + * @return Value of the concentration. */ - double getValue(); + T getValue() const { return this->value; } private: BC_TYPE type{BC_TYPE_CLOSED}; - double value{-1}; + T value{-1}; }; /** * This class implements the functionality and management of the boundary * conditions in the grid to be simulated. + * + * @tparam Data type of the boundary condition value */ -class Boundary { +template class Boundary { public: /** * @brief Creates a boundary object based on the passed grid object and @@ -98,7 +116,27 @@ public: * @param grid Grid object on the basis of which the simulation takes place * and from which the dimensions (in 2D case) are taken. */ - Boundary(Grid grid); + Boundary(const Grid &grid) + : dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) { + if (this->dim == 1) { + this->boundaries = std::vector>>( + 2); // in 1D only left and right boundary + + this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); + } else if (this->dim == 2) { + this->boundaries = std::vector>>(4); + + this->boundaries[BC_SIDE_LEFT] = + std::vector>(this->rows, BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT] = + std::vector>(this->rows, BoundaryElement()); + this->boundaries[BC_SIDE_TOP] = + std::vector>(this->cols, BoundaryElement()); + this->boundaries[BC_SIDE_BOTTOM] = + std::vector>(this->cols, BoundaryElement()); + } + } /** * @brief Sets all elements of the specified boundary side to the boundary @@ -106,7 +144,21 @@ public: * * @param side Side to be set to closed, e.g. BC_SIDE_LEFT. */ - void setBoundarySideClosed(BC_SIDE side); + void setBoundarySideClosed(BC_SIDE side) { + if (this->dim == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw std::invalid_argument( + "For the one-dimensional case, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } + + const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; + const int n = is_vertical ? this->rows : this->cols; + + this->boundaries[side] = + std::vector>(n, BoundaryElement()); + } /** * @brief Sets all elements of the specified boundary side to the boundary @@ -117,7 +169,21 @@ public: * @param value Concentration to be set for all elements of the specified * page. */ - void setBoundarySideConstant(BC_SIDE side, double value); + void setBoundarySideConstant(BC_SIDE side, double value) { + if (this->dim == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw std::invalid_argument( + "For the one-dimensional case, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } + + const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; + const int n = is_vertical ? this->rows : this->cols; + + this->boundaries[side] = + std::vector>(n, BoundaryElement(value)); + } /** * @brief Specifically sets the boundary element of the specified side @@ -128,7 +194,14 @@ public: * boundary side. Must index an element of the corresponding * side. */ - void setBoundaryElementClosed(BC_SIDE side, int index); + void setBoundaryElemenClosed(BC_SIDE side, int index) { + // tests whether the index really points to an element of the boundary side. + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument( + "Index is selected either too large or too small."); + } + this->boundaries[side][index].setType(BC_TYPE_CLOSED); + } /** * @brief Specifically sets the boundary element of the specified side @@ -142,7 +215,15 @@ public: * @param value Concentration value to which the boundary element should be set. */ - void setBoundaryElementConstant(BC_SIDE side, int index, double value); + void setBoundaryElementConstant(BC_SIDE side, int index, double value) { + // tests whether the index really points to an element of the boundary side. + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument( + "Index is selected either too large or too small."); + } + this->boundaries[side][index].setType(BC_TYPE_CONSTANT); + this->boundaries[side][index].setValue(value); + } /** * @brief Returns the boundary condition of a specified side as a vector @@ -150,19 +231,41 @@ public: * * @param side Boundary side from which the boundary conditions are to be * returned. - * @return vector Contains the boundary conditions as - * BoundaryElement objects. + * @return Contains the boundary conditions as + * BoundaryElement objects. */ - const std::vector getBoundarySide(BC_SIDE side); + const std::vector> &getBoundarySide(BC_SIDE side) const { + if (this->dim == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw std::invalid_argument( + "For the one-dimensional trap, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } + return this->boundaries[side]; + } /** * @brief Get thes Boundary Side Values as a vector. Value is -1 in case some specific boundary is closed. * * @param side Boundary side for which the values are to be returned. - * @return VectorXd Vector with values as doubles. + * @return Vector with values as T. */ - Eigen::VectorXd getBoundarySideValues(BC_SIDE side); + Eigen::VectorX getBoundarySideValues(BC_SIDE side) const { + const std::size_t length = boundaries[side].size(); + Eigen::VectorX values(length); + + for (int i = 0; i < length; i++) { + if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { + values(i) = -1; + continue; + } + values(i) = getBoundaryElementValue(side, i); + } + + return values; + } /** * @brief Returns the boundary condition of a specified element on a given @@ -172,9 +275,16 @@ public: * @param index Index of the boundary element on the corresponding * boundary side. Must index an element of the corresponding * side. - * @return BoundaryElement Boundary condition as a BoundaryElement object. + * @return Boundary condition as a BoundaryElement + * object. */ - BoundaryElement getBoundaryElement(BC_SIDE side, int index); + BoundaryElement getBoundaryElement(BC_SIDE side, int index) const { + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument( + "Index is selected either too large or too small."); + } + return this->boundaries[side][index]; + } /** * @brief Returns the type of a boundary condition, i.e. either BC_TYPE_CLOSED @@ -184,28 +294,47 @@ public: * @param index Index of the boundary element on the corresponding * boundary side. Must index an element of the corresponding side. - * @return BC_TYPE Boundary Type of the corresponding boundary condition. + * @return Boundary Type of the corresponding boundary condition. */ - BC_TYPE getBoundaryElementType(BC_SIDE side, int index); + BC_TYPE getBoundaryElementType(BC_SIDE side, int index) const { + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument( + "Index is selected either too large or too small."); + } + return this->boundaries[side][index].getType(); + } /** * @brief Returns the concentration value of a corresponding - * BoundaryElement object if it is a constant boundary condition. + * BoundaryElement object if it is a constant boundary condition. * * @param side Boundary side in which the boundary condition value is * located. * @param index Index of the boundary element on the corresponding * boundary side. Must index an element of the corresponding * side. - * @return double Concentration of the corresponding BoundaryElement object. + * @return Concentration of the corresponding BoundaryElement + * object. */ - double getBoundaryElementValue(BC_SIDE side, int index); + T getBoundaryElementValue(BC_SIDE side, int index) const { + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument( + "Index is selected either too large or too small."); + } + if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) { + throw std::invalid_argument( + "A value can only be output if it is a constant boundary condition."); + } + return this->boundaries[side][index].getValue(); + } private: - Grid grid; // Boundary is directly dependent on the dimensions of a predefined + const std::uint8_t dim; + const std::uint32_t cols; + const std::uint32_t rows; - std::vector> + std::vector>> boundaries; // Vector with Boundary Element information }; - +} // namespace tug #endif // BOUNDARY_H_ diff --git a/src/BTCS.cpp b/include/tug/Core/BTCS.hpp similarity index 74% rename from src/BTCS.cpp rename to include/tug/Core/BTCS.hpp index 3b314e7..b844ba4 100644 --- a/src/BTCS.cpp +++ b/include/tug/Core/BTCS.hpp @@ -1,5 +1,5 @@ /** - * @file BTCSv2.cpp + * @file BTCS.hpp * @brief Implementation of heterogenous BTCS (backward time-centered space) * solution of diffusion equation in 1D and 2D space. Internally the * alternating-direction implicit (ADI) method is used. Version 2, because @@ -7,10 +7,11 @@ * */ -#include "Schemes.hpp" +#ifndef BTCS_H_ +#define BTCS_H_ + #include "TugUtils.hpp" -#include #include #include #include @@ -22,6 +23,8 @@ #define omp_get_thread_num() 0 #endif +namespace tug { + // calculates coefficient for boundary in constant case template constexpr std::pair calcBoundaryCoeffConstant(T alpha_center, @@ -45,13 +48,15 @@ constexpr std::pair calcBoundaryCoeffClosed(T alpha_center, T alpha_side, } // creates coefficient matrix for next time step from alphas in x-direction -static Eigen::SparseMatrix -createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector &bcLeft, - std::vector &bcRight, int numCols, - int rowIndex, double sx) { +template +static Eigen::SparseMatrix +createCoeffMatrix(const Eigen::MatrixX &alpha, + const std::vector> &bcLeft, + const std::vector> &bcRight, int numCols, + int rowIndex, T sx) { // square matrix of column^2 dimension for the coefficients - Eigen::SparseMatrix cm(numCols, numCols); + Eigen::SparseMatrix cm(numCols, numCols); cm.reserve(Eigen::VectorXi::Constant(numCols, 3)); // left column @@ -142,14 +147,18 @@ constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc, // creates a solution vector for next time step from the current state of // concentrations -static Eigen::VectorXd createSolutionVector( - Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alphaX, - Eigen::MatrixXd &alphaY, std::vector &bcLeft, - std::vector &bcRight, std::vector &bcTop, - std::vector &bcBottom, int length, int rowIndex, double sx, - double sy) { +template +static Eigen::VectorX +createSolutionVector(const Eigen::MatrixX &concentrations, + const Eigen::MatrixX &alphaX, + const Eigen::MatrixX &alphaY, + const std::vector> &bcLeft, + const std::vector> &bcRight, + const std::vector> &bcTop, + const std::vector> &bcBottom, + int length, int rowIndex, T sx, T sy) { - Eigen::VectorXd sv(length); + Eigen::VectorX sv(length); const std::size_t numRows = concentrations.rows(); // inner rows @@ -235,10 +244,11 @@ static Eigen::VectorXd createSolutionVector( // solver for linear equation system; A corresponds to coefficient matrix, // b to the solution vector // use of EigenLU solver -static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix &A, - Eigen::VectorXd &b) { +template +static Eigen::VectorX EigenLUAlgorithm(Eigen::SparseMatrix &A, + Eigen::VectorX &b) { - Eigen::SparseLU> solver; + Eigen::SparseLU> solver; solver.analyzePattern(A); solver.factorize(A); @@ -248,14 +258,15 @@ static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix &A, // solver for linear equation system; A corresponds to coefficient matrix, // b to the solution vector // implementation of Thomas Algorithm -static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix &A, - Eigen::VectorXd &b) { +template +static Eigen::VectorX ThomasAlgorithm(Eigen::SparseMatrix &A, + Eigen::VectorX &b) { Eigen::Index n = b.size(); - Eigen::VectorXd a_diag(n); - Eigen::VectorXd b_diag(n); - Eigen::VectorXd c_diag(n); - Eigen::VectorXd x_vec = b; + Eigen::VectorX a_diag(n); + Eigen::VectorX b_diag(n); + Eigen::VectorX c_diag(n); + Eigen::VectorX x_vec = b; // Fill diagonals vectors b_diag[0] = A.coeff(0, 0); @@ -269,6 +280,29 @@ static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix &A, a_diag[n - 1] = A.coeff(n - 1, n - 2); b_diag[n - 1] = A.coeff(n - 1, n - 1); + // HACK: write CSV to file +#ifdef WRITE_THOMAS_CSV +#include +#include + static std::uint32_t file_index = 0; + std::string file_name = "Thomas_" + std::to_string(file_index++) + ".csv"; + + std::ofstream out_file; + + out_file.open(file_name, std::ofstream::trunc | std::ofstream::out); + + // print header + out_file << "Aa, Ab, Ac, b\n"; + + // iterate through all elements + for (std::size_t i = 0; i < n; i++) { + out_file << a_diag[i] << ", " << b_diag[i] << ", " << c_diag[i] << ", " + << b[i] << "\n"; + } + + out_file.close(); +#endif + // start solving - c_diag and x_vec are overwritten n--; c_diag[0] /= b_diag[0]; @@ -291,23 +325,24 @@ static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix &A, } // BTCS solution for 1D grid -static void -BTCS_1D(Grid &grid, Boundary &bc, double timestep, - Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix &A, - Eigen::VectorXd &b)) { +template +static void BTCS_1D(Grid &grid, Boundary &bc, T timestep, + Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, + Eigen::VectorX &b)) { int length = grid.getLength(); - double sx = timestep / (grid.getDelta() * grid.getDelta()); + T sx = timestep / (grid.getDelta() * grid.getDelta()); - Eigen::VectorXd concentrations_t1(length); + Eigen::VectorX concentrations_t1(length); - Eigen::SparseMatrix A; - Eigen::VectorXd b(length); + Eigen::SparseMatrix A; + Eigen::VectorX b(length); - Eigen::MatrixXd alpha = grid.getAlpha(); - std::vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - std::vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + const auto &alpha = grid.getAlpha(); - Eigen::MatrixXd concentrations = grid.getConcentrations(); + const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + + Eigen::MatrixX concentrations = grid.getConcentrations(); int rowIndex = 0; A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx); // this is exactly same as in 2D @@ -331,31 +366,32 @@ BTCS_1D(Grid &grid, Boundary &bc, double timestep, } // BTCS solution for 2D grid -static void -BTCS_2D(Grid &grid, Boundary &bc, double timestep, - Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix &A, - Eigen::VectorXd &b), - int numThreads) { +template +static void BTCS_2D(Grid &grid, Boundary &bc, T timestep, + Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, + Eigen::VectorX &b), + int numThreads) { int rowMax = grid.getRow(); int colMax = grid.getCol(); - double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); - double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); + T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); + T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); - Eigen::MatrixXd concentrations_t1 = - Eigen::MatrixXd::Constant(rowMax, colMax, 0); - Eigen::VectorXd row_t1(colMax); + Eigen::MatrixX concentrations_t1 = + Eigen::MatrixX::Constant(rowMax, colMax, 0); + Eigen::VectorX row_t1(colMax); - Eigen::SparseMatrix A; - Eigen::VectorXd b; + Eigen::SparseMatrix A; + Eigen::VectorX b; - Eigen::MatrixXd alphaX = grid.getAlphaX(); - Eigen::MatrixXd alphaY = grid.getAlphaY(); - std::vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - std::vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); - std::vector bcTop = bc.getBoundarySide(BC_SIDE_TOP); - std::vector bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); + auto alphaX = grid.getAlphaX(); + auto alphaY = grid.getAlphaY(); - Eigen::MatrixXd concentrations = grid.getConcentrations(); + const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + const auto &bcTop = bc.getBoundarySide(BC_SIDE_TOP); + const auto &bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); + + Eigen::MatrixX concentrations = grid.getConcentrations(); #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) for (int i = 0; i < rowMax; i++) { @@ -364,8 +400,6 @@ BTCS_2D(Grid &grid, Boundary &bc, double timestep, b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, colMax, i, sx, sy); - Eigen::SparseLU> solver; - row_t1 = solverFunc(A, b); concentrations_t1.row(i) = row_t1; @@ -394,7 +428,8 @@ BTCS_2D(Grid &grid, Boundary &bc, double timestep, } // entry point for EigenLU solver; differentiate between 1D and 2D grid -void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { +template +void BTCS_LU(Grid &grid, Boundary &bc, T timestep, int numThreads) { if (grid.getDim() == 1) { BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); } else if (grid.getDim() == 2) { @@ -406,7 +441,8 @@ void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { } // entry point for Thomas algorithm solver; differentiate 1D and 2D grid -void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) { +template +void BTCS_Thomas(Grid &grid, Boundary &bc, T timestep, int numThreads) { if (grid.getDim() == 1) { BTCS_1D(grid, bc, timestep, ThomasAlgorithm); } else if (grid.getDim() == 2) { @@ -416,3 +452,6 @@ void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) { "Error: Only 1- and 2-dimensional grids are defined!"); } } +} // namespace tug + +#endif // BTCS_H_ diff --git a/src/FTCS.cpp b/include/tug/Core/FTCS.hpp similarity index 80% rename from src/FTCS.cpp rename to include/tug/Core/FTCS.hpp index 375fe50..7f12780 100644 --- a/src/FTCS.cpp +++ b/include/tug/Core/FTCS.hpp @@ -1,11 +1,13 @@ /** - * @file FTCS.cpp + * @file FTCS.hpp * @brief Implementation of heterogenous FTCS (forward time-centered space) * solution of diffusion equation in 1D and 2D space. * */ -#include "Schemes.hpp" +#ifndef FTCS_H_ +#define FTCS_H_ + #include "TugUtils.hpp" #include @@ -18,8 +20,11 @@ #define omp_get_thread_num() 0 #endif +namespace tug { + // calculates horizontal change on one cell independent of boundary type -static inline double calcHorizontalChange(Grid &grid, int &row, int &col) { +template +static inline T calcHorizontalChange(Grid &grid, int &row, int &col) { return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), grid.getAlphaX()(row, col)) * @@ -35,7 +40,8 @@ static inline double calcHorizontalChange(Grid &grid, int &row, int &col) { } // calculates vertical change on one cell independent of boundary type -static inline double calcVerticalChange(Grid &grid, int &row, int &col) { +template +static inline T calcVerticalChange(Grid &grid, int &row, int &col) { return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), grid.getAlphaY()(row, col)) * @@ -51,10 +57,10 @@ static inline double calcVerticalChange(Grid &grid, int &row, int &col) { } // calculates horizontal change on one cell with a constant left boundary -static inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, - int &col) { +template +static inline T calcHorizontalChangeLeftBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, int &col) { return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), grid.getAlphaX()(row, col)) * @@ -68,8 +74,9 @@ static inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, } // calculates horizontal change on one cell with a closed left boundary -static inline double -calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) { +template +static inline T calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, + int &col) { return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), grid.getAlphaX()(row, col)) * @@ -78,8 +85,9 @@ calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) { } // checks boundary condition type for a cell on the left edge of grid -static inline double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, + int &row, int &col) { if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) { return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col); } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CLOSED) { @@ -90,10 +98,10 @@ static inline double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, } // calculates horizontal change on one cell with a constant right boundary -static inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, - int &col) { +template +static inline T calcHorizontalChangeRightBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, int &col) { return 2 * grid.getAlphaX()(row, col) * bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) - @@ -107,8 +115,9 @@ static inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid, } // calculates horizontal change on one cell with a closed right boundary -static inline double -calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) { +template +static inline T calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, + int &col) { return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1), grid.getAlphaX()(row, col)) * @@ -117,8 +126,10 @@ calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) { } // checks boundary condition type for a cell on the right edge of grid -static inline double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T calcHorizontalChangeRightBoundary(Grid &grid, + Boundary &bc, int &row, + int &col) { if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) { return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col); } else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CLOSED) { @@ -129,9 +140,10 @@ static inline double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, } // calculates vertical change on one cell with a constant top boundary -static inline double calcVerticalChangeTopBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, int &col) { +template +static inline T calcVerticalChangeTopBoundaryConstant(Grid &grid, + Boundary &bc, int &row, + int &col) { return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), grid.getAlphaY()(row, col)) * @@ -145,8 +157,9 @@ static inline double calcVerticalChangeTopBoundaryConstant(Grid &grid, } // calculates vertical change on one cell with a closed top boundary -static inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, - int &col) { +template +static inline T calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, + int &col) { return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), grid.getConcentrations()(row, col)) * @@ -155,8 +168,9 @@ static inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, } // checks boundary condition type for a cell on the top edge of grid -static inline double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, + int &row, int &col) { if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) { return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col); } else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) { @@ -167,10 +181,10 @@ static inline double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, } // calculates vertical change on one cell with a constant bottom boundary -static inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, - int &col) { +template +static inline T calcVerticalChangeBottomBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, int &col) { return 2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) - @@ -184,8 +198,9 @@ static inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid, } // calculates vertical change on one cell with a closed bottom boundary -static inline double -calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) { +template +static inline T calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, + int &col) { return -(calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row - 1, col)) * @@ -194,8 +209,9 @@ calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) { } // checks boundary condition type for a cell on the bottom edge of grid -static inline double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, + int &row, int &col) { if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) { return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col); } else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) { @@ -206,12 +222,14 @@ static inline double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, } // FTCS solution for 1D grid -static void FTCS_1D(Grid &grid, Boundary &bc, double ×tep) { +template +static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { int colMax = grid.getCol(); - double deltaCol = grid.getDeltaCol(); + T deltaCol = grid.getDeltaCol(); // matrix for concentrations at time t+1 - Eigen::MatrixXd concentrations_t1 = Eigen::MatrixXd::Constant(1, colMax, 0); + Eigen::MatrixX concentrations_t1 = + Eigen::MatrixX::Constant(1, colMax, 0); // only one row in 1D case -> row constant at index 0 int row = 0; @@ -243,16 +261,17 @@ static void FTCS_1D(Grid &grid, Boundary &bc, double ×tep) { } // FTCS solution for 2D grid -static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep, +template +static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, int numThreads) { int rowMax = grid.getRow(); int colMax = grid.getCol(); - double deltaRow = grid.getDeltaRow(); - double deltaCol = grid.getDeltaCol(); + T deltaRow = grid.getDeltaRow(); + T deltaCol = grid.getDeltaCol(); // matrix for concentrations at time t+1 - Eigen::MatrixXd concentrations_t1 = - Eigen::MatrixXd::Constant(rowMax, colMax, 0); + Eigen::MatrixX concentrations_t1 = + Eigen::MatrixX::Constant(rowMax, colMax, 0); // inner cells // these are independent of the boundary condition type @@ -369,7 +388,8 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep, } // entry point; differentiate between 1D and 2D grid -void FTCS(Grid &grid, Boundary &bc, double ×tep, int &numThreads) { +template +void FTCS(Grid &grid, Boundary &bc, T timestep, int &numThreads) { if (grid.getDim() == 1) { FTCS_1D(grid, bc, timestep); } else if (grid.getDim() == 2) { @@ -379,3 +399,6 @@ void FTCS(Grid &grid, Boundary &bc, double ×tep, int &numThreads) { "Error: Only 1- and 2-dimensional grids are defined!"); } } +} // namespace tug + +#endif // FTCS_H_ diff --git a/src/TugUtils.hpp b/include/tug/Core/TugUtils.hpp similarity index 100% rename from src/TugUtils.hpp rename to include/tug/Core/TugUtils.hpp diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index 6a6f973..dee004a 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -10,8 +10,17 @@ #include #include +#include -class Grid { +namespace tug { + +/** + * @brief Holds a matrix with concenctration and respective matrix/matrices of + * alpha coefficients. + * + * @tparam T Type to be used for matrices, e.g. double or float + */ +template class Grid { public: /** * @brief Constructs a new 1D-Grid object of a given length, which holds a @@ -21,7 +30,19 @@ public: * * @param length Length of the 1D-Grid. Must be greater than 3. */ - Grid(int length); + Grid(int length) : col(length), domainCol(length) { + if (length <= 3) { + throw std::invalid_argument( + "Given grid length too small. Must be greater than 3."); + } + + this->dim = 1; + this->deltaCol = + static_cast(this->domainCol) / static_cast(this->col); // -> 1 + + this->concentrations = Eigen::MatrixX::Constant(1, col, MAT_INIT_VAL); + this->alphaX = Eigen::MatrixX::Constant(1, col, MAT_INIT_VAL); + } /** * @brief Constructs a new 2D-Grid object of given dimensions, which holds a @@ -34,146 +55,282 @@ public: * @param row Length of the 2D-Grid in y-direction. Must be greater than 3. * @param col Length of the 2D-Grid in x-direction. Must be greater than 3. */ - Grid(int row, int col); + Grid(int _row, int _col) + : row(_row), col(_col), domainRow(_row), domainCol(_col) { + if (row <= 3 || col <= 3) { + throw std::invalid_argument( + "Given grid dimensions too small. Must each be greater than 3."); + } + + this->dim = 2; + this->deltaRow = + static_cast(this->domainRow) / static_cast(this->row); // -> 1 + this->deltaCol = + static_cast(this->domainCol) / static_cast(this->col); // -> 1 + + this->concentrations = Eigen::MatrixX::Constant(row, col, MAT_INIT_VAL); + this->alphaX = Eigen::MatrixX::Constant(row, col, MAT_INIT_VAL); + this->alphaY = Eigen::MatrixX::Constant(row, col, MAT_INIT_VAL); + } /** * @brief Sets the concentrations matrix for a 1D or 2D-Grid. * - * @param concentrations An Eigen3 MatrixXd holding the concentrations. Matrix - * must have correct dimensions as defined in row and col. (Or length, in 1D - * case). + * @param concentrations An Eigen3 MatrixX holding the concentrations. + * Matrix must have correct dimensions as defined in row and col. (Or length, + * in 1D case). */ - void setConcentrations(const Eigen::MatrixXd &concentrations); + void setConcentrations(const Eigen::MatrixX &concentrations) { + if (concentrations.rows() != this->row || + concentrations.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of concentrations mismatch with Grid dimensions!"); + } + + this->concentrations = concentrations; + } /** * @brief Gets the concentrations matrix for a Grid. * - * @return MatrixXd An Eigen3 matrix holding the concentrations and having the - * same dimensions as the grid. + * @return An Eigen3 matrix holding the concentrations and having + * the same dimensions as the grid. */ - const Eigen::MatrixXd &getConcentrations() { return this->concentrations; } + const Eigen::MatrixX &getConcentrations() { return this->concentrations; } /** * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one * dimensional. * - * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients. - * Matrix columns must have same size as length of grid. + * @param alpha An Eigen3 MatrixX with 1 row holding the alpha + * coefficients. Matrix columns must have same size as length of grid. */ - void setAlpha(const Eigen::MatrixXd &alpha); + void setAlpha(const Eigen::MatrixX &alpha) { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use 2D setter function!"); + } + if (alpha.rows() != 1 || alpha.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of alpha coefficients mismatch with Grid dimensions!"); + } + + this->alphaX = alpha; + } /** * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two * dimensional. * - * @param alphaX An Eigen3 MatrixXd holding the alpha coefficients in + * @param alphaX An Eigen3 MatrixX holding the alpha coefficients in * x-direction. Matrix must be of same size as the grid. - * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in + * @param alphaY An Eigen3 MatrixX holding the alpha coefficients in * y-direction. Matrix must be of same size as the grid. */ - void setAlpha(const Eigen::MatrixXd &alphaX, const Eigen::MatrixXd &alphaY); + void setAlpha(const Eigen::MatrixX &alphaX, + const Eigen::MatrixX &alphaY) { + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably " + "use 1D setter function!"); + } + if (alphaX.rows() != this->row || alphaX.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of alpha coefficients in x-direction " + "mismatch with GRid dimensions!"); + } + if (alphaY.rows() != this->row || alphaY.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of alpha coefficients in y-direction " + "mismatch with GRid dimensions!"); + } + + this->alphaX = alphaX; + this->alphaY = alphaY; + } /** * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one * dimensional. * - * @return MatrixXd A matrix with 1 row holding the alpha coefficients. + * @return A matrix with 1 row holding the alpha coefficients. */ - const Eigen::MatrixXd &getAlpha(); + const Eigen::MatrixX &getAlpha() const { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use either getAlphaX() or getAlphaY()!"); + } + + return this->alphaX; + } /** * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. * Grid must be two dimensional. * - * @return MatrixXd A matrix holding the alpha coefficients in x-direction. + * @return A matrix holding the alpha coefficients in x-direction. */ - const Eigen::MatrixXd &getAlphaX(); + const Eigen::MatrixX &getAlphaX() const { + + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably use getAlpha()!"); + } + + return this->alphaX; + } /** * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. * Grid must be two dimensional. * - * @return MatrixXd A matrix holding the alpha coefficients in y-direction. + * @return A matrix holding the alpha coefficients in y-direction. */ - const Eigen::MatrixXd &getAlphaY(); + const Eigen::MatrixX &getAlphaY() const { + + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably use getAlpha()!"); + } + + return this->alphaY; + } /** * @brief Gets the dimensions of the grid. * - * @return int Dimensions, either 1 or 2. + * @return Dimensions, either 1 or 2. */ - int getDim(); + int getDim() const { return this->dim; } /** * @brief Gets length of 1D grid. Must be one dimensional grid. * - * @return int Length of 1D grid. + * @return Length of 1D grid. */ - int getLength(); + int getLength() const { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use getRow() or getCol()!"); + } + + return col; + } /** * @brief Gets the number of rows of the grid. * - * @return int Number of rows. + * @return Number of rows. */ - int getRow(); + int getRow() const { return this->row; } /** * @brief Gets the number of columns of the grid. * - * @return int Number of columns. + * @return Number of columns. */ - int getCol(); + int getCol() const { return this->col; } /** * @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional. * * @param domainLength A double value of the domain length. Must be positive. */ - void setDomain(double domainLength); + void setDomain(double domainLength) { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probaly " + "use the 2D domain setter!"); + } + if (domainLength <= 0) { + throw std::invalid_argument("Given domain length is not positive!"); + } + + this->domainCol = domainLength; + this->deltaCol = double(this->domainCol) / double(this->col); + } /** * @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional. * - * @param domainRow A double value of the domain size in y-direction. Must be - * positive. - * @param domainCol A double value of the domain size in x-direction. Must be - * positive. + * @param domainRow A double value of the domain size in y-direction. Must + * be positive. + * @param domainCol A double value of the domain size in x-direction. Must + * be positive. */ - void setDomain(double domainRow, double domainCol); + void setDomain(double domainRow, double domainCol) { + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably " + "use the 1D domain setter!"); + } + if (domainRow <= 0 || domainCol <= 0) { + throw std::invalid_argument("Given domain size is not positive!"); + } + + this->domainRow = domainRow; + this->domainCol = domainCol; + this->deltaRow = double(this->domainRow) / double(this->row); + this->deltaCol = double(this->domainCol) / double(this->col); + } /** * @brief Gets the delta value for 1D-Grid. Grid must be one dimensional. * - * @return double Delta value. + * @return Delta value. */ - double getDelta(); + T getDelta() const { + + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use the 2D delta getters"); + } + + return this->deltaCol; + } /** * @brief Gets the delta value in x-direction. * - * @return double Delta value in x-direction. + * @return Delta value in x-direction. */ - double getDeltaCol(); + T getDeltaCol() const { return this->deltaCol; } /** * @brief Gets the delta value in y-direction. Must be two dimensional grid. * - * @return double Delta value in y-direction. + * @return Delta value in y-direction. */ - double getDeltaRow(); + T getDeltaRow() const { + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, meaning there is no " + "delta in y-direction!"); + } + + return this->deltaRow; + } private: - int col; // number of grid columns - int row{1}; // number of grid rows - int dim; // 1D or 2D - double domainCol; // number of domain columns - double domainRow{0}; // number of domain rows - double deltaCol; // delta in x-direction (between columns) - double deltaRow{0}; // delta in y-direction (between rows) - Eigen::MatrixXd concentrations; // Matrix holding grid concentrations - Eigen::MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction - Eigen::MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction + int col; // number of grid columns + int row{1}; // number of grid rows + int dim; // 1D or 2D + T domainCol; // number of domain columns + T domainRow{0}; // number of domain rows + T deltaCol; // delta in x-direction (between columns) + T deltaRow{0}; // delta in y-direction (between rows) + Eigen::MatrixX concentrations; // Matrix holding grid concentrations + Eigen::MatrixX alphaX; // Matrix holding alpha coefficients in x-direction + Eigen::MatrixX alphaY; // Matrix holding alpha coefficients in y-direction + + static constexpr T MAT_INIT_VAL = 0; }; +using Grid64 = Grid; +using Grid32 = Grid; +} // namespace tug #endif // GRID_H_ diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 0c09a7e..4c0f5c8 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -11,23 +11,34 @@ #include "Boundary.hpp" #include "Grid.hpp" +#include +#include +#include +#include +#include #include #include +#include "Core/BTCS.hpp" +#include "Core/FTCS.hpp" +#include "Core/TugUtils.hpp" + #ifdef _OPENMP #include #else #define omp_get_num_procs() 1 #endif +namespace tug { + /** - * @brief Enum defining the two implemented solution approaches. + * @brief Enum defining the implemented solution approaches. * */ enum APPROACH { - FTCS_APPROACH, // Forward Time-Centered Space - BTCS_APPROACH, // Backward Time-Centered Space solved with EigenLU solver - CRANK_NICOLSON_APPROACH + FTCS_APPROACH, /*!< Forward Time-Centered Space */ + BTCS_APPROACH, /*!< Backward Time-Centered Space */ + CRANK_NICOLSON_APPROACH /*!< Crank-Nicolson method */ }; /** @@ -35,9 +46,9 @@ enum APPROACH { * */ enum SOLVER { - EIGEN_LU_SOLVER, // EigenLU solver - THOMAS_ALGORITHM_SOLVER // Thomas Algorithm solver; more efficient for - // tridiagonal matrices + EIGEN_LU_SOLVER, /*!< EigenLU solver */ + THOMAS_ALGORITHM_SOLVER /*!< Thomas Algorithm solver; more efficient for + tridiagonal matrices */ }; /** @@ -45,11 +56,11 @@ enum SOLVER { * */ 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 + 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 */ }; /** @@ -57,9 +68,9 @@ enum CSV_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 + 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 */ }; /** @@ -67,8 +78,8 @@ enum CONSOLE_OUTPUT { * */ enum TIME_MEASURE { - TIME_MEASURE_OFF, // do not print any time measures - TIME_MEASURE_ON // print time measure after last iteration + TIME_MEASURE_OFF, /*!< do not print any time measures */ + TIME_MEASURE_ON /*!< print time measure after last iteration */ }; /** @@ -76,7 +87,13 @@ enum TIME_MEASURE { * and contains all the methods for controlling the desired parameters, such as * time step, number of simulations, etc. * + * @tparam T the type of the internal data structures for grid, boundary + * condition and timestep + * @tparam approach Set the SLE scheme to be used + * @tparam solver Set the solver to be used */ +template class Simulation { public: /** @@ -91,7 +108,7 @@ public: * @param bc Valid boundary condition object * @param approach Approach to solving the problem. Either FTCS or BTCS. */ - Simulation(Grid &grid, Boundary &bc, APPROACH approach); + Simulation(Grid &_grid, Boundary &_bc) : grid(_grid), bc(_bc){}; /** * @brief Set the option to output the results to a CSV file. Off by default. @@ -107,7 +124,13 @@ public: * - CSV_OUTPUT_XTREME: produce csv output with all * concentration matrices and simulation environment */ - void setOutputCSV(CSV_OUTPUT csv_output); + void setOutputCSV(CSV_OUTPUT csv_output) { + if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { + throw std::invalid_argument("Invalid CSV output option given!"); + } + + this->csv_output = csv_output; + } /** * @brief Set the options for outputting information to the console. Off by @@ -122,7 +145,14 @@ public: * - CONSOLE_OUTPUT_VERBOSE: print all concentration * matrices to console */ - void setOutputConsole(CONSOLE_OUTPUT console_output); + void setOutputConsole(CONSOLE_OUTPUT console_output) { + if (console_output < CONSOLE_OUTPUT_OFF && + console_output > CONSOLE_OUTPUT_VERBOSE) { + throw std::invalid_argument("Invalid console output option given!"); + } + + this->console_output = console_output; + } /** * @brief Set the Time Measure option. Off by default. @@ -133,7 +163,13 @@ public: * - TIME_MEASURE_ON: Time of simulation run is printed to * console */ - void setTimeMeasure(TIME_MEASURE time_measure); + void setTimeMeasure(TIME_MEASURE time_measure) { + if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) { + throw std::invalid_argument("Invalid time measure option given!"); + } + + this->time_measure = time_measure; + } /** * @brief Setting the time step for each iteration step. Time step must be @@ -141,14 +177,72 @@ public: * * @param timestep Valid timestep greater than zero. */ - void setTimestep(double timestep); + void setTimestep(T timestep) { + if (timestep <= 0) { + throw_invalid_argument("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(); + + // Courant-Friedrichs-Lewy condition + cfl = deltaSquare / (4 * maxAlpha); + } else if (grid.getDim() == 2) { + const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + // will be 0 if 1D, else ... + const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); + + const T maxAlpha = + std::min(grid.getAlphaX().maxCoeff(), grid.getAlphaY().maxCoeff()); + + cfl = minDeltaSquare / (4 * maxAlpha); + } + const std::string dim = std::to_string(grid.getDim()) + "D"; + + const std::string &approachPrefix = this->approach_names[approach]; + std::cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl + << std::endl; + std::cout << approachPrefix << "_" << dim + << " :: required dt=" << timestep << std::endl; + + if (timestep > cfl) { + + this->innerIterations = (int)ceil(timestep / cfl); + this->timestep = timestep / (double)innerIterations; + + std::cerr << "Warning :: Timestep was adjusted, because of stability " + "conditions. Time duration was approximately preserved by " + "adjusting internal number of iterations." + << std::endl; + std::cout << approachPrefix << "_" << dim << " :: Required " + << this->innerIterations + << " inner iterations with dt=" << this->timestep + << std::endl; + + } else { + + this->timestep = timestep; + std::cout << approachPrefix << "_" << dim + << " :: No inner iterations required, dt=" << timestep + << std::endl; + } + } else { + this->timestep = timestep; + } + } /** * @brief Currently set time step is returned. * * @return double timestep */ - double getTimestep(); + T getTimestep() const { return this->timestep; } /** * @brief Set the desired iterations to be calculated. A value greater @@ -156,16 +250,13 @@ public: * * @param iterations Number of iterations to be simulated. */ - void setIterations(int iterations); - - /** - * @brief Set the desired linear equation solver to be used for BTCS approach. - * Without effect in case of FTCS approach. - * - * @param solver Solver to be used. Default is Thomas Algorithm as it is more - * efficient for tridiagonal Matrices. - */ - void setSolver(SOLVER solver); + 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. @@ -175,20 +266,32 @@ public: * maximum number of processors is set as the default case during Simulation * construction. */ - void setNumberThreads(int num_threads); + void setNumberThreads(int num_threads) { + if (numThreads > 0 && numThreads <= omp_get_num_procs()) { + this->numThreads = numThreads; + } else { + int maxThreadNumber = omp_get_num_procs(); + throw std::invalid_argument( + "Number of threads exceeds the number of processor cores (" + + std::to_string(maxThreadNumber) + ") or is less than 1."); + } + } /** * @brief Return the currently set iterations to be calculated. * * @return int Number of iterations. */ - int getIterations(); + int getIterations() const { return this->iterations; } /** * @brief Outputs the current concentrations of the grid on the console. * */ - void printConcentrationsConsole(); + inline void printConcentrationsConsole() const { + std::cout << grid.getConcentrations() << std::endl; + std::cout << std::endl; + } /** * @brief Creates a CSV file with a name containing the current simulation @@ -199,7 +302,52 @@ public: * * @return string Filename with configured simulation parameters. */ - std::string createCSVfile(); + std::string createCSVfile() const { + std::ofstream file; + int appendIdent = 0; + std::string appendIdentString; + + // string approachString = (approach == 0) ? "FTCS" : "BTCS"; + const std::string &approachString = this->approach_names[approach]; + std::string row = std::to_string(grid.getRow()); + std::string col = std::to_string(grid.getCol()); + std::string numIterations = std::to_string(iterations); + + std::string filename = + approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; + + while (std::filesystem::exists(filename)) { + appendIdent += 1; + appendIdentString = std::to_string(appendIdent); + filename = approachString + "_" + row + "_" + col + "_" + numIterations + + "-" + appendIdentString + ".csv"; + } + + file.open(filename); + if (!file) { + exit(1); + } + + // adds lines at the beginning of verbose output csv that represent the + // boundary conditions and their values -1 in case of closed boundary + if (csv_output == CSV_OUTPUT_XTREME) { + Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "", + " "); + file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) + << std::endl; // boundary left + file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) + << std::endl; // boundary right + file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) + << std::endl; // boundary top + file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) + << std::endl; // boundary bottom + file << std::endl << std::endl; + } + + file.close(); + + return filename; + } /** * @brief Writes the currently calculated concentration values of the grid @@ -208,16 +356,138 @@ public: * @param filename Name of the file to which the concentration values are * to be written. */ - void printConcentrationsCSV(const std::string &filename); + void printConcentrationsCSV(const std::string &filename) const { + std::ofstream file; + + file.open(filename, std::ios_base::app); + if (!file) { + exit(1); + } + + Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); + file << grid.getConcentrations().format(do_not_align) << std::endl; + file << std::endl << std::endl; + file.close(); + } /** * @brief Method starts the simulation process with the previously set * parameters. */ - void run(); + 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!"); + } + + std::string filename; + if (this->console_output > CONSOLE_OUTPUT_OFF) { + printConcentrationsConsole(); + } + if (this->csv_output > CSV_OUTPUT_OFF) { + filename = createCSVfile(); + } + + 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(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + FTCS(this->grid, this->bc, this->timestep, this->numThreads); + + // if (i % (iterations * innerIterations / 100) == 0) { + // double percentage = (double)i / ((double)iterations * + // (double)innerIterations) * 100; if ((int)percentage % 10 == 0) { + // cout << "Progress: " << percentage << "%" << endl; + // } + // } + } + + } else if constexpr (approach == BTCS_APPROACH) { // BTCS case + + if constexpr (solver == EIGEN_LU_SOLVER) { + for (int i = 0; i < iterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads); + } + } else if constexpr (solver == THOMAS_ALGORITHM_SOLVER) { + for (int i = 0; i < iterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); + } + } + + } else if constexpr (approach == + CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case + + constexpr T beta = 0.5; + + // TODO this implementation is very inefficient! + // a separate implementation that sets up a specific tridiagonal matrix + // for Crank-Nicolson would be better + Eigen::MatrixX concentrations; + Eigen::MatrixX concentrationsFTCS; + Eigen::MatrixX concentrationsResult; + for (int i = 0; i < iterations * innerIterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + concentrations = grid.getConcentrations(); + FTCS(this->grid, this->bc, this->timestep, this->numThreads); + concentrationsFTCS = grid.getConcentrations(); + grid.setConcentrations(concentrations); + BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); + concentrationsResult = + beta * concentrationsFTCS + (1 - beta) * grid.getConcentrations(); + grid.setConcentrations(concentrationsResult); + } + } + + auto end = std::chrono::high_resolution_clock::now(); + auto milliseconds = + std::chrono::duration_cast(end - begin); + + if (this->console_output > CONSOLE_OUTPUT_OFF) { + printConcentrationsConsole(); + } + if (this->csv_output > CSV_OUTPUT_OFF) { + printConcentrationsCSV(filename); + } + if (this->time_measure > TIME_MEASURE_OFF) { + const std::string &approachString = this->approach_names[approach]; + const std::string dimString = std::to_string(grid.getDim()) + "D"; + std::cout << approachString << dimString << ":: run() finished in " + << milliseconds.count() << "ms" << std::endl; + } + } private: - double timestep{-1}; + T timestep{-1}; int iterations{-1}; int innerIterations{1}; int numThreads{omp_get_num_procs()}; @@ -225,12 +495,10 @@ private: CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT_OFF}; TIME_MEASURE time_measure{TIME_MEASURE_OFF}; - Grid &grid; - Boundary &bc; - APPROACH approach; - SOLVER solver{THOMAS_ALGORITHM_SOLVER}; + 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 new file mode 100644 index 0000000..9ed7ccb --- /dev/null +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -0,0 +1,169 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "files.hpp" + +using namespace tug; + +/** + * Try to parse an input string into a given template type. + */ +template inline T parseString(const std::string &str) { + T result; + std::istringstream iss(str); + + if (!(iss >> result)) { + throw std::invalid_argument("Invalid input for parsing."); + } + + return result; +} + +/** + * Splits a given string into a vector by using a delimiter character. + */ +template +std::vector tokenize(const std::string &input, char delimiter) { + std::vector tokens; + std::istringstream tokenStream(input); + std::string token; + + while (std::getline(tokenStream, token, delimiter)) { + tokens.push_back(parseString(token)); + } + + return tokens; +} + +/** + * Opens a file containing CSV and transform it into row-major 2D STL vector. + */ +template +std::vector> CSVToVector(const char *filename) { + std::ifstream in_file(filename); + if (!in_file.is_open()) { + throw std::runtime_error("Error opening file \'" + std::string(filename) + + "\'."); + } + + std::vector> csv_data; + + std::string line; + while (std::getline(in_file, line)) { + csv_data.push_back(tokenize(line, ',')); + } + + in_file.close(); + + return csv_data; +} + +/** + * Converts a 2D STL vector, where values are stored row-major into a + * column-major Eigen::Matrix. + */ +template +Eigen::MatrixXd rmVecTocmMatrix(const std::vector> &vec, + std::uint32_t exp_rows, + std::uint32_t exp_cols) { + if (exp_rows != vec.size()) { + throw std::runtime_error( + "Mismatch in y dimension while converting to Eigen::Matrix."); + } + + Eigen::MatrixXd out_mat(exp_rows, exp_cols); + + for (std::uint32_t ri = 0; ri < exp_rows; ri++) { + const auto &vec_row = vec[ri]; + if (vec[ri].size() != exp_cols) { + throw std::runtime_error( + "Mismatch in x dimension while converting to Eigen::Matrix."); + } + for (std::uint32_t cj = 0; cj < exp_cols; cj++) { + out_mat(ri, cj) = vec_row[cj]; + } + } + + return out_mat; +} + +int main(int argc, char *argv[]) { + // EASY_PROFILER_ENABLE; + // profiler::startListen(); + // ************** + // **** GRID **** + // ************** + // profiler::startListen(); + // 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); + + const double sum_init = concentrations.sum(); + + // // (optional) set alphas of the grid, e.g.: + const auto alphax_vec = CSVToVector(INPUT_ALPHAX_FILE); + Eigen::MatrixXd alphax = rmVecTocmMatrix(alphax_vec, row, col); + + constexpr double alphay_val = 5e-10; + Eigen::MatrixXd alphay = + Eigen::MatrixXd::Constant(row, col, alphay_val); // row,col,value + grid.setAlpha(alphax, alphay); + + // // ****************** + // // **** BOUNDARY **** + // // ****************** + + // create a boundary with constant values + 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 ENV **** + // // ************************ + + // set up a simulation environment + Simulation simulation = Simulation(grid, bc); // grid,boundary + + // set the timestep of the simulation + simulation.setTimestep(360); // timestep + + // set the number of iterations + simulation.setIterations(1); + + // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_ON); + + // set output to the console to 'ON' + simulation.setOutputConsole(CONSOLE_OUTPUT_ON); + + // // **** RUN SIMULATION **** + simulation.run(); + + const double sum_after = grid.getConcentrations().sum(); + + std::cout << "Sum of init field: " << std::to_string(sum_init) + << "\nSum after iteration: " << std::to_string(sum_after) + << std::endl; + + return 0; +} diff --git a/naaice/CMakeLists.txt b/naaice/CMakeLists.txt new file mode 100644 index 0000000..209dc91 --- /dev/null +++ b/naaice/CMakeLists.txt @@ -0,0 +1,14 @@ +add_executable(naaice BTCS_2D_NAAICE.cpp) +add_executable(NAAICE_dble_vs_float NAAICE_dble_vs_float.cpp) + +get_filename_component(IN_CONC_FILE "init_conc.csv" REALPATH) +get_filename_component(IN_ALPHAX_FILE "alphax.csv" REALPATH) + +configure_file(files.hpp.in files.hpp) +target_include_directories(naaice PRIVATE ${CMAKE_CURRENT_BINARY_DIR}) +target_include_directories(NAAICE_dble_vs_float PRIVATE ${CMAKE_CURRENT_BINARY_DIR}) + +target_compile_definitions(naaice PRIVATE WRITE_THOMAS_CSV) + +target_link_libraries(naaice PUBLIC tug) +target_link_libraries(NAAICE_dble_vs_float PUBLIC tug) diff --git a/naaice/NAAICE_dble_vs_float.cpp b/naaice/NAAICE_dble_vs_float.cpp new file mode 100644 index 0000000..ac3718c --- /dev/null +++ b/naaice/NAAICE_dble_vs_float.cpp @@ -0,0 +1,193 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "files.hpp" + +using namespace tug; + +/** + * Try to parse an input string into a given template type. + */ +template inline T parseString(const std::string &str) { + T result; + std::istringstream iss(str); + + if (!(iss >> result)) { + throw std::invalid_argument("Invalid input for parsing."); + } + + return result; +} + +/** + * Splits a given string into a vector by using a delimiter character. + */ +template +std::vector tokenize(const std::string &input, char delimiter) { + std::vector tokens; + std::istringstream tokenStream(input); + std::string token; + + while (std::getline(tokenStream, token, delimiter)) { + tokens.push_back(parseString(token)); + } + + return tokens; +} + +/** + * Opens a file containing CSV and transform it into row-major 2D STL vector. + */ +template +std::vector> CSVToVector(const char *filename) { + std::ifstream in_file(filename); + if (!in_file.is_open()) { + throw std::runtime_error("Error opening file \'" + std::string(filename) + + "\'."); + } + + std::vector> csv_data; + + std::string line; + while (std::getline(in_file, line)) { + csv_data.push_back(tokenize(line, ',')); + } + + in_file.close(); + + return csv_data; +} + +/** + * Converts a 2D STL vector, where values are stored row-major into a + * column-major Eigen::Matrix. + */ +template +Eigen::MatrixXd rmVecTocmMatrix(const std::vector> &vec, + std::uint32_t exp_rows, + std::uint32_t exp_cols) { + if (exp_rows != vec.size()) { + throw std::runtime_error( + "Mismatch in y dimension while converting to Eigen::Matrix."); + } + + Eigen::MatrixXd out_mat(exp_rows, exp_cols); + + for (std::uint32_t ri = 0; ri < exp_rows; ri++) { + const auto &vec_row = vec[ri]; + if (vec[ri].size() != exp_cols) { + throw std::runtime_error( + "Mismatch in x dimension while converting to Eigen::Matrix."); + } + for (std::uint32_t cj = 0; cj < exp_cols; cj++) { + out_mat(ri, cj) = vec_row[cj]; + } + } + + return out_mat; +} + +template int doWork(int ngrid) { + + constexpr T dt = 10; + + // create a grid + std::string name; + + if constexpr (std::is_same_v) { + name = "DOUBLE"; + } else if constexpr (std::is_same_v) { + name = "FLOAT"; + } else { + name = "unknown"; + } + + 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); + + const T sum_init64 = initConc64.sum(); + + constexpr T alphax_val = 5e-10; + Eigen::MatrixX alphax = + Eigen::MatrixX::Constant(ngrid, ngrid, alphax_val); // row,col,value + constexpr T alphay_val = 1e-10; + Eigen::MatrixX alphay = + Eigen::MatrixX::Constant(ngrid, ngrid, alphay_val); // row,col,value + grid.setAlpha(alphax, alphay); + + // create a boundary with constant values + Boundary bc = Boundary(grid); + bc.setBoundarySideClosed(BC_SIDE_LEFT); + bc.setBoundarySideClosed(BC_SIDE_RIGHT); + bc.setBoundarySideClosed(BC_SIDE_TOP); + bc.setBoundarySideClosed(BC_SIDE_BOTTOM); + + // set up a simulation environment + Simulation Sim = + Simulation(grid, bc); // grid_64,boundary,simulation-approach + + // Sim64.setSolver(THOMAS_ALGORITHM_SOLVER); + + // set the timestep of the simulation + Sim.setTimestep(dt); // timestep + + // set the number of iterations + Sim.setIterations(2); + + // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + Sim.setOutputCSV(CSV_OUTPUT_ON); + + // set output to the console to 'ON' + Sim.setOutputConsole(CONSOLE_OUTPUT_OFF); + + // // **** RUN SIM64 **** + auto begin_t = std::chrono::high_resolution_clock::now(); + + Sim.run(); + + auto end_t = std::chrono::high_resolution_clock::now(); + auto ms_t = + std::chrono::duration_cast(end_t - begin_t); + + const double sum_after64 = grid.getConcentrations().sum(); + + std::cout << "Sum of init field: " << std::setprecision(15) << sum_init64 + << "\nSum after 2 iterations: " << sum_after64 + << "\nMilliseconds: " << ms_t.count() << std::endl + << std::endl; + return 0; +} + +int main(int argc, char *argv[]) { + + int n[] = {101, 201, 501, 1001, 2001}; + + for (int i = 0; i < std::size(n); i++) { + doWork(n[i]); + doWork(n[i]); + doWork(n[i]); + doWork(n[i]); + } + + return 0; +} diff --git a/naaice/README.md b/naaice/README.md new file mode 100644 index 0000000..b4f1f9d --- /dev/null +++ b/naaice/README.md @@ -0,0 +1,70 @@ +This directory contains a concise benchmark designed for validating FPGA +offloading of the Thomas algorithm, primarily employed for solving linear +equation systems structured within a tridiagonal matrix. + + +# Benchmark Setup + +The benchmark defines a domain measuring $1 \text{cm} \times 0.5 \text{cm}$ (easting $\times$ northing), +discretized in a $10 \times 5$ grid. Each grid cell initially +contains a specific concentration. The concentration in the left domain half is set to $6.92023 \times 10^{-7}$, while in the right half to +$2.02396 \times 10^{-8}$, creating an horizontal concentration discontinuity at +the center of the grid. These initial concentrations are read from headerless csv file [init_conc.csv](./init_conc.csv). + +A diffusion time step is simulated with the +heterogeneous 2D-ADI approach detailed in the +[ADI_scheme.pdf](../doc/ADI_scheme.pdf) file. The x component of the +diffusion coefficients, read from headerless csv file [alphax.csv](./alphax.csv) ranges from $\alpha = 10^{-9}$ to $10^{-10}$ (distributed randomly), while the +y-component is held constant at $5 \times 10^{-10}$. Closed +boundary conditions are enforced at all domain boundaries, meaning that concentration cannot enter or exit +the system, or in other terms, that the sum of concentrations over the domain must stay constant. The benchmark simulates a single iteration with a +time step ($\Delta t$) of 360 seconds. + + +# Usage + +To generate new makefiles using the `-DTUG_NAAICE_EXAMPLE=ON` option in CMake, +compile the executable, and run it to generate the benchmark output, follow +these steps: + +1. Navigate to your project's build directory. +2. Run the following CMake command with the `-DTUG_NAAICE_EXAMPLE=ON` option to + generate the makefiles: + + cmake -DTUG_NAAICE_EXAMPLE=ON .. + +3. After CMake configuration is complete, build the `naaice` executable by running `make`: + + make naaice + +4. Once the compilation is successful, navigate to the build directory by `cd + /naaice` + +5. Finally, run the `naaice` executable to generate the benchmark output: + + ./naaice + + +## Output Files + + +### `Thomas_.csv` + +These files contain the values of the tridiagonal coefficient matrix $A$, where: + +- $Aa$ represents the leftmost value, +- $Ab$ represents the middle value, and +- $Ac$ represents the rightmost value of one row of the matrix. + +Additionally, the corresponding values of the right-hand-side vector $b$ are +provided. + +Since the 2D-ADI BTCS scheme processes each row first and then proceeds +column-wise through the grid, each iteration is saved separately in +consecutively numbered files. + + +### `BTCS_5_10_1.csv` + +The result of the simulation, **separated by whitespaces**! + diff --git a/naaice/alphax.csv b/naaice/alphax.csv new file mode 100644 index 0000000..17db4c6 --- /dev/null +++ b/naaice/alphax.csv @@ -0,0 +1,5 @@ +9.35863547772169e-10,4.1475358910393e-10,7.40997499064542e-10,4.63898729439825e-10,5.50232032872736e-10,1.83670640387572e-10,5.84096802584827e-10,4.72976535512134e-10,2.92526249657385e-10,5.8247926668264e-10 +8.08492869278416e-10,3.5943602763582e-10,6.87254608655348e-10,1.52116895909421e-10,5.95404988038354e-10,8.90064929216169e-10,6.13143724552356e-10,1.23507722397335e-10,2.898759260308e-10,9.90528965741396e-10 +3.12359050870873e-10,6.34084914484993e-10,8.28234328492545e-10,1.60584925650619e-10,1.31777092232369e-10,7.34155903453939e-10,9.86383097898215e-10,5.42474469379522e-10,1.23030534153804e-10,2.33146838657558e-10 +7.76231796317734e-10,1.69479642040096e-10,6.74477553134784e-10,2.4903867142275e-10,2.70820381003432e-10,1.67315319390036e-10,7.1631961625535e-10,4.51548034301959e-10,2.41610987577587e-10,4.98075650655665e-10 +7.33895266707987e-10,5.82150935800746e-10,1.49088777275756e-10,5.00345975253731e-10,8.26257051993161e-10,5.28838745504618e-10,9.94136832957156e-10,7.44971914449707e-10,7.53557282453403e-10,7.54089470859617e-10 diff --git a/naaice/files.hpp.in b/naaice/files.hpp.in new file mode 100644 index 0000000..c7dd813 --- /dev/null +++ b/naaice/files.hpp.in @@ -0,0 +1,7 @@ +#ifndef FILES_H_ +#define FILES_H_ + +const char *INPUT_CONC_FILE = "@IN_CONC_FILE@"; +const char *INPUT_ALPHAX_FILE = "@IN_ALPHAX_FILE@"; + +#endif // FILES_H_ diff --git a/naaice/init_conc.csv b/naaice/init_conc.csv new file mode 100644 index 0000000..ad55817 --- /dev/null +++ b/naaice/init_conc.csv @@ -0,0 +1,5 @@ +6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08 +6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08 +6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08 +6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08 +6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08 diff --git a/src/Boundary.cpp b/src/Boundary.cpp deleted file mode 100644 index def9bbd..0000000 --- a/src/Boundary.cpp +++ /dev/null @@ -1,149 +0,0 @@ -#include "TugUtils.hpp" - -#include -#include -#include - -BoundaryElement::BoundaryElement() {} - -BoundaryElement::BoundaryElement(double _value) - : value(_value), type(BC_TYPE_CONSTANT) {} - -void BoundaryElement::setType(BC_TYPE type) { this->type = type; } - -void BoundaryElement::setValue(double value) { - if (value < 0) { - throw_invalid_argument("No negative concentration allowed."); - } - if (type == BC_TYPE_CLOSED) { - throw_invalid_argument( - "No constant boundary concentrations can be set for closed " - "boundaries. Please change type first."); - } - this->value = value; -} - -BC_TYPE BoundaryElement::getType() { return this->type; } - -double BoundaryElement::getValue() { return this->value; } - -Boundary::Boundary(Grid grid) : grid(grid) { - if (grid.getDim() == 1) { - this->boundaries = std::vector>( - 2); // in 1D only left and right boundary - - this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); - this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); - } else if (grid.getDim() == 2) { - this->boundaries = std::vector>(4); - - this->boundaries[BC_SIDE_LEFT] = - std::vector(grid.getRow(), BoundaryElement()); - this->boundaries[BC_SIDE_RIGHT] = - std::vector(grid.getRow(), BoundaryElement()); - this->boundaries[BC_SIDE_TOP] = - std::vector(grid.getCol(), BoundaryElement()); - this->boundaries[BC_SIDE_BOTTOM] = - std::vector(grid.getCol(), BoundaryElement()); - } -} - -void Boundary::setBoundarySideClosed(BC_SIDE side) { - if (grid.getDim() == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw_invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } - - const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; - const int n = is_vertical ? grid.getRow() : grid.getCol(); - - this->boundaries[side] = std::vector(n, BoundaryElement()); -} - -void Boundary::setBoundarySideConstant(BC_SIDE side, double value) { - if (grid.getDim() == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw_invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } - - const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; - const int n = is_vertical ? grid.getRow() : grid.getCol(); - - this->boundaries[side] = - std::vector(n, BoundaryElement(value)); -} - -void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) { - // tests whether the index really points to an element of the boundary side. - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - this->boundaries[side][index].setType(BC_TYPE_CLOSED); -} - -void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, - double value) { - // tests whether the index really points to an element of the boundary side. - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - this->boundaries[side][index].setType(BC_TYPE_CONSTANT); - this->boundaries[side][index].setValue(value); -} - -const std::vector Boundary::getBoundarySide(BC_SIDE side) { - if (grid.getDim() == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw_invalid_argument( - "For the one-dimensional trap, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } - return this->boundaries[side]; -} - -Eigen::VectorXd Boundary::getBoundarySideValues(BC_SIDE side) { - const std::size_t length = boundaries[side].size(); - Eigen::VectorXd values(length); - - for (int i = 0; i < length; i++) { - if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { - values(i) = -1; - continue; - } - values(i) = getBoundaryElementValue(side, i); - } - - return values; -} - -BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) { - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - return this->boundaries[side][index]; -} - -BC_TYPE Boundary::getBoundaryElementType(BC_SIDE side, int index) { - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - return this->boundaries[side][index].getType(); -} - -double Boundary::getBoundaryElementValue(BC_SIDE side, int index) { - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) { - throw_invalid_argument( - "A value can only be output if it is a constant boundary condition."); - } - return this->boundaries[side][index].getValue(); -} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt deleted file mode 100644 index 4bbb32f..0000000 --- a/src/CMakeLists.txt +++ /dev/null @@ -1,9 +0,0 @@ -add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCS.cpp) - -target_link_libraries(tug Eigen3::Eigen) - -if(TUG_USE_OPENMP AND OpenMP_CXX_FOUND) - target_link_libraries(tug OpenMP::OpenMP_CXX) -endif() - -target_include_directories(tug PUBLIC ${PROJECT_SOURCE_DIR}/include) diff --git a/src/Grid.cpp b/src/Grid.cpp deleted file mode 100644 index c6573e4..0000000 --- a/src/Grid.cpp +++ /dev/null @@ -1,167 +0,0 @@ -#include "TugUtils.hpp" - -#include -#include - -constexpr double MAT_INIT_VAL = 0; - -Grid::Grid(int length) : col(length), domainCol(length) { - if (length <= 3) { - throw_invalid_argument( - "Given grid length too small. Must be greater than 3."); - } - - this->dim = 1; - this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 - - this->concentrations = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); - this->alphaX = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); -} - -Grid::Grid(int _row, int _col) - : row(_row), col(_col), domainRow(_row), domainCol(_col) { - if (row <= 3 || col <= 3) { - throw_invalid_argument( - "Given grid dimensions too small. Must each be greater than 3."); - } - - this->dim = 2; - this->deltaRow = double(this->domainRow) / double(this->row); // -> 1 - this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 - - this->concentrations = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); - this->alphaX = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); - this->alphaY = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); -} - -void Grid::setConcentrations(const Eigen::MatrixXd &concentrations) { - if (concentrations.rows() != this->row || - concentrations.cols() != this->col) { - throw_invalid_argument( - "Given matrix of concentrations mismatch with Grid dimensions!"); - } - - this->concentrations = concentrations; -} - -void Grid::setAlpha(const Eigen::MatrixXd &alpha) { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use 2D setter function!"); - } - if (alpha.rows() != 1 || alpha.cols() != this->col) { - throw_invalid_argument( - "Given matrix of alpha coefficients mismatch with Grid dimensions!"); - } - - this->alphaX = alpha; -} - -void Grid::setAlpha(const Eigen::MatrixXd &alphaX, - const Eigen::MatrixXd &alphaY) { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably " - "use 1D setter function!"); - } - if (alphaX.rows() != this->row || alphaX.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients in x-direction " - "mismatch with GRid dimensions!"); - } - if (alphaY.rows() != this->row || alphaY.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients in y-direction " - "mismatch with GRid dimensions!"); - } - - this->alphaX = alphaX; - this->alphaY = alphaY; -} - -const Eigen::MatrixXd &Grid::getAlpha() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use either getAlphaX() or getAlphaY()!"); - } - - return this->alphaX; -} - -const Eigen::MatrixXd &Grid::getAlphaX() { - if (dim != 2) { - throw_invalid_argument( - "Grid is not two dimensional, you should probably use getAlpha()!"); - } - - return this->alphaX; -} - -const Eigen::MatrixXd &Grid::getAlphaY() { - if (dim != 2) { - throw_invalid_argument( - "Grid is not two dimensional, you should probably use getAlpha()!"); - } - - return this->alphaY; -} - -int Grid::getDim() { return dim; } - -int Grid::getLength() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use getRow() or getCol()!"); - } - - return col; -} - -int Grid::getRow() { return row; } - -int Grid::getCol() { return col; } - -void Grid::setDomain(double domainLength) { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probaly " - "use the 2D domain setter!"); - } - if (domainLength <= 0) { - throw_invalid_argument("Given domain length is not positive!"); - } - - this->domainCol = domainLength; - this->deltaCol = double(this->domainCol) / double(this->col); -} - -void Grid::setDomain(double domainRow, double domainCol) { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably " - "use the 1D domain setter!"); - } - if (domainRow <= 0 || domainCol <= 0) { - throw_invalid_argument("Given domain size is not positive!"); - } - - this->domainRow = domainRow; - this->domainCol = domainCol; - this->deltaRow = double(this->domainRow) / double(this->row); - this->deltaCol = double(this->domainCol) / double(this->col); -} - -double Grid::getDelta() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use the 2D delta getters"); - } - - return this->deltaCol; -} - -double Grid::getDeltaCol() { return this->deltaCol; } - -double Grid::getDeltaRow() { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, meaning there is no " - "delta in y-direction!"); - } - - return this->deltaRow; -} diff --git a/src/README.md b/src/README.md deleted file mode 100644 index 82cd165..0000000 --- a/src/README.md +++ /dev/null @@ -1,36 +0,0 @@ -

src-Directory

-This is the src-directory that holds the source code to the implementation of the TUG framework. - -
-src/  
-  ├── CMakeFiles/
-  ├── Boundary.cpp  
-  ├── BTCSv2.cpp
-  ├── CMakeLists.txt  
-  ├── FTCS.cpp
-  ├── Grid.cpp 
-  ├── README.md
-  ├── Simulation.cpp
-  └── TugUtils.hpp 
-
- -**src/** Directory with the source code. - -**CMakeFiles/** Automatically generated directory by CMake. - -**Boundary.cpp** Implementation of Boundary class, that holds the boundary conditions. - -**BTCSv2.cpp** Implementation of BTCS solution to heterogeneous diffusion in 1D and 2D. - -**CMakeLists.txt** CMakeLists for this directory. - -**FTCS.cpp** Implementation of FTCS solution to heterogeneous diffusion in 1D and 2D. - -**Grid.cpp** Implementation of Grid class, that holds all of the concentrations alpha coefficient in x- and y-direction. - -**README.md** This file. - -**Simulation.cpp** Implementation of Simulation class, that holds all of the information for a specific simulation run, as well as the Boundary and Grid objects. - -**TugUtils.hpp** Helper functions for other source files. - diff --git a/src/Schemes.hpp b/src/Schemes.hpp deleted file mode 100644 index d10e0c5..0000000 --- a/src/Schemes.hpp +++ /dev/null @@ -1,28 +0,0 @@ -/** - * @file BTCSv2.cpp - * @brief Implementation of heterogenous BTCS (backward time-centered space) - * solution of diffusion equation in 1D and 2D space. Internally the - * alternating-direction implicit (ADI) method is used. Version 2, because - * Version 1 was an implementation for the homogeneous BTCS solution. - * - */ - -#ifndef SCHEMES_H_ -#define SCHEMES_H_ - -#include "TugUtils.hpp" - -#include -#include - -// entry point; differentiate between 1D and 2D grid -extern void FTCS(Grid &grid, Boundary &bc, double ×tep, int &numThreads); - -// entry point for EigenLU solver; differentiate between 1D and 2D grid -extern void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads); - -// entry point for Thomas algorithm solver; differentiate 1D and 2D grid -extern void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, - int numThreads); - -#endif // SCHEMES_H_ diff --git a/src/Simulation.cpp b/src/Simulation.cpp deleted file mode 100644 index 0a5a961..0000000 --- a/src/Simulation.cpp +++ /dev/null @@ -1,326 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include "Schemes.hpp" -#include "TugUtils.hpp" - -Simulation::Simulation(Grid &_grid, Boundary &_bc, APPROACH _approach) - : grid(_grid), bc(_bc), approach(_approach) {} - -void Simulation::setOutputCSV(CSV_OUTPUT csv_output) { - if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { - throw_invalid_argument("Invalid CSV output option given!"); - } - - this->csv_output = csv_output; -} - -void Simulation::setOutputConsole(CONSOLE_OUTPUT console_output) { - if (console_output < CONSOLE_OUTPUT_OFF && - console_output > CONSOLE_OUTPUT_VERBOSE) { - throw_invalid_argument("Invalid console output option given!"); - } - - this->console_output = console_output; -} - -void Simulation::setTimeMeasure(TIME_MEASURE time_measure) { - if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) { - throw_invalid_argument("Invalid time measure option given!"); - } - - this->time_measure = time_measure; -} - -void Simulation::setTimestep(double timestep) { - if (timestep <= 0) { - throw_invalid_argument("Timestep has to be greater than zero."); - } - - if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { - - const double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - // will be 0 if 1D, else ... - const double deltaRowSquare = grid.getDim() != 1 - ? grid.getDeltaRow() * grid.getDeltaRow() - : deltaColSquare; - const double minDeltaSquare = - (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - - double maxAlpha = std::numeric_limits::quiet_NaN(); - - // determine maximum alpha - if (grid.getDim() == 2) { - - const double maxAlphaX = grid.getAlphaX().maxCoeff(); - const double maxAlphaY = grid.getAlphaY().maxCoeff(); - maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; - - } else if (grid.getDim() == 1) { - maxAlpha = grid.getAlpha().maxCoeff(); - - } else { - throw_invalid_argument("Critical error: Undefined number of dimensions!"); - } - - const std::string dim = std::to_string(grid.getDim()) + "D"; - - // Courant-Friedrichs-Lewy condition - double cfl = minDeltaSquare / (4 * maxAlpha); - - // stability equation from Wikipedia; might be useful if applied cfl does - // not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha * - // ((1/deltaRowSquare) + (1/deltaColSquare))); - - const std::string &approachPrefix = this->approach_names[approach]; - std::cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl - << std::endl; - std::cout << approachPrefix << "_" << dim << " :: required dt=" << timestep - << std::endl; - - if (timestep > cfl) { - - this->innerIterations = (int)ceil(timestep / cfl); - this->timestep = timestep / (double)innerIterations; - - std::cerr << "Warning :: Timestep was adjusted, because of stability " - "conditions. Time duration was approximately preserved by " - "adjusting internal number of iterations." - << std::endl; - std::cout << approachPrefix << "_" << dim << " :: Required " - << this->innerIterations - << " inner iterations with dt=" << this->timestep << std::endl; - - } else { - - this->timestep = timestep; - std::cout << approachPrefix << "_" << dim - << " :: No inner iterations required, dt=" << timestep - << std::endl; - } - - } else { - this->timestep = timestep; - } -} - -double Simulation::getTimestep() { return this->timestep; } - -void Simulation::setIterations(int iterations) { - if (iterations <= 0) { - throw_invalid_argument("Number of iterations must be greater than zero."); - } - this->iterations = iterations; -} - -void Simulation::setSolver(SOLVER solver) { - if (this->approach == FTCS_APPROACH) { - std::cerr - << "Warning: Solver was set, but FTCS approach initialized. Setting " - "the solver " - "is thus without effect." - << std::endl; - } - - this->solver = solver; -} - -void Simulation::setNumberThreads(int numThreads) { - if (numThreads > 0 && numThreads <= omp_get_num_procs()) { - this->numThreads = numThreads; - } else { - int maxThreadNumber = omp_get_num_procs(); - std::string outputMessage = - "Number of threads exceeds the number of processor cores (" + - std::to_string(maxThreadNumber) + ") or is less than 1."; - - throw_invalid_argument(outputMessage); - } -} - -int Simulation::getIterations() { return this->iterations; } - -void Simulation::printConcentrationsConsole() { - std::cout << grid.getConcentrations() << std::endl; - std::cout << std::endl; -} - -std::string Simulation::createCSVfile() { - std::ofstream file; - int appendIdent = 0; - std::string appendIdentString; - - // string approachString = (approach == 0) ? "FTCS" : "BTCS"; - const std::string &approachString = this->approach_names[approach]; - std::string row = std::to_string(grid.getRow()); - std::string col = std::to_string(grid.getCol()); - std::string numIterations = std::to_string(iterations); - - std::string filename = - approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; - - while (std::filesystem::exists(filename)) { - appendIdent += 1; - appendIdentString = std::to_string(appendIdent); - filename = approachString + "_" + row + "_" + col + "_" + numIterations + - "-" + appendIdentString + ".csv"; - } - - file.open(filename); - if (!file) { - exit(1); - } - - // adds lines at the beginning of verbose output csv that represent the - // boundary conditions and their values -1 in case of closed boundary - if (csv_output == CSV_OUTPUT_XTREME) { - Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "", - " "); - file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) - << std::endl; // boundary left - file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) - << std::endl; // boundary right - file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) - << std::endl; // boundary top - file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) - << std::endl; // boundary bottom - file << std::endl << std::endl; - } - - file.close(); - - return filename; -} - -void Simulation::printConcentrationsCSV(const std::string &filename) { - std::ofstream file; - - file.open(filename, std::ios_base::app); - if (!file) { - exit(1); - } - - Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); - file << grid.getConcentrations().format(do_not_align) << std::endl; - file << std::endl << std::endl; - file.close(); -} - -void Simulation::run() { - if (this->timestep == -1) { - throw_invalid_argument("Timestep is not set!"); - } - if (this->iterations == -1) { - throw_invalid_argument("Number of iterations are not set!"); - } - - std::string filename; - if (this->console_output > CONSOLE_OUTPUT_OFF) { - printConcentrationsConsole(); - } - if (this->csv_output > CSV_OUTPUT_OFF) { - filename = createCSVfile(); - } - - auto begin = std::chrono::high_resolution_clock::now(); - - if (approach == FTCS_APPROACH) { // FTCS case - - for (int i = 0; i < iterations * innerIterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } - - FTCS(this->grid, this->bc, this->timestep, this->numThreads); - - // if (i % (iterations * innerIterations / 100) == 0) { - // double percentage = (double)i / ((double)iterations * - // (double)innerIterations) * 100; if ((int)percentage % 10 == 0) { - // cout << "Progress: " << percentage << "%" << endl; - // } - // } - } - - } else if (approach == BTCS_APPROACH) { // BTCS case - - if (solver == EIGEN_LU_SOLVER) { - for (int i = 0; i < iterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } - - BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads); - } - } else if (solver == THOMAS_ALGORITHM_SOLVER) { - for (int i = 0; i < iterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } - - BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); - } - } - - } else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case - - double beta = 0.5; - - // TODO this implementation is very inefficient! - // a separate implementation that sets up a specific tridiagonal matrix for - // Crank-Nicolson would be better - Eigen::MatrixXd concentrations; - Eigen::MatrixXd concentrationsFTCS; - Eigen::MatrixXd concentrationsResult; - for (int i = 0; i < iterations * innerIterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } - - concentrations = grid.getConcentrations(); - FTCS(this->grid, this->bc, this->timestep, this->numThreads); - concentrationsFTCS = grid.getConcentrations(); - grid.setConcentrations(concentrations); - BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); - concentrationsResult = - beta * concentrationsFTCS + (1 - beta) * grid.getConcentrations(); - grid.setConcentrations(concentrationsResult); - } - } - - auto end = std::chrono::high_resolution_clock::now(); - auto milliseconds = - std::chrono::duration_cast(end - begin); - - if (this->console_output > CONSOLE_OUTPUT_OFF) { - printConcentrationsConsole(); - } - if (this->csv_output > CSV_OUTPUT_OFF) { - printConcentrationsCSV(filename); - } - if (this->time_measure > TIME_MEASURE_OFF) { - const std::string &approachString = this->approach_names[approach]; - const std::string dimString = std::to_string(grid.getDim()) + "D"; - std::cout << approachString << dimString << ":: run() finished in " - << milliseconds.count() << "ms" << std::endl; - } -} diff --git a/test/testBoundary.cpp b/test/testBoundary.cpp index 1d2bff2..a7b2b5e 100644 --- a/test/testBoundary.cpp +++ b/test/testBoundary.cpp @@ -6,12 +6,13 @@ #include using namespace std; +using namespace tug; TEST_CASE("BoundaryElement") { SUBCASE("Closed case") { - BoundaryElement boundaryElementClosed = BoundaryElement(); - CHECK_NOTHROW(BoundaryElement()); + BoundaryElement boundaryElementClosed = BoundaryElement(); + CHECK_NOTHROW(BoundaryElement()); CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED); CHECK_EQ(boundaryElementClosed.getValue(), -1); CHECK_THROWS(boundaryElementClosed.setValue(0.2)); @@ -28,11 +29,11 @@ TEST_CASE("BoundaryElement") { } TEST_CASE("Boundary Class") { - Grid grid1D = Grid(10); - Grid grid2D = Grid(10, 12); + Grid grid1D = Grid64(10); + Grid grid2D = Grid64(10, 12); Boundary boundary1D = Boundary(grid1D); Boundary boundary2D = Boundary(grid2D); - vector boundary1DVector(1, BoundaryElement(1.0)); + vector> boundary1DVector(1, BoundaryElement(1.0)); SUBCASE("Boundaries 1D case") { CHECK_NOTHROW(Boundary boundary(grid1D)); diff --git a/test/testFTCS.cpp b/test/testFTCS.cpp index 209ac8a..8e577ad 100644 --- a/test/testFTCS.cpp +++ b/test/testFTCS.cpp @@ -1,4 +1,4 @@ -#include +#include #include #include diff --git a/test/testGrid.cpp b/test/testGrid.cpp index a1f2886..d0179dd 100644 --- a/test/testGrid.cpp +++ b/test/testGrid.cpp @@ -4,26 +4,26 @@ using namespace Eigen; using namespace std; - +using namespace tug; TEST_CASE("1D Grid, too small length") { int l = 2; - CHECK_THROWS(Grid(l)); + CHECK_THROWS(Grid64(l)); } -TEST_CASE("2D Grid, too small side") { +TEST_CASE("2D Grid64, too small side") { int r = 2; int c = 4; - CHECK_THROWS(Grid(r, c)); + CHECK_THROWS(Grid64(r, c)); r = 4; c = 2; - CHECK_THROWS(Grid(r, c)); + CHECK_THROWS(Grid64(r, c)); } -TEST_CASE("1D Grid") { +TEST_CASE("1D Grid64") { int l = 12; - Grid grid(l); + Grid64 grid(l); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 1); @@ -89,9 +89,9 @@ TEST_CASE("1D Grid") { } } -TEST_CASE("2D Grid quadratic") { +TEST_CASE("2D Grid64 quadratic") { int rc = 12; - Grid grid(rc, rc); + Grid64 grid(rc, rc); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 2); @@ -170,10 +170,10 @@ TEST_CASE("2D Grid quadratic") { } } -TEST_CASE("2D Grid non-quadratic") { +TEST_CASE("2D Grid64 non-quadratic") { int r = 12; int c = 15; - Grid grid(r, c); + Grid64 grid(r, c); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 2); diff --git a/test/testSimulation.cpp b/test/testSimulation.cpp index 004799a..ca6888d 100644 --- a/test/testSimulation.cpp +++ b/test/testSimulation.cpp @@ -10,16 +10,16 @@ using namespace Eigen; using namespace std; +using namespace tug; -static Grid setupSimulation(APPROACH approach, double timestep, - int iterations) { +Grid64 setupSimulation(double timestep, int iterations) { int row = 11; int col = 11; int domain_row = 10; int domain_col = 10; // Grid - Grid grid = Grid(row, col); + Grid grid = Grid64(row, col); grid.setDomain(domain_row, domain_col); MatrixXd concentrations = MatrixXd::Constant(row, col, 0); @@ -44,26 +44,29 @@ static Grid setupSimulation(APPROACH approach, double timestep, } grid.setAlpha(alpha, alpha); - // Boundary - Boundary bc = Boundary(grid); - - // Simulation - Simulation sim = Simulation(grid, bc, approach); - // sim.setOutputConsole(CONSOLE_OUTPUT_ON); - sim.setTimestep(timestep); - sim.setIterations(iterations); - sim.run(); - - // RUN return grid; } +constexpr double timestep = 0.001; +constexpr double iterations = 7000; + TEST_CASE("equality to reference matrix with FTCS") { // set string from the header file string test_path = testSimulationCSVDir; MatrixXd reference = CSV2Eigen(test_path); cout << "FTCS Test: " << endl; - Grid grid = setupSimulation(FTCS_APPROACH, 0.001, 7000); + + Grid grid = setupSimulation(timestep, iterations); // Boundary + Boundary bc = Boundary(grid); + + // Simulation + + Simulation sim = Simulation(grid, bc); + // sim.setOutputConsole(CONSOLE_OUTPUT_ON); + sim.setTimestep(timestep); + sim.setIterations(iterations); + sim.run(); + cout << endl; CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true); } @@ -73,24 +76,34 @@ TEST_CASE("equality to reference matrix with BTCS") { string test_path = testSimulationCSVDir; MatrixXd reference = CSV2Eigen(test_path); cout << "BTCS Test: " << endl; - Grid grid = setupSimulation(BTCS_APPROACH, 1, 7); + + Grid grid = setupSimulation(timestep, iterations); // Boundary + Boundary bc = Boundary(grid); + + // Simulation + Simulation sim = Simulation(grid, bc); + // sim.setOutputConsole(CONSOLE_OUTPUT_ON); + sim.setTimestep(timestep); + sim.setIterations(iterations); + sim.run(); + cout << endl; CHECK(checkSimilarityV2(reference, grid.getConcentrations(), 0.01) == true); } TEST_CASE("Initialize environment") { int rc = 12; - Grid grid(rc, rc); + Grid64 grid(rc, rc); Boundary boundary(grid); - CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH)); + CHECK_NOTHROW(Simulation sim(grid, boundary)); } TEST_CASE("Simulation environment") { int rc = 12; - Grid grid(rc, rc); + Grid64 grid(rc, rc); Boundary boundary(grid); - Simulation sim(grid, boundary, FTCS_APPROACH); + Simulation sim(grid, boundary); SUBCASE("default paremeters") { CHECK_EQ(sim.getIterations(), -1); }