From 8456f2212d058eeec118f2d65a034fcb8bc55a3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 17 Oct 2023 12:00:35 +0200 Subject: [PATCH] BREAKING CHANGE: tug as header-only library build: installation of library is now possible --- CMakeLists.txt | 102 ++++++--- cmake/tugConfig.cmake.in | 4 + docs_sphinx/Boundary.rst | 6 +- docs_sphinx/Grid.rst | 2 +- docs_sphinx/Simulation.rst | 8 +- docs_sphinx/developper.rst | 2 - docs_sphinx/images/class_diagram.svg | 2 +- docs_sphinx/theory.rst | 2 +- docs_sphinx/user.rst | 3 +- examples/BTCS_1D_proto_example.cpp | 4 +- examples/BTCS_2D_proto_example.cpp | 3 +- examples/CRNI_2D_proto_example.cpp | 4 +- examples/FTCS_1D_proto_example.cpp | 3 +- examples/FTCS_2D_proto_closed_mdl.cpp | 3 +- examples/FTCS_2D_proto_example.cpp | 3 +- examples/FTCS_2D_proto_example_mdl.cpp | 3 +- examples/profiling_openmp.cpp | 4 +- examples/profiling_speedup.cpp | 4 +- examples/reference-FTCS_2D_closed.cpp | 3 +- include/tug/Boundary.hpp | 18 +- src/BTCS.cpp => include/tug/Core/BTCS.hpp | 19 +- src/FTCS.cpp => include/tug/Core/FTCS.hpp | 13 +- {src => include/tug/Core}/TugUtils.hpp | 0 include/tug/Grid.hpp | 43 ++-- include/tug/Simulation.hpp | 245 ++++++++++++++++++---- naaice/BTCS_2D_NAAICE.cpp | 6 +- src/CMakeLists.txt | 13 -- src/README.md | 36 ---- src/Schemes.hpp | 36 ---- src/Simulation.cpp | 205 ------------------ test/testFTCS.cpp | 2 +- test/testSimulation.cpp | 46 ++-- 32 files changed, 396 insertions(+), 451 deletions(-) create mode 100644 cmake/tugConfig.cmake.in delete mode 100644 docs_sphinx/developper.rst rename src/BTCS.cpp => include/tug/Core/BTCS.hpp (96%) rename src/FTCS.cpp => include/tug/Core/FTCS.hpp (98%) rename {src => include/tug/Core}/TugUtils.hpp (100%) delete mode 100644 src/CMakeLists.txt delete mode 100644 src/README.md delete mode 100644 src/Schemes.hpp delete mode 100644 src/Simulation.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index afed992..78e45c9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,49 +1,93 @@ -#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_HANNESPHILIPP_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}) -option(TUG_NAAICE_EXAMPLE - "Enables NAAICE examples with export of CSV files" - OFF) +include(CMakePackageConfigHelpers) +write_basic_package_version_file( + "tugConfigVersion.cmake" + VERSION ${PROJECT_VERSION} + COMPATIBILITY SameMajorVersion) -add_subdirectory(src) +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) @@ -53,6 +97,6 @@ if(TUG_HANNESPHILIPP_EXAMPLES) add_subdirectory(examples) endif() -if (TUG_NAAICE_EXAMPLE) +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/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 4bfd1b4..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[]) { // ************** @@ -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 f14e0f5..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; @@ -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 d605038..fc3f046 100644 --- a/examples/CRNI_2D_proto_example.cpp +++ b/examples/CRNI_2D_proto_example.cpp @@ -2,6 +2,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { int row = 20; @@ -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 74b6f3e..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[]) { // ************** @@ -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 e17f92b..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[]) { @@ -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 89f821d..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); @@ -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 d77710d..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[]) { @@ -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 c781900..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[]) { @@ -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 c781900..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[]) { @@ -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 fa16063..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; @@ -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 4870ae2..3a7f9b8 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -30,6 +30,8 @@ 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 */ template class BoundaryElement { public: @@ -82,7 +84,7 @@ public: * @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() const { return this->type; } @@ -90,7 +92,7 @@ public: /** * @brief Return the concentration value for the constant boundary condition. * - * @return double Value of the concentration. + * @return Value of the concentration. */ T getValue() const { return this->value; } @@ -102,6 +104,8 @@ private: /** * 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 */ template class Boundary { public: @@ -227,7 +231,7 @@ public: * * @param side Boundary side from which the boundary conditions are to be * returned. - * @return vector> Contains the boundary conditions as + * @return Contains the boundary conditions as * BoundaryElement objects. */ const std::vector> &getBoundarySide(BC_SIDE side) const { @@ -246,7 +250,7 @@ public: specific boundary is closed. * * @param side Boundary side for which the values are to be returned. - * @return VectorX Vector with values as T. + * @return Vector with values as T. */ Eigen::VectorX getBoundarySideValues(BC_SIDE side) const { const std::size_t length = boundaries[side].size(); @@ -271,7 +275,7 @@ 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 + * @return Boundary condition as a BoundaryElement * object. */ BoundaryElement getBoundaryElement(BC_SIDE side, int index) const { @@ -290,7 +294,7 @@ 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) const { if ((boundaries[side].size() < index) || index < 0) { @@ -309,7 +313,7 @@ public: * @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 + * @return Concentration of the corresponding BoundaryElement * object. */ T getBoundaryElementValue(BC_SIDE side, int index) const { diff --git a/src/BTCS.cpp b/include/tug/Core/BTCS.hpp similarity index 96% rename from src/BTCS.cpp rename to include/tug/Core/BTCS.hpp index e166519..cea830a 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 @@ -428,14 +429,6 @@ void BTCS_Thomas(Grid &grid, Boundary &bc, T timestep, int numThreads) { "Error: Only 1- and 2-dimensional grids are defined!"); } } - -template void BTCS_Thomas(Grid &grid, Boundary &bc, - double timestep, int numThreads); -template void BTCS_Thomas(Grid &grid, Boundary &bc, - float timestep, int numThreads); - -template void BTCS_LU(Grid &grid, Boundary &bc, double timestep, - int numThreads); -template void BTCS_LU(Grid &grid, Boundary &bc, float timestep, - int numThreads); } // namespace tug + +#endif // BTCS_H_ diff --git a/src/FTCS.cpp b/include/tug/Core/FTCS.hpp similarity index 98% rename from src/FTCS.cpp rename to include/tug/Core/FTCS.hpp index 914cc41..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 @@ -397,9 +399,6 @@ void FTCS(Grid &grid, Boundary &bc, T timestep, int &numThreads) { "Error: Only 1- and 2-dimensional grids are defined!"); } } - -template void FTCS(Grid &grid, Boundary &bc, double timestep, - int &numThreads); -template void FTCS(Grid &grid, Boundary &bc, float timestep, - int &numThreads); } // 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 0aa9f58..dee004a 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -14,6 +14,12 @@ 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: /** @@ -31,10 +37,11 @@ public: } this->dim = 1; - this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 + this->deltaCol = + static_cast(this->domainCol) / static_cast(this->col); // -> 1 - this->concentrations = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); - this->alphaX = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); + this->concentrations = Eigen::MatrixX::Constant(1, col, MAT_INIT_VAL); + this->alphaX = Eigen::MatrixX::Constant(1, col, MAT_INIT_VAL); } /** @@ -56,8 +63,10 @@ public: } this->dim = 2; - this->deltaRow = double(this->domainRow) / double(this->row); // -> 1 - this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 + 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); @@ -84,7 +93,7 @@ public: /** * @brief Gets the concentrations matrix for a Grid. * - * @return MatrixX An Eigen3 matrix holding the concentrations and having + * @return An Eigen3 matrix holding the concentrations and having * the same dimensions as the grid. */ const Eigen::MatrixX &getConcentrations() { return this->concentrations; } @@ -145,7 +154,7 @@ public: * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one * dimensional. * - * @return MatrixX A matrix with 1 row holding the alpha coefficients. + * @return A matrix with 1 row holding the alpha coefficients. */ const Eigen::MatrixX &getAlpha() const { if (dim != 1) { @@ -161,7 +170,7 @@ public: * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. * Grid must be two dimensional. * - * @return MatrixX A matrix holding the alpha coefficients in x-direction. + * @return A matrix holding the alpha coefficients in x-direction. */ const Eigen::MatrixX &getAlphaX() const { @@ -177,7 +186,7 @@ public: * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. * Grid must be two dimensional. * - * @return MatrixX A matrix holding the alpha coefficients in y-direction. + * @return A matrix holding the alpha coefficients in y-direction. */ const Eigen::MatrixX &getAlphaY() const { @@ -192,14 +201,14 @@ public: /** * @brief Gets the dimensions of the grid. * - * @return int Dimensions, either 1 or 2. + * @return Dimensions, either 1 or 2. */ 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() const { if (dim != 1) { @@ -214,14 +223,14 @@ public: /** * @brief Gets the number of rows of the grid. * - * @return int Number of rows. + * @return Number of rows. */ 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() const { return this->col; } @@ -271,7 +280,7 @@ public: /** * @brief Gets the delta value for 1D-Grid. Grid must be one dimensional. * - * @return double Delta value. + * @return Delta value. */ T getDelta() const { @@ -287,14 +296,14 @@ public: /** * @brief Gets the delta value in x-direction. * - * @return double Delta value in x-direction. + * @return Delta value in x-direction. */ T getDeltaCol() const { return this->deltaCol; } /** * @brief Gets the delta value in y-direction. Must be two dimensional grid. * - * @return double Delta value in y-direction. + * @return Delta value in y-direction. */ T getDeltaRow() const { if (dim != 2) { @@ -318,7 +327,7 @@ private: Eigen::MatrixX alphaX; // Matrix holding alpha coefficients in x-direction Eigen::MatrixX alphaY; // Matrix holding alpha coefficients in y-direction - static constexpr double MAT_INIT_VAL = 0; + static constexpr T MAT_INIT_VAL = 0; }; using Grid64 = Grid; diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index e83d8f4..4c0f5c8 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -11,7 +11,7 @@ #include "Boundary.hpp" #include "Grid.hpp" -#include +#include #include #include #include @@ -19,6 +19,10 @@ #include #include +#include "Core/BTCS.hpp" +#include "Core/FTCS.hpp" +#include "Core/TugUtils.hpp" + #ifdef _OPENMP #include #else @@ -28,13 +32,13 @@ 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 */ }; /** @@ -42,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 */ }; /** @@ -52,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 */ }; /** @@ -64,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 */ }; /** @@ -74,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 */ }; /** @@ -83,9 +87,14 @@ 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 { +template +class Simulation { public: /** * @brief Set up a simulation environment. The timestep and number of @@ -99,8 +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) - : grid(_grid), bc(_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. @@ -169,7 +177,65 @@ public: * * @param timestep Valid timestep greater than zero. */ - void setTimestep(T 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. @@ -192,25 +258,6 @@ public: this->iterations = 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) { - if (this->approach == FTCS_APPROACH) { - std::cerr - << "Warning: Solver was set, but FTCS approach initialized. Setting " - "the solver " - "is thus without effect." - << std::endl; - } - - this->solver = solver; - } - /** * @brief Set the number of desired openMP Threads. * @@ -327,7 +374,117 @@ public: * @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: T timestep{-1}; @@ -340,8 +497,6 @@ private: Grid &grid; Boundary &bc; - APPROACH approach; - SOLVER solver{THOMAS_ALGORITHM_SOLVER}; const std::vector approach_names = {"FTCS", "BTCS", "CRNI"}; }; diff --git a/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp index c356aea..9ed7ccb 100644 --- a/naaice/BTCS_2D_NAAICE.cpp +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -121,7 +121,8 @@ int main(int argc, char *argv[]) { 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 + Eigen::MatrixXd alphay = + Eigen::MatrixXd::Constant(row, col, alphay_val); // row,col,value grid.setAlpha(alphax, alphay); // // ****************** @@ -140,8 +141,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(360); // timestep diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt deleted file mode 100644 index 91d141d..0000000 --- a/src/CMakeLists.txt +++ /dev/null @@ -1,13 +0,0 @@ -add_library(tug Simulation.cpp FTCS.cpp BTCS.cpp) - -IF(TUG_NAAICE_EXAMPLE) - target_compile_definitions(tug PRIVATE WRITE_THOMAS_CSV) -endif() - -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/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 09e6d9e..0000000 --- a/src/Schemes.hpp +++ /dev/null @@ -1,36 +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 - -namespace tug { - -// entry point; differentiate between 1D and 2D grid -template -extern void FTCS(tug::Grid &grid, tug::Boundary &bc, T timestep, - int &numThreads); - -// entry point for EigenLU solver; differentiate between 1D and 2D grid -template -extern void BTCS_LU(tug::Grid &grid, tug::Boundary &bc, T timestep, - int numThreads); - -// entry point for Thomas algorithm solver; differentiate 1D and 2D grid -template -extern void BTCS_Thomas(tug::Grid &grid, tug::Boundary &bc, T timestep, - int numThreads); -} // namespace tug - -#endif // SCHEMES_H_ diff --git a/src/Simulation.cpp b/src/Simulation.cpp deleted file mode 100644 index ca623d4..0000000 --- a/src/Simulation.cpp +++ /dev/null @@ -1,205 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include "Schemes.hpp" -#include "TugUtils.hpp" - -namespace tug { - -template void Simulation::setTimestep(T timestep) { - if (timestep <= 0) { - throw_invalid_argument("Timestep has to be greater than zero."); - } - - if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { - - const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - // will be 0 if 1D, else ... - const T deltaRowSquare = grid.getDim() != 1 - ? grid.getDeltaRow() * grid.getDeltaRow() - : deltaColSquare; - const T minDeltaSquare = - (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - - T maxAlpha = std::numeric_limits::quiet_NaN(); - - // determine maximum alpha - if (grid.getDim() == 2) { - - const T maxAlphaX = grid.getAlphaX().maxCoeff(); - const T 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 - T cfl = minDeltaSquare / (4 * maxAlpha); - - // stability equation from Wikipedia; might be useful if applied cfl does - // not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha * - // ((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; - } -} - -template 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 - - 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; - } -} - -template void Simulation::setTimestep(double timestep); -template void Simulation::setTimestep(float timestep); - -template void Simulation::run(); -template void Simulation::run(); -} // namespace tug 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/testSimulation.cpp b/test/testSimulation.cpp index a7600c5..ca6888d 100644 --- a/test/testSimulation.cpp +++ b/test/testSimulation.cpp @@ -12,8 +12,7 @@ using namespace Eigen; using namespace std; using namespace tug; -static Grid64 setupSimulation(APPROACH approach, double timestep, - int iterations) { +Grid64 setupSimulation(double timestep, int iterations) { int row = 11; int col = 11; int domain_row = 10; @@ -45,26 +44,29 @@ static Grid64 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); } @@ -74,7 +76,17 @@ 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); } @@ -84,14 +96,14 @@ TEST_CASE("Initialize environment") { 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; 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); }