Merge branch 'naaice' into 'main'

BREAKING CHANGE: tug as header-only library

See merge request naaice/tug!18

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")
+include(GNUInstallDirs)
-option(TUG_USE_OPENMP "Compile with OpenMP support" ON)
 
-set(CMAKE_CXX_FLAGS_GENERICOPT "-O3 -march=native" CACHE STRING
+# find_package(easy_profiler) option(EASY_OPTION_LOG "Verbose easy_profiler" 1)
- "Flags used by the C++ compiler during opt builds."
-  FORCE)
 
-set(CMAKE_BUILD_TYPE "${CMAKE_BUILD_TYPE}" CACHE STRING
+# SET(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -O2 -mfma")
-  "Choose the type of build, options are: None Debug Release RelWithDebInfo MinSizeRel GenericOpt."
+option(TUG_USE_OPENMP "Compile tug with OpenMP support" ON)
-  FORCE)
 
-option(TUG_USE_UNSAFE_MATH_OPT
+set(CMAKE_CXX_FLAGS_GENERICOPT
-  "Use compiler options to break IEEE compliances by
+    "-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."
+    points." OFF)
-  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 $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
+                $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
+
+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
+if(TUG_USE_UNSAFE_MATH_OPT)
-  "Run tests after succesfull compilation"
+  target_compile_options(tug INTERFACE -ffast-math)
-  OFF)
+endif()
 
-option(TUG_ENABLE_EXAMPLES
+install(
-  "Compile example applications"
+  TARGETS tug
-  OFF)
+  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()

cmake/tugConfig.cmake.in

@@ -0,0 +1,4 @@
+@PACKAGE_INIT@
+
+include("${CMAKE_CURRENT_LIST_DIR}/@PROJECT_NAME@Targets.cmake")
+check_required_components("@PROJECT_NAME@")

docs_sphinx/Boundary.rst

@@ -1,4 +1,8 @@
 Boundary
 ========
 
-.. doxygenclass:: Boundary
+.. doxygenenum:: tug::BC_TYPE
+.. doxygenenum:: tug::BC_SIDE
+
+.. doxygenclass:: tug::Boundary
+.. doxygenclass:: tug::BoundaryElement

docs_sphinx/Grid.rst

@@ -1,4 +1,4 @@
 Grid
 ====
 
-.. doxygenclass:: Grid
+.. doxygenclass:: tug::Grid

docs_sphinx/Simulation.rst

@@ -1,4 +1,10 @@
 Simulation
 ==========
 
-.. doxygenclass:: Simulation
+.. doxygenenum:: tug::APPROACH
+.. doxygenenum:: tug::SOLVER
+.. doxygenenum:: tug::CSV_OUTPUT
+.. doxygenenum:: tug::CONSOLE_OUTPUT
+.. doxygenenum:: tug::TIME_MEASURE
+
+.. doxygenclass:: tug::Simulation

docs_sphinx/developper.rst

@@ -1,2 +0,0 @@
-Developper API
-==============

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 
+For example, the initalization of a grid with 20 by 20 cells using double values
-also 20 by 20 length units can be done as follows. The setting of the domain is optional here and is set to the 
+and a domain size (physical extent of the grid) of also 20 by 20 length units
-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 
+can be done as follows. The setting of the domain is optional here and is set to
-initial value of 0 and only in the upper left corner (0,0) the starting concentration is set to the value 20.
+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<double> grid(row,col);
     grid.setDomain(row, col);
     MatrixXd concentrations = MatrixXd::Constant(row,col,0);
+    // or MatrixX<double> concentrations = MatrixX<double>::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
+In the last block, a simulation class is created and the objects of the grid and
-parameter specifies the output of the simulated results in a CSV file. In the present case, the result of each iteration step is written
+the boundary conditions are passed. The solution method is also specified
-one below the other into the corresponding CSV file.
+(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<double, FTCS_APPROACH> 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
--------------------------------------------------------
+-------------------------------------------------------

docs_sphinx/theory.rst

@@ -12,4 +12,4 @@ Numerical Solver
 **Backward Time-Centered Space (BTCS) Method**
 
 
-**Forward Time-Centered Space (BTCS) Method**
+**Forward Time-Centered Space (FTCS) Method**

docs_sphinx/user.rst

@@ -3,7 +3,6 @@ User API
 
 .. toctree::    
   
-    :maxdepth: 2
     Grid
     Boundary
-    Simulation
+    Simulation

examples/BTCS_1D_proto_example.cpp

@@ -2,6 +2,7 @@
 #include <tug/Simulation.hpp>
 
 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 simulation = Simulation(grid, bc); // grid,boundary
-      Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach
 
   // set the timestep of the simulation
   simulation.setTimestep(0.1); // timestep

examples/BTCS_2D_proto_example.cpp

@@ -2,6 +2,7 @@
 #include <tug/Simulation.hpp>
 
 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

examples/CRNI_2D_proto_example.cpp

@@ -2,11 +2,12 @@
 #include <tug/Simulation.hpp>
 
 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<double, tug::CRANK_NICOLSON_APPROACH>(grid, bc);
   simulation.setTimestep(0.1);
   simulation.setIterations(50);
   simulation.setOutputCSV(CSV_OUTPUT_XTREME);

examples/FTCS_1D_proto_example.cpp

@@ -3,6 +3,7 @@
 #include <tug/Simulation.hpp>
 
 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<double, tug::FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach
 
   // (optional) set the timestep of the simulation
   simulation.setTimestep(0.1); // timestep

examples/FTCS_2D_proto_closed_mdl.cpp

@@ -12,6 +12,7 @@
 #include <tug/Simulation.hpp>
 
 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<double, FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach
 
   // set the timestep of the simulation
   simulation.setTimestep(10000); // timestep

examples/FTCS_2D_proto_example.cpp

@@ -10,6 +10,7 @@
 #include <tug/Simulation.hpp>
 
 using namespace Eigen;
+using namespace tug;
 // #include <easy/profiler.h>
 // #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<double, tug::FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach
 
   // set the timestep of the simulation
   simulation.setTimestep(0.1); // timestep

examples/FTCS_2D_proto_example_mdl.cpp

@@ -10,6 +10,7 @@
 #include <tug/Simulation.hpp>
 
 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<double, tug::FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach
 
   // (optional) set the timestep of the simulation
   simulation.setTimestep(1000); // timestep

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);
+          Simulation sim = Simulation(grid, bc);
-          sim.setSolver(THOMAS_ALGORITHM_SOLVER);
 
           if (argc == 2) {
             int numThreads = atoi(argv[1]);

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);
+          Simulation sim = Simulation(grid, bc);
-          sim.setSolver(THOMAS_ALGORITHM_SOLVER);
 
           if (argc == 2) {
             int numThreads = atoi(argv[1]);

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<double, tug::FTCS_APPROACH>(grid, bc);
   sim.setTimestep(0.001);
   sim.setIterations(10000);
   sim.setOutputCSV(CSV_OUTPUT_OFF);

include/tug/Boundary.hpp

@@ -9,6 +9,10 @@
 
 #include "Grid.hpp"
 #include <cstddef>
+#include <cstdint>
+#include <exception>
+
+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 T> 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 T> 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<T> &grid)
+      : dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) {
+    if (this->dim == 1) {
+      this->boundaries = std::vector<std::vector<BoundaryElement<T>>>(
+          2); // in 1D only left and right boundary
+
+      this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement<T>());
+      this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement<T>());
+    } else if (this->dim == 2) {
+      this->boundaries = std::vector<std::vector<BoundaryElement<T>>>(4);
+
+      this->boundaries[BC_SIDE_LEFT] =
+          std::vector<BoundaryElement<T>>(this->rows, BoundaryElement<T>());
+      this->boundaries[BC_SIDE_RIGHT] =
+          std::vector<BoundaryElement<T>>(this->rows, BoundaryElement<T>());
+      this->boundaries[BC_SIDE_TOP] =
+          std::vector<BoundaryElement<T>>(this->cols, BoundaryElement<T>());
+      this->boundaries[BC_SIDE_BOTTOM] =
+          std::vector<BoundaryElement<T>>(this->cols, BoundaryElement<T>());
+    }
+  }
 
   /**
    * @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<BoundaryElement<T>>(n, BoundaryElement<T>());
+  }
 
   /**
    * @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<BoundaryElement<T>>(n, BoundaryElement<T>(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<BoundaryElement> Contains the boundary conditions as
+   * @return Contains the boundary conditions as
-   * BoundaryElement objects.
+   * BoundaryElement<T> objects.
    */
-  const std::vector<BoundaryElement> getBoundarySide(BC_SIDE side);
+  const std::vector<BoundaryElement<T>> &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<T> getBoundarySideValues(BC_SIDE side) const {
+    const std::size_t length = boundaries[side].size();
+    Eigen::VectorX<T> 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<T>
+   * object.
    */
-  BoundaryElement getBoundaryElement(BC_SIDE side, int index);
+  BoundaryElement<T> 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<T> 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<T>
+   * 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<BoundaryElement>>
+  std::vector<std::vector<BoundaryElement<T>>>
       boundaries; // Vector with Boundary Element information
 };
-
+} // namespace tug
 #endif // BOUNDARY_H_

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 <Eigen/src/Core/util/Meta.h>
 #include <cstddef>
 #include <tug/Boundary.hpp>
 #include <tug/Grid.hpp>
@@ -22,6 +23,8 @@
 #define omp_get_thread_num() 0
 #endif
 
+namespace tug {
+
 // calculates coefficient for boundary in constant case
 template <class T>
 constexpr std::pair<T, T> calcBoundaryCoeffConstant(T alpha_center,
@@ -45,13 +48,15 @@ constexpr std::pair<T, T> calcBoundaryCoeffClosed(T alpha_center, T alpha_side,
 }
 
 // creates coefficient matrix for next time step from alphas in x-direction
-static Eigen::SparseMatrix<double>
+template <class T>
-createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector<BoundaryElement> &bcLeft,
+static Eigen::SparseMatrix<T>
-                  std::vector<BoundaryElement> &bcRight, int numCols,
+createCoeffMatrix(const Eigen::MatrixX<T> &alpha,
-                  int rowIndex, double sx) {
+                  const std::vector<BoundaryElement<T>> &bcLeft,
+                  const std::vector<BoundaryElement<T>> &bcRight, int numCols,
+                  int rowIndex, T sx) {
 
   // square matrix of column^2 dimension for the coefficients
-  Eigen::SparseMatrix<double> cm(numCols, numCols);
+  Eigen::SparseMatrix<T> 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(
+template <class T>
-    Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alphaX,
+static Eigen::VectorX<T>
-    Eigen::MatrixXd &alphaY, std::vector<BoundaryElement> &bcLeft,
+createSolutionVector(const Eigen::MatrixX<T> &concentrations,
-    std::vector<BoundaryElement> &bcRight, std::vector<BoundaryElement> &bcTop,
+                     const Eigen::MatrixX<T> &alphaX,
-    std::vector<BoundaryElement> &bcBottom, int length, int rowIndex, double sx,
+                     const Eigen::MatrixX<T> &alphaY,
-    double sy) {
+                     const std::vector<BoundaryElement<T>> &bcLeft,
+                     const std::vector<BoundaryElement<T>> &bcRight,
+                     const std::vector<BoundaryElement<T>> &bcTop,
+                     const std::vector<BoundaryElement<T>> &bcBottom,
+                     int length, int rowIndex, T sx, T sy) {
 
-  Eigen::VectorXd sv(length);
+  Eigen::VectorX<T> 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<double> &A,
+template <class T>
-                                        Eigen::VectorXd &b) {
+static Eigen::VectorX<T> EigenLUAlgorithm(Eigen::SparseMatrix<T> &A,
+                                          Eigen::VectorX<T> &b) {
 
-  Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
+  Eigen::SparseLU<Eigen::SparseMatrix<T>> solver;
   solver.analyzePattern(A);
   solver.factorize(A);
 
@@ -248,14 +258,15 @@ static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix<double> &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<double> &A,
+template <class T>
-                                       Eigen::VectorXd &b) {
+static Eigen::VectorX<T> ThomasAlgorithm(Eigen::SparseMatrix<T> &A,
+                                         Eigen::VectorX<T> &b) {
   Eigen::Index n = b.size();
 
-  Eigen::VectorXd a_diag(n);
+  Eigen::VectorX<T> a_diag(n);
-  Eigen::VectorXd b_diag(n);
+  Eigen::VectorX<T> b_diag(n);
-  Eigen::VectorXd c_diag(n);
+  Eigen::VectorX<T> c_diag(n);
-  Eigen::VectorXd x_vec = b;
+  Eigen::VectorX<T> x_vec = b;
 
   // Fill diagonals vectors
   b_diag[0] = A.coeff(0, 0);
@@ -269,6 +280,29 @@ static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix<double> &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 <fstream>
+#include <string>
+  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<double> &A,
 }
 
 // BTCS solution for 1D grid
-static void
+template <class T>
-BTCS_1D(Grid &grid, Boundary &bc, double timestep,
+static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
-        Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix<double> &A,
+                    Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
-                                      Eigen::VectorXd &b)) {
+                                                    Eigen::VectorX<T> &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<T> concentrations_t1(length);
 
-  Eigen::SparseMatrix<double> A;
+  Eigen::SparseMatrix<T> A;
-  Eigen::VectorXd b(length);
+  Eigen::VectorX<T> b(length);
 
-  Eigen::MatrixXd alpha = grid.getAlpha();
+  const auto &alpha = grid.getAlpha();
-  std::vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
-  std::vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
 
-  Eigen::MatrixXd concentrations = grid.getConcentrations();
+  const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
+  const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
+
+  Eigen::MatrixX<T> 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
+template <class T>
-BTCS_2D(Grid &grid, Boundary &bc, double timestep,
+static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
-        Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix<double> &A,
+                    Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
-                                      Eigen::VectorXd &b),
+                                                    Eigen::VectorX<T> &b),
-        int numThreads) {
+                    int numThreads) {
   int rowMax = grid.getRow();
   int colMax = grid.getCol();
-  double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
+  T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
-  double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
+  T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
 
-  Eigen::MatrixXd concentrations_t1 =
+  Eigen::MatrixX<T> concentrations_t1 =
-      Eigen::MatrixXd::Constant(rowMax, colMax, 0);
+      Eigen::MatrixX<T>::Constant(rowMax, colMax, 0);
-  Eigen::VectorXd row_t1(colMax);
+  Eigen::VectorX<T> row_t1(colMax);
 
-  Eigen::SparseMatrix<double> A;
+  Eigen::SparseMatrix<T> A;
-  Eigen::VectorXd b;
+  Eigen::VectorX<T> b;
 
-  Eigen::MatrixXd alphaX = grid.getAlphaX();
+  auto alphaX = grid.getAlphaX();
-  Eigen::MatrixXd alphaY = grid.getAlphaY();
+  auto alphaY = grid.getAlphaY();
-  std::vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
-  std::vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
-  std::vector<BoundaryElement> bcTop = bc.getBoundarySide(BC_SIDE_TOP);
-  std::vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);
 
-  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<T> 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<Eigen::SparseMatrix<double>> 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 <class T>
+void BTCS_LU(Grid<T> &grid, Boundary<T> &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 <class T>
+void BTCS_Thomas(Grid<T> &grid, Boundary<T> &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_

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 <cstddef>
@@ -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 <class T>
+static inline T calcHorizontalChange(Grid<T> &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 <class T>
+static inline T calcVerticalChange(Grid<T> &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,
+template <class T>
-                                                              Boundary &bc,
+static inline T calcHorizontalChangeLeftBoundaryConstant(Grid<T> &grid,
-                                                              int &row,
+                                                         Boundary<T> &bc,
-                                                              int &col) {
+                                                         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
+template <class T>
-calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) {
+static inline T calcHorizontalChangeLeftBoundaryClosed(Grid<T> &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,
+template <class T>
-                                                      int &row, int &col) {
+static inline T calcHorizontalChangeLeftBoundary(Grid<T> &grid, Boundary<T> &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,
+template <class T>
-                                                               Boundary &bc,
+static inline T calcHorizontalChangeRightBoundaryConstant(Grid<T> &grid,
-                                                               int &row,
+                                                          Boundary<T> &bc,
-                                                               int &col) {
+                                                          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
+template <class T>
-calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) {
+static inline T calcHorizontalChangeRightBoundaryClosed(Grid<T> &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,
+template <class T>
-                                                       int &row, int &col) {
+static inline T calcHorizontalChangeRightBoundary(Grid<T> &grid,
+                                                  Boundary<T> &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,
+template <class T>
-                                                           Boundary &bc,
+static inline T calcVerticalChangeTopBoundaryConstant(Grid<T> &grid,
-                                                           int &row, int &col) {
+                                                      Boundary<T> &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,
+template <class T>
-                                                         int &col) {
+static inline T calcVerticalChangeTopBoundaryClosed(Grid<T> &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,
+template <class T>
-                                                   int &row, int &col) {
+static inline T calcVerticalChangeTopBoundary(Grid<T> &grid, Boundary<T> &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,
+template <class T>
-                                                              Boundary &bc,
+static inline T calcVerticalChangeBottomBoundaryConstant(Grid<T> &grid,
-                                                              int &row,
+                                                         Boundary<T> &bc,
-                                                              int &col) {
+                                                         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
+template <class T>
-calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) {
+static inline T calcVerticalChangeBottomBoundaryClosed(Grid<T> &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,
+template <class T>
-                                                      int &row, int &col) {
+static inline T calcVerticalChangeBottomBoundary(Grid<T> &grid, Boundary<T> &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 &timestep) {
+template <class T>
+static void FTCS_1D(Grid<T> &grid, Boundary<T> &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<T> concentrations_t1 =
+      Eigen::MatrixX<T>::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 &timestep) {
 }
 
 // FTCS solution for 2D grid
-static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep,
+template <class T>
+static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
                     int numThreads) {
   int rowMax = grid.getRow();
   int colMax = grid.getCol();
-  double deltaRow = grid.getDeltaRow();
+  T deltaRow = grid.getDeltaRow();
-  double deltaCol = grid.getDeltaCol();
+  T deltaCol = grid.getDeltaCol();
 
   // matrix for concentrations at time t+1
-  Eigen::MatrixXd concentrations_t1 =
+  Eigen::MatrixX<T> concentrations_t1 =
-      Eigen::MatrixXd::Constant(rowMax, colMax, 0);
+      Eigen::MatrixX<T>::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 &timestep,
 }
 
 // entry point; differentiate between 1D and 2D grid
-void FTCS(Grid &grid, Boundary &bc, double &timestep, int &numThreads) {
+template <class T>
+void FTCS(Grid<T> &grid, Boundary<T> &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 &timestep, int &numThreads) {
         "Error: Only 1- and 2-dimensional grids are defined!");
   }
 }
+} // namespace tug
+
+#endif // FTCS_H_

include/tug/Grid.hpp

@@ -10,8 +10,17 @@
 
 #include <Eigen/Core>
 #include <Eigen/Sparse>
+#include <stdexcept>
 
-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 T> 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<T>(this->domainCol) / static_cast<T>(this->col); // -> 1
+
+    this->concentrations = Eigen::MatrixX<T>::Constant(1, col, MAT_INIT_VAL);
+    this->alphaX = Eigen::MatrixX<T>::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<T>(this->domainRow) / static_cast<T>(this->row); // -> 1
+    this->deltaCol =
+        static_cast<T>(this->domainCol) / static_cast<T>(this->col); // -> 1
+
+    this->concentrations = Eigen::MatrixX<T>::Constant(row, col, MAT_INIT_VAL);
+    this->alphaX = Eigen::MatrixX<T>::Constant(row, col, MAT_INIT_VAL);
+    this->alphaY = Eigen::MatrixX<T>::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
+   * @param concentrations An Eigen3 MatrixX<T> holding the concentrations.
-   * must have correct dimensions as defined in row and col. (Or length, in 1D
+   * Matrix must have correct dimensions as defined in row and col. (Or length,
-   * case).
+   * in 1D case).
    */
-  void setConcentrations(const Eigen::MatrixXd &concentrations);
+  void setConcentrations(const Eigen::MatrixX<T> &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
+   * @return An Eigen3 matrix holding the concentrations and having
-   *                  same dimensions as the grid.
+   * the same dimensions as the grid.
    */
-  const Eigen::MatrixXd &getConcentrations() { return this->concentrations; }
+  const Eigen::MatrixX<T> &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.
+   * @param alpha An Eigen3 MatrixX<T> with 1 row holding the alpha
-   * Matrix columns must have same size as length of grid.
+   * coefficients. Matrix columns must have same size as length of grid.
    */
-  void setAlpha(const Eigen::MatrixXd &alpha);
+  void setAlpha(const Eigen::MatrixX<T> &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<T> 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<T> 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<T> &alphaX,
+                const Eigen::MatrixX<T> &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<T> &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<T> &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<T> &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
+   * @param domainRow A double value of the domain size in y-direction. Must
-   * positive.
+   * be positive.
-   * @param domainCol A double value of the domain size in x-direction. Must be
+   * @param domainCol A double value of the domain size in x-direction. Must
-   * positive.
+   * 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 col;                          // number of grid columns
-  int row{1};                     // number of grid rows
+  int row{1};                       // number of grid rows
-  int dim;                        // 1D or 2D
+  int dim;                          // 1D or 2D
-  double domainCol;               // number of domain columns
+  T domainCol;                      // number of domain columns
-  double domainRow{0};            // number of domain rows
+  T domainRow{0};                   // number of domain rows
-  double deltaCol;                // delta in x-direction (between columns)
+  T deltaCol;                       // delta in x-direction (between columns)
-  double deltaRow{0};             // delta in y-direction (between rows)
+  T deltaRow{0};                    // delta in y-direction (between rows)
-  Eigen::MatrixXd concentrations; // Matrix holding grid concentrations
+  Eigen::MatrixX<T> concentrations; // Matrix holding grid concentrations
-  Eigen::MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction
+  Eigen::MatrixX<T> alphaX; // Matrix holding alpha coefficients in x-direction
-  Eigen::MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction
+  Eigen::MatrixX<T> alphaY; // Matrix holding alpha coefficients in y-direction
+
+  static constexpr T MAT_INIT_VAL = 0;
 };
 
+using Grid64 = Grid<double>;
+using Grid32 = Grid<float>;
+} // namespace tug
 #endif // GRID_H_

include/tug/Simulation.hpp

@@ -11,23 +11,34 @@
 
 #include "Boundary.hpp"
 #include "Grid.hpp"
+#include <algorithm>
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+#include <stdexcept>
 #include <string>
 #include <vector>
 
+#include "Core/BTCS.hpp"
+#include "Core/FTCS.hpp"
+#include "Core/TugUtils.hpp"
+
 #ifdef _OPENMP
 #include <omp.h>
 #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
+  FTCS_APPROACH,          /*!< Forward Time-Centered Space */
-  BTCS_APPROACH, // Backward Time-Centered Space solved with EigenLU solver
+  BTCS_APPROACH,          /*!< Backward Time-Centered Space */
-  CRANK_NICOLSON_APPROACH
+  CRANK_NICOLSON_APPROACH /*!< Crank-Nicolson method */
 };
 
 /**
@@ -35,9 +46,9 @@ enum APPROACH {
  *
  */
 enum SOLVER {
-  EIGEN_LU_SOLVER,        // EigenLU solver
+  EIGEN_LU_SOLVER,        /*!<  EigenLU solver */
-  THOMAS_ALGORITHM_SOLVER // Thomas Algorithm solver; more efficient for
+  THOMAS_ALGORITHM_SOLVER /*!< Thomas Algorithm solver; more efficient for
-                          // tridiagonal matrices
+                             tridiagonal matrices */
 };
 
 /**
@@ -45,11 +56,11 @@ enum SOLVER {
  *
  */
 enum CSV_OUTPUT {
-  CSV_OUTPUT_OFF,     // do not produce csv output
+  CSV_OUTPUT_OFF,     /*!< do not produce csv output */
-  CSV_OUTPUT_ON,      // produce csv output with last concentration matrix
+  CSV_OUTPUT_ON,      /*!< produce csv output with last concentration matrix */
-  CSV_OUTPUT_VERBOSE, // produce csv output with all concentration matrices
+  CSV_OUTPUT_VERBOSE, /*!< produce csv output with all concentration matrices */
-  CSV_OUTPUT_XTREME   // csv output like VERBOSE but additional boundary
+  CSV_OUTPUT_XTREME   /*!< csv output like VERBOSE but additional boundary
-                      // conditions at beginning
+                         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_OFF, /*!< do not print any output to console */
-  CONSOLE_OUTPUT_ON,     // print before and after concentrations to console
+  CONSOLE_OUTPUT_ON,  /*!< print before and after concentrations to console */
-  CONSOLE_OUTPUT_VERBOSE // print all concentration matrices 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_OFF, /*!< do not print any time measures */
-  TIME_MEASURE_ON   // print time measure after last iteration
+  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 T, APPROACH approach = BTCS_APPROACH,
+          SOLVER solver = THOMAS_ALGORITHM_SOLVER>
 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<T> &_grid, Boundary<T> &_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);
+  void setIterations(int iterations) {
-
+    if (iterations <= 0) {
-  /**
+      throw std::invalid_argument(
-   * @brief Set the desired linear equation solver to be used for BTCS approach.
+          "Number of iterations must be greater than zero.");
-   * Without effect in case of FTCS approach.
+    }
-   *
+    this->iterations = iterations;
-   * @param solver Solver to be used. Default is Thomas Algorithm as it is more
+  }
-   * efficient for tridiagonal Matrices.
-   */
-  void setSolver(SOLVER solver);
 
   /**
    * @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<T> concentrations;
+      Eigen::MatrixX<T> concentrationsFTCS;
+      Eigen::MatrixX<T> 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<std::chrono::milliseconds>(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;
+  Grid<T> &grid;
-  Boundary &bc;
+  Boundary<T> &bc;
-  APPROACH approach;
-  SOLVER solver{THOMAS_ALGORITHM_SOLVER};
 
   const std::vector<std::string> approach_names = {"FTCS", "BTCS", "CRNI"};
 };
-
+} // namespace tug
 #endif // SIMULATION_H_

naaice/BTCS_2D_NAAICE.cpp

@@ -0,0 +1,169 @@
+#include <Eigen/Eigen>
+#include <cstdint>
+#include <fstream>
+#include <iostream>
+#include <ostream>
+#include <stdexcept>
+#include <string>
+#include <string_view>
+#include <tug/Simulation.hpp>
+#include <vector>
+
+#include "files.hpp"
+
+using namespace tug;
+
+/**
+ * Try to parse an input string into a given template type.
+ */
+template <typename T> 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 <typename T>
+std::vector<T> tokenize(const std::string &input, char delimiter) {
+  std::vector<T> tokens;
+  std::istringstream tokenStream(input);
+  std::string token;
+
+  while (std::getline(tokenStream, token, delimiter)) {
+    tokens.push_back(parseString<T>(token));
+  }
+
+  return tokens;
+}
+
+/**
+ * Opens a file containing CSV and transform it into row-major 2D STL vector.
+ */
+template <typename T>
+std::vector<std::vector<T>> 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<std::vector<T>> csv_data;
+
+  std::string line;
+  while (std::getline(in_file, line)) {
+    csv_data.push_back(tokenize<T>(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 <typename T>
+Eigen::MatrixXd rmVecTocmMatrix(const std::vector<std::vector<T>> &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<double>(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<double>(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;
+}

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)

naaice/NAAICE_dble_vs_float.cpp

@@ -0,0 +1,193 @@
+#include <Eigen/Eigen>
+#include <chrono>
+#include <cstdint>
+#include <fstream>
+#include <iostream>
+#include <ostream>
+#include <stdexcept>
+#include <string>
+#include <string_view>
+#include <tug/Simulation.hpp>
+#include <type_traits>
+#include <vector>
+
+#include "files.hpp"
+
+using namespace tug;
+
+/**
+ * Try to parse an input string into a given template type.
+ */
+template <typename T> 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 <typename T>
+std::vector<T> tokenize(const std::string &input, char delimiter) {
+  std::vector<T> tokens;
+  std::istringstream tokenStream(input);
+  std::string token;
+
+  while (std::getline(tokenStream, token, delimiter)) {
+    tokens.push_back(parseString<T>(token));
+  }
+
+  return tokens;
+}
+
+/**
+ * Opens a file containing CSV and transform it into row-major 2D STL vector.
+ */
+template <typename T>
+std::vector<std::vector<T>> 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<std::vector<T>> csv_data;
+
+  std::string line;
+  while (std::getline(in_file, line)) {
+    csv_data.push_back(tokenize<T>(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 <typename T>
+Eigen::MatrixXd rmVecTocmMatrix(const std::vector<std::vector<T>> &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 <class T, tug::APPROACH app> int doWork(int ngrid) {
+
+  constexpr T dt = 10;
+
+  // create a grid
+  std::string name;
+
+  if constexpr (std::is_same_v<T, double>) {
+    name = "DOUBLE";
+  } else if constexpr (std::is_same_v<T, float>) {
+    name = "FLOAT";
+  } else {
+    name = "unknown";
+  }
+
+  std::cout << name << " grid: " << ngrid << std::endl;
+
+  Grid<T> grid(ngrid, ngrid);
+  // Grid64 grid(ngrid, ngrid);
+
+  // (optional) set the domain, e.g.:
+  grid.setDomain(0.1, 0.1);
+
+  Eigen::MatrixX<T> initConc64 = Eigen::MatrixX<T>::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<T> alphax =
+      Eigen::MatrixX<T>::Constant(ngrid, ngrid, alphax_val); // row,col,value
+  constexpr T alphay_val = 1e-10;
+  Eigen::MatrixX<T> alphay =
+      Eigen::MatrixX<T>::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<T, app>(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<std::chrono::milliseconds>(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<float, tug::BTCS_APPROACH>(n[i]);
+    doWork<double, tug::BTCS_APPROACH>(n[i]);
+    doWork<float, tug::FTCS_APPROACH>(n[i]);
+    doWork<double, tug::FTCS_APPROACH>(n[i]);
+  }
+
+  return 0;
+}

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
+       <build_dir>/naaice`
+
+5.  Finally, run the `naaice` executable to generate the benchmark output:
+    
+        ./naaice
+
+
+## Output Files
+
+
+### `Thomas_<n>.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**!
+

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

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_

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

src/Boundary.cpp

@@ -1,149 +0,0 @@
-#include "TugUtils.hpp"
-
-#include <cstddef>
-#include <stdexcept>
-#include <tug/Boundary.hpp>
-
-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<std::vector<BoundaryElement>>(
-        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<std::vector<BoundaryElement>>(4);
-
-    this->boundaries[BC_SIDE_LEFT] =
-        std::vector<BoundaryElement>(grid.getRow(), BoundaryElement());
-    this->boundaries[BC_SIDE_RIGHT] =
-        std::vector<BoundaryElement>(grid.getRow(), BoundaryElement());
-    this->boundaries[BC_SIDE_TOP] =
-        std::vector<BoundaryElement>(grid.getCol(), BoundaryElement());
-    this->boundaries[BC_SIDE_BOTTOM] =
-        std::vector<BoundaryElement>(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<BoundaryElement>(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<BoundaryElement>(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<BoundaryElement> 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();
-}

src/CMakeLists.txt

@@ -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)

src/Grid.cpp

@@ -1,167 +0,0 @@
-#include "TugUtils.hpp"
-
-#include <iostream>
-#include <tug/Grid.hpp>
-
-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;
-}

src/README.md

@@ -1,36 +0,0 @@
-<h1>src-Directory</h1>
-This is the src-directory that holds the source code to the implementation of the TUG framework.
-
-<pre>
-src/  
-  ├── CMakeFiles/
-  ├── Boundary.cpp  
-  ├── BTCSv2.cpp
-  ├── CMakeLists.txt  
-  ├── FTCS.cpp
-  ├── Grid.cpp 
-  ├── README.md
-  ├── Simulation.cpp
-  └── TugUtils.hpp 
-</pre>
-
-**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** <i>This</i> 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. 
-

src/Schemes.hpp

@@ -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 <tug/Boundary.hpp>
-#include <tug/Grid.hpp>
-
-// entry point; differentiate between 1D and 2D grid
-extern void FTCS(Grid &grid, Boundary &bc, double &timestep, int &numThreads);
-
-// entry point for EigenLU solver; differentiate between 1D and 2D grid
-extern void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads);
-
-// 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_

src/Simulation.cpp

@@ -1,326 +0,0 @@
-#include <cmath>
-#include <cstddef>
-#include <filesystem>
-#include <fstream>
-#include <iostream>
-#include <limits>
-#include <stdexcept>
-#include <string>
-
-#include <tug/Simulation.hpp>
-
-#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<double>::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<std::chrono::milliseconds>(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;
-  }
-}

test/testBoundary.cpp

@@ -6,12 +6,13 @@
 #include <typeinfo>
 
 using namespace std;
+using namespace tug;
 
 TEST_CASE("BoundaryElement") {
 
   SUBCASE("Closed case") {
-    BoundaryElement boundaryElementClosed = BoundaryElement();
+    BoundaryElement boundaryElementClosed = BoundaryElement<double>();
-    CHECK_NOTHROW(BoundaryElement());
+    CHECK_NOTHROW(BoundaryElement<double>());
     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 grid1D = Grid64(10);
-  Grid grid2D = Grid(10, 12);
+  Grid grid2D = Grid64(10, 12);
   Boundary boundary1D = Boundary(grid1D);
   Boundary boundary2D = Boundary(grid2D);
-  vector<BoundaryElement> boundary1DVector(1, BoundaryElement(1.0));
+  vector<BoundaryElement<double>> boundary1DVector(1, BoundaryElement(1.0));
 
   SUBCASE("Boundaries 1D case") {
     CHECK_NOTHROW(Boundary boundary(grid1D));

test/testFTCS.cpp

@@ -1,4 +1,4 @@
-#include <TugUtils.hpp>
+#include <tug/Core/TugUtils.hpp>
 
 #include <doctest/doctest.h>
 #include <limits>

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);

test/testSimulation.cpp

@@ -10,16 +10,16 @@
 
 using namespace Eigen;
 using namespace std;
+using namespace tug;
 
-static Grid setupSimulation(APPROACH approach, double timestep,
+Grid64 setupSimulation(double timestep, int iterations) {
-                            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<double, tug::FTCS_APPROACH>(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<double, tug::FTCS_APPROACH>(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<double, tug::FTCS_APPROACH> sim(grid, boundary);
 
   SUBCASE("default paremeters") { CHECK_EQ(sim.getIterations(), -1); }