feat: Add implementation of Advection from Christopher Eschenbach (does not work!)

examples/Advection.cpp

@@ -0,0 +1,59 @@
+#include <Eigen/Eigen>
+#include <iostream>
+#include <tug/Advection.hpp>
+#include <tug/Core/Matrix.hpp>
+
+using namespace Eigen;
+using namespace tug;
+
+int main(int argc, char *argv[]) {
+  int row = 21;
+  int col = 21;
+  // create two grids of equal size, grid1 for hydraulics heads, grid2 for
+  // Concentrations
+
+  RowMajMat<double> hydHeads = RowMajMat<double>::Constant(row, col, 1);
+  RowMajMat<double> concentrations = RowMajMat<double>::Constant(row, col, 0);
+
+  Grid64 gridHeads(hydHeads);
+  Grid64 gridConc(concentrations);
+  gridHeads.setDomain(100, 100);
+  gridConc.setDomain(100, 100);
+
+  // create boundaries
+  Boundary bcH = Boundary(gridHeads);
+  bcH.setBoundarySideConstant(BC_SIDE_LEFT, 10);
+  bcH.setBoundarySideConstant(BC_SIDE_RIGHT, 1);
+  bcH.setBoundarySideClosed(BC_SIDE_TOP);
+  bcH.setBoundarySideClosed(BC_SIDE_BOTTOM);
+  Boundary bcC = Boundary(gridConc);
+  bcC.setBoundarySideConstant(BC_SIDE_LEFT, 0.1);
+  bcC.setBoundarySideConstant(BC_SIDE_RIGHT, 1);
+  bcC.setBoundarySideClosed(BC_SIDE_TOP);
+  bcC.setBoundarySideClosed(BC_SIDE_BOTTOM);
+
+  Velocities velocities = Velocities(gridHeads, bcH);
+  velocities.setOutputCSV(CSV_OUTPUT_ON);
+  velocities.setK(1E-2);
+  velocities.setEpsilon(1E-8);
+  velocities.setInjh(10);
+  velocities.setIterations(0);
+  // calculate steady hydraulic heads
+  velocities.run();
+
+  std::cout << "Velocities simulation finished." << std::endl;
+  std::cout << hydHeads << std::endl;
+
+  // set true for steady case
+  Advection advection = Advection(velocities, gridConc, bcC, true);
+  advection.setPorosity(0.2);
+  advection.setIterations(21);
+  // set timestep close almost exactly to clf to test advection
+  advection.setTimestep(5039.05);
+  //   velocities.setOutputCSV(CSV_OUTPUT_VERBOSE);
+  advection.run();
+
+  std::cout << "Advection simulation finished." << std::endl;
+
+  std::cout << concentrations << std::endl;
+}

examples/CMakeLists.txt

@@ -1,7 +1,9 @@
 add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp)
 add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp)
 add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp)
+add_executable(Advection Advection.cpp)
 
 target_link_libraries(BTCS_2D_proto_example tug)
 target_link_libraries(FTCS_2D_proto_closed_mdl tug)
 target_link_libraries(FTCS_2D_proto_example_mdl tug)
+target_link_libraries(Advection tug)

include/tug/Advection.hpp

@@ -0,0 +1,368 @@
+/**
+ * @file Advection.hpp
+ * @brief API of Advection class, holding information for a simulation of
+ * advection. Holds a predifined Grid object, Boundary object and Velocities
+ * object
+ */
+
+#pragma once
+
+#include "tug/Core/Matrix.hpp"
+#include <algorithm>
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+#include <stdexcept>
+#include <stdlib.h>
+#include <string>
+#include <tug/Boundary.hpp>
+#include <tug/Grid.hpp>
+
+#include <Eigen/Eigen>
+#include <Eigen/Sparse>
+#include <tug/Core/Numeric/BTCS.hpp>
+#include <tug/Core/Numeric/FTCS.hpp>
+#include <tug/Core/TugUtils.hpp>
+#include <tug/Diffusion.hpp>
+
+#include <tug/Core/Velocities.hpp>
+
+using namespace Eigen;
+namespace tug {
+template <class T> class Advection {
+public:
+  /**
+   * @brief Construct a new Advection object, used to calculate material
+   * transport. A timestep and number of iterations must be set. A transient
+   * case can be selected by initializing Steady=false. With each timestep the
+   * Velocities object will also be updated.
+   * A steady case can be selected by initializing Steady=true. The velocities
+   * object will not be updated. Velocities can be simulated to convergence
+   * beforehand. Porosity can be set, the default is 1. CSV Output is off by
+   * default.
+   *
+   * @param grid Valid grid object
+   * @param bc Valid Boundary object
+   * @param Steady Used to choose between Steady and Transient case. Either true
+   * or false
+   */
+  Advection(Velocities<T> &_velocities, Grid<T> &_grid, Boundary<T> &_bc,
+            bool Steady)
+      : velocities(_velocities), grid(_grid), bc(_bc),
+        outx(_velocities.getOutx()), outy(_velocities.getOuty()),
+        Steady(Steady) {};
+
+  /**
+   * @brief Sets the porosity of the medium
+   *
+   * @param porosity new porosity value
+   */
+  void setPorosity(T porosity) {
+    if (porosity < 0 || porosity > 1) {
+      throw std::invalid_argument(
+          "Porosity must be a value between 0 and 1 (inclusive)");
+    }
+    this->porosity = porosity;
+  }
+
+  /**
+   * @brief Set the desired iterations to be calculated. A value greater
+   *        than zero must be specified here. Setting iterations is required.
+   *
+   * @param iterations Number of iterations to be simulated.
+   */
+  void setIterations(int iterations) {
+    if (iterations <= 0) {
+      throw std::invalid_argument(
+          "Number of iterations must be greater than zero. Provided value: " +
+          std::to_string(iterations));
+    }
+    this->iterations = iterations;
+  };
+
+  /**
+   * @brief Set the size of the timestep. Must be greater than zero
+   *
+   * @param timestep Size of the timestep
+   */
+  void setTimestep(T timestep) {
+    if (timestep <= 0) {
+      throw std::invalid_argument(
+          "Timestep must be greater than zero. Provided value: " +
+          std::to_string(timestep));
+    } else {
+      this->timestep = timestep;
+    }
+  }
+
+  /**
+   * @brief Set the number of desired openMP Threads.
+   *
+   * @param num_threads Number of desired threads. Must have a value between
+   *                    1 and the maximum available number of processors. The
+   * maximum number of processors is set as the default case during Advection
+   * construction.
+   */
+  void setNumberThreads(int num_threads) {
+    if (num_threads > 0 && num_threads <= omp_get_num_procs()) {
+      this->numThreads = num_threads;
+    } else {
+      int maxThreadNumber = omp_get_num_procs();
+      if (num_threads > maxThreadNumber) {
+        throw std::invalid_argument(
+            "Number of threads exceeds the number of processor cores (" +
+            std::to_string(maxThreadNumber) + ").");
+      } else {
+        throw std::invalid_argument("Number of threads is less than 1.");
+      }
+    }
+  }
+
+  /**
+   * @brief Set the option to output the results to a CSV file. Off by default.
+   *
+   *
+   * @param csv_output Valid output option. The following options can be set
+   *                   here:
+   *                     - CSV_OUTPUT_OFF: do not produce csv output
+   *                     - CSV_OUTPUT_ON: produce csv output with last
+   *                       concentration matrix
+    if (csv_output < CSV_OUTPUT_OFF || csv_output > CSV_OUTPUT_VERBOSE) {
+      throw std::invalid_argument("Invalid CSV output option given: " +
+  std::to_string(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 Creates a CSV file with a name containing the current simulation
+   *        parameters. If the data name already exists, an additional counter
+   * is appended to the name. The name of the file is built up as follows:
+   *        <Information contained in file> + <number rows> + <number columns> +
+   * <number of iterations>+<counter>.csv
+   *
+   * @return string Filename with configured simulation parameters.
+   */
+  std::string createCSVfile(std::string Type) const {
+    std::ofstream file;
+    int appendIdent = 0;
+    std::string appendIdentString;
+
+    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 =
+        Type + "_" + row + "_" + col + "_" + numIterations + ".csv";
+
+    while (std::filesystem::exists(filename)) {
+      appendIdent += 1;
+      appendIdentString = std::to_string(appendIdent);
+      filename = Type + "_" + row + "_" + col + "_" + numIterations + "-" +
+                 appendIdentString + ".csv";
+    }
+
+    file.open(filename);
+    if (!file) {
+      throw std::runtime_error("Failed to open file: " + filename);
+    }
+
+    file.close();
+
+    return filename;
+  }
+  /**
+   * @brief Writes the currently calculated Concentration values of the grid
+   *        into the CSV file with the passed filename.
+   *
+   * @param filename Name of the file to which the Concentration values are
+   *                 to be written.
+   */
+  void printConcentrationCSV(const std::string &filename) const {
+    std::ofstream file;
+
+    file.open(filename, std::ios_base::app);
+    if (!file) {
+      throw std::runtime_error("Failed to open file: " + filename);
+    }
+
+    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 function calculating material transport for one timestep
+   */
+  void adv() {
+    int rows = grid.getRow();
+    int cols = grid.getCol();
+    T volume = grid.getDeltaRow() * grid.getDeltaCol();
+
+    RowMajMat<T> &newConcentrations = grid.getConcentrations();
+
+    // Calculate Courant-Levy-Frederich condition
+    T maxFx = std::max(abs(outx.maxCoeff()), abs(outx.minCoeff()));
+    T maxFy = std::max(abs(outy.maxCoeff()), abs(outy.minCoeff()));
+    T maxF = std::max(maxFx, maxFy);
+
+    if (maxF == 0) {
+      throw std::runtime_error("Division by zero: maxF is zero.");
+    }
+
+    T cvf = abs((volume * porosity) / maxF);
+    int innerSteps = (int)ceil(timestep / cvf);
+    T innerTimestep = timestep / innerSteps;
+
+    for (int k = 0; k < innerSteps; k++) {
+      const Eigen::MatrixX<T> oldConcentrations = newConcentrations;
+// Calculate sum of incoming/outgoing Flow*Concentration in x-direction in each
+// cell
+#pragma omp parallel for num_threads(numThreads) schedule(static)
+      for (int i = 0; i < rows; i++) {
+        for (int j = 0; j < cols + 1; j++) {
+          if (j == 0) {
+            if (bc.getBoundaryElementType(BC_SIDE_LEFT, i) != BC_TYPE_CLOSED) {
+              if (outx(i, j) > 0) {
+                // outx positive -> flow from border to cell i,j
+                newConcentrations(i, j) +=
+                    outx(i, j) * bc.getBoundaryElementValue(BC_SIDE_LEFT, i);
+              } else if (outx(i, j) < 0) {
+                // outx negative -> flow from i,j towards border
+                newConcentrations(i, j) += outx(i, j) * oldConcentrations(i, j);
+              }
+            }
+          } else if (j == cols) {
+            if (bc.getBoundaryElementType(BC_SIDE_RIGHT, i) != BC_TYPE_CLOSED) {
+              if (outx(i, j) > 0) {
+                // outx positive-> flow from i,j-1 towards border
+                newConcentrations(i, j - 1) -=
+                    outx(i, j) * oldConcentrations(i, j - 1);
+              } else if (outx(i, j) < 0) {
+                // outx negative -> flow from border to cell i,j-1
+                newConcentrations(i, j - 1) -=
+                    outx(i, j) * bc.getBoundaryElementValue(BC_SIDE_LEFT, i);
+              }
+            }
+          }
+          // flow between inner cells
+          else {
+            // outx positive -> flow from cell i,j-1 towards cell i,j
+            if (outx(i, j) > 0) {
+              newConcentrations(i, j - 1) -=
+                  outx(i, j) * oldConcentrations(i, j - 1);
+              newConcentrations(i, j) +=
+                  outx(i, j) * oldConcentrations(i, j - 1);
+            }
+            // outx negative -> flow from cell i,j toward cell i,j-1
+            else if (outx(i, j) < 0) {
+              newConcentrations(i, j - 1) -=
+                  outx(i, j) * oldConcentrations(i, j);
+              newConcentrations(i, j) += outx(i, j) * oldConcentrations(i, j);
+            }
+          }
+        }
+      }
+// calculate sum in y-direction
+// parallelize outer loop over columns to ensure thread-safety, each thread only
+// modifies cells within its column
+#pragma omp parallel for num_threads(numThreads)
+      for (int j = 0; j < cols; j++) {
+        for (int i = 0; i < rows + 1; i++) {
+          if (i == 0) {
+            if (bc.getBoundaryElementType(BC_SIDE_TOP, j) != BC_TYPE_CLOSED) {
+              if (outy(i, j) > 0) {
+                // outy positive -> flow from border to cell i,j
+                newConcentrations(i, j) +=
+                    outy(i, j) * bc.getBoundaryElementValue(BC_SIDE_TOP, j);
+              } else if (outy(i, j) < 0) {
+                // outy negative -> flow from i,j towards border
+                newConcentrations(i, j) += outy(i, j) * oldConcentrations(i, j);
+              }
+            }
+          } else if (i == rows) {
+            if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, j) !=
+                BC_TYPE_CLOSED) {
+              if (outy(i, j) > 0) {
+                // outy positive-> flow from i-1,j towards border
+                newConcentrations(i - 1, j) -=
+                    outy(i, j) * oldConcentrations(i - 1, j);
+              } else if (outy(i, j) < 0) {
+                // outy negative -> flow from border to cell i,j-1
+                newConcentrations(i - 1, j) -=
+                    outy(i, j) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, j);
+              }
+            }
+          }
+          // flow between inner cells
+          else {
+            // outy positive -> flow from cell i-1,j towards cell i,j
+            if (outy(i, j) > 0) {
+              newConcentrations(i - 1, j) -=
+                  outy(i, j) * oldConcentrations(i - 1, j);
+              newConcentrations(i, j) +=
+                  outy(i, j) * oldConcentrations(i - 1, j);
+            }
+            // outy negative -> flow from cell i,j toward cell i-1,j
+            else if (outy(i, j) < 0) {
+              newConcentrations(i - 1, j) -=
+                  outy(i, j) * oldConcentrations(i, j);
+              newConcentrations(i, j) += outy(i, j) * oldConcentrations(i, j);
+            }
+          }
+        }
+      }
+      newConcentrations =
+          oldConcentrations +
+          newConcentrations * (innerTimestep / (porosity * volume));
+    }
+  }
+
+  void run() {
+    std::string filename;
+    if (csv_output >= CSV_OUTPUT_ON) {
+      filename = createCSVfile("Concentrations");
+    }
+
+    if (Steady == false) {
+      velocities.setTimestep(timestep);
+      velocities.setIterations(1);
+    }
+
+    for (int i = 0; i < iterations; i++) {
+      if (csv_output >= CSV_OUTPUT_VERBOSE) {
+        printConcentrationCSV(filename);
+      }
+      // if steady==false update charge and velocities with equal timestep
+      if (Steady == false) {
+        velocities.run();
+      }
+      adv();
+    }
+
+    if (csv_output >= CSV_OUTPUT_ON) {
+      printConcentrationCSV(filename);
+    }
+  }
+
+private:
+  Grid<T> &grid;
+  Boundary<T> &bc;
+  Velocities<T> &velocities;
+  bool Steady{true};
+  int iterations{-1};
+  int innerIterations{1};
+  T timestep{-1};
+  int numThreads{omp_get_num_procs()};
+  T porosity{1};
+  const Eigen::MatrixX<T> &outx;
+  const Eigen::MatrixX<T> &outy;
+  CSV_OUTPUT csv_output{CSV_OUTPUT_OFF};
+};
+} // namespace tug

include/tug/Core/Velocities.hpp

@@ -0,0 +1,488 @@
+/**
+ * @file Velocities.hpp
+ * @brief API of Velocities class, holding information for a simulation of
+ * Hydraulic Charge and Darcy-Velocities. Holds a predifined Grid object and
+ * Boundary object.
+ *
+ */
+
+#pragma once
+
+#include "tug/Core/Matrix.hpp"
+#include <algorithm>
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+#include <stdexcept>
+#include <stdlib.h>
+#include <string>
+#include <tug/Boundary.hpp>
+#include <tug/Grid.hpp>
+#include <vector>
+
+#include <Eigen/Eigen>
+#include <Eigen/Sparse>
+#include <tug/Core/Numeric/BTCS.hpp>
+#include <tug/Core/Numeric/FTCS.hpp>
+#include <tug/Core/TugUtils.hpp>
+#include <tug/Diffusion.hpp>
+
+using namespace Eigen;
+namespace tug {
+template <class T> class Velocities {
+public:
+  /**
+   * @brief Construct a new Velocities object, used to calculate Hydraulic
+   * Charge and Darcy-Velocities. A timestep and a number of iterations can be
+   * set. By default iterations is set to -1. If the number of iterations is set
+   * to a value below 1 the simulation will run until the Hydraulic Charge
+   * converges. The Epsilon value to check convergence can be set, the default
+   * is 1E-5. CSV Output is off by default.
+   *
+   * @param grid Valid grid object
+   * @param bc Valid boundary condition object
+   */
+  Velocities(Grid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc) {
+    outx = MatrixX<T>::Constant(grid.getRow(), grid.getCol() + 1, 0);
+    outy = MatrixX<T>::Constant(grid.getRow() + 1, grid.getCol(), 0);
+    center = std::make_pair(grid.getRow() / 2, grid.getCol() / 2);
+  };
+
+  /**
+   * @brief Sets a fixed, constant hydraulic charge at domain center.
+   *
+   * @param inj_h fixed hydraulic charge at domain center.
+   */
+  void setInjh(T inj_h) {
+    if (inj_h >= 0) {
+      this->inj_h = inj_h;
+      RowMajMat<T> &concentrations = grid.getConcentrations();
+      concentrations(center.first, center.second) = inj_h;
+      injhIsSet = true;
+    } else {
+      throw std::invalid_argument("Fixed hydraulic charge can not be negative");
+    }
+  };
+
+  /**
+   * @brief Sets a constant permeability coefficient
+   * @param K constant permeability coefficient
+   */
+  void setK(T K) {
+    this->K = K;
+    MatrixXd alphax = MatrixXd::Constant(grid.getRow(), grid.getCol(), K);
+    MatrixXd alphay = MatrixXd::Constant(grid.getRow(), grid.getCol(), K);
+    grid.setAlpha(alphax, alphay);
+  };
+
+  /**
+   * @brief Set the epsilon value, the relativ error allowed for convergence
+   *
+   * @param epsilon the new epsilon value
+   */
+  void setEpsilon(T epsilon) {
+    if (0 <= epsilon && epsilon < 1) {
+      this->epsilon = epsilon;
+    } else {
+      throw std::invalid_argument(
+          "Relative Error epsilon must be between 0 and 1");
+    }
+  }
+
+  /**
+   * @brief Set the timestep per iteration
+   *
+   * @param timestep timestep per iteration
+   */
+  void setTimestep(T timestep) {
+    if (timestep <= 0) {
+      throw std::invalid_argument("Timestep must be greater than zero");
+    }
+    this->timestep = timestep;
+    const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
+    const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow();
+    const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare);
+    T cfl = minDeltaSquare / (4 * K);
+    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 << "FTCS" << "_" << "2D" << " :: Required "
+                << this->innerIterations
+                << " inner iterations with dt=" << this->timestep << std::endl;
+    } else {
+      this->innerIterations = 1;
+    }
+  };
+
+  /**
+   * @brief Set the number of iterations. If set to a number smaller than 1,
+   * calculation will terminate at convergence
+   *
+   * @param iterations Number of desired iterations
+   */
+
+  void setIterations(int iterations) { this->iterations = iterations; }
+
+  /**
+   * @brief Set the number of desired openMP Threads.
+   *
+   * @param num_threads Number of desired threads. Must have a value between
+   *                    1 and the maximum available number of processors. The
+   * maximum number of processors is set as the default case during Velocities
+   * construction.
+   */
+  void setNumberThreads(int num_threads) {
+    if (num_threads > 0 && num_threads <= omp_get_num_procs()) {
+      this->numThreads = num_threads;
+    } 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 Set the option to output the results to a CSV file. Off by default.
+   *
+   *
+   * @param csv_output Valid output option. The following options can be set
+   *                   here:
+   *                     - CSV_OUTPUT_OFF: do not produce csv output
+   *                     - CSV_OUTPUT_ON: produce csv output with last
+   *                       charge and velocities matrizes
+   *                     - CSV_OUTPUT_VERBOSE: produce csv output with all
+   *                       charge matrizes and and last velocities matrix
+   */
+  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;
+    if (csv_output >= CSV_OUTPUT_ON) {
+      filename1 = createCSVfile("Charge");
+      filename2 = createCSVfile("outx");
+      filename3 = createCSVfile("outy");
+    }
+  }
+
+  /**
+   * @brief Creates a CSV file with a name containing the current simulation
+   *        parameters. If the data name already exists, an additional counter
+   * is appended to the name. The name of the file is built up as follows:
+   *        <Information contained in file> + <number rows> + <number columns> +
+   * <number of iterations>+<counter>.csv
+   *
+   * @return string Filename with configured simulation parameters.
+   */
+  std::string createCSVfile(std::string Type) const {
+    std::ofstream file;
+    int appendIdent = 0;
+    std::string appendIdentString;
+
+    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 =
+        Type + "_" + row + "_" + col + "_" + numIterations + ".csv";
+
+    while (std::filesystem::exists(filename)) {
+      appendIdent += 1;
+      appendIdentString = std::to_string(appendIdent);
+      filename = Type + "_" + row + "_" + col + "_" + numIterations + "-" +
+                 appendIdentString + ".csv";
+    }
+
+    file.open(filename);
+    if (!file) {
+      exit(1);
+    }
+
+    file.close();
+
+    return filename;
+  }
+
+  /**
+   * @brief Writes the currently calculated Hydraulic Charge values of the grid
+   *        into the CSV file with the passed filename.
+   *
+   * @param filename Name of the file to which the Hydraulic Charge values are
+   *                 to be written.
+   */
+  void printChargeCSV(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 Reads a matrix stored in a CSV file and uses it for values of
+   * Hydraulic Heads, Matrix and grid must be of equal size
+   *
+   * @param filename name of the CSV file
+   */
+  void readChargeCSV(std::string filename) {
+    std::ifstream file(filename);
+    std::string line;
+    Eigen::MatrixXd matrix;
+
+    if (file.is_open()) {
+      while (std::getline(file, line)) {
+        std::istringstream iss(line);
+        double value;
+        std::vector<double> row;
+        while (iss >> value) {
+          row.push_back(value);
+        }
+        if (!row.empty()) {
+          if (matrix.rows() == 0) {
+            matrix.resize(1, row.size());
+          } else {
+            matrix.conservativeResize(matrix.rows() + 1, Eigen::NoChange);
+          }
+          matrix.row(matrix.rows() - 1) =
+              Eigen::VectorXd::Map(row.data(), row.size());
+        }
+      }
+      file.close();
+    } else {
+      std::cerr << "Unable to open file: " << filename << std::endl;
+    }
+    if (matrix.rows() == grid.getRow() && matrix.cols() == grid.getCol()) {
+      grid.setConcentrations(matrix);
+      velocities();
+      if (csv_output > CSV_OUTPUT_OFF) {
+        printChargeCSV(filename1);
+        printVelocitiesCSV(filename2, filename3);
+      }
+    } else {
+      std::cerr << "gridsize and size of stored matrix dont align\n";
+    }
+  }
+
+  /**
+   * @brief Writes the current Darcy-velocities into a CSV file
+   *
+   * @param filenamex Name of the file to which velocities in direction x are
+   * written
+   * @param filenamey Name of the file to which velocities in direction y are
+   * written
+   */
+  void printVelocitiesCSV(const std::string &filenamex,
+                          const std::string &filenamey) const {
+    std::ofstream filex;
+    std::ofstream filey;
+
+    filex.open(filenamex, std::ios_base::app);
+    if (!filex) {
+      exit(1);
+    }
+    filey.open(filenamey, std::ios_base::app);
+    if (!filey) {
+      exit(1);
+    }
+
+    Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols);
+    filex << outx.format(do_not_align) << std::endl;
+    filex << std::endl << std::endl;
+    filex.close();
+
+    filey << outy.format(do_not_align) << std::endl;
+    filey << std::endl << std::endl;
+    filey.close();
+  }
+
+  /**
+   * @brief Calculate the new hydraulic charge using FTCS
+   */
+  void hydraulic_charge() {
+    FTCS_2D(this->grid, this->bc, this->timestep, this->numThreads);
+    if (injhIsSet == true) {
+      RowMajMat<T> &concentrations = grid.getConcentrations();
+      concentrations(center.first, center.second) = inj_h;
+    }
+  };
+
+  /**
+   * @brief checks if the matrix of Hydraulic Heads has converged
+   *
+   * @return bool true if for all corresponding cells of the matrices,
+   * containing old and new Charge values, the relative error is below the
+   * selected Epsilon
+   */
+  bool checkConvergance(Eigen::MatrixX<T> oldHeads,
+                        Eigen::MatrixX<T> newHeads) {
+    for (int i = 0; i < grid.getRow(); i++) {
+      for (int j = 0; j < grid.getCol(); j++) {
+        if (newHeads(i, j) != 0) {
+          if (abs((oldHeads(i, j) - newHeads(i, j)) / newHeads(i, j)) >
+              epsilon) {
+            return false;
+          }
+        }
+      }
+    }
+    return true;
+  }
+
+  /**
+   * @brief Update the matrices containing Darcy velocities in x and y
+   * directions
+   */
+  void velocities() {
+    int rows = grid.getRow();
+    int cols = grid.getCol();
+    float dx = grid.getDeltaRow();
+    float dy = grid.getDeltaCol();
+    Eigen::MatrixX<T> concentrations = grid.getConcentrations();
+// calculate outx
+#pragma omp parallel for num_threads(numThreads)
+    for (int i = 0; i < rows; i++) {
+      for (int j = 0; j < cols + 1; j++) {
+        if (j == 0) {
+          if (bc.getBoundaryElementType(BC_SIDE_LEFT, i) == BC_TYPE_CLOSED) {
+            outx(i, j) = 0;
+          } else {
+            outx(i, j) = -K *
+                         (concentrations(i, j) -
+                          bc.getBoundaryElementValue(BC_SIDE_LEFT, i)) /
+                         (dx / 2);
+          }
+        } else if (j == cols) {
+          if (bc.getBoundaryElementType(BC_SIDE_RIGHT, i) == BC_TYPE_CLOSED) {
+            outx(i, j) = 0;
+          } else {
+            outx(i, j) = -K *
+                         (bc.getBoundaryElementValue(BC_SIDE_RIGHT, i) -
+                          concentrations(i, j - 1)) /
+                         (dx / 2);
+          }
+
+        } else {
+          outx(i, j) =
+              -K * (concentrations(i, j) - concentrations(i, j - 1)) / dx;
+        }
+      }
+    }
+// calculate outy
+#pragma omp parallel for num_threads(numThreads)
+    for (int i = 0; i < rows + 1; i++) {
+      for (int j = 0; j < cols; j++) {
+        if (i == 0) {
+          if (bc.getBoundaryElementType(BC_SIDE_TOP, j) == BC_TYPE_CLOSED) {
+            outy(i, j) = 0;
+          } else {
+            outy(i, j) = -K *
+                         (concentrations(i, j) -
+                          bc.getBoundaryElementValue(BC_SIDE_TOP, j)) /
+                         (dy / 2);
+          }
+        } else if (i == rows) {
+          if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, j) == BC_TYPE_CLOSED) {
+            outy(i, j) = 0;
+          } else {
+            outy(i, j) = -K *
+                         (bc.getBoundaryElementValue(BC_SIDE_BOTTOM, j) -
+                          concentrations(i - 1, j)) /
+                         (dy / 2);
+          }
+        } else {
+          outy(i, j) =
+              -K * (concentrations(i, j) - concentrations(i - 1, j)) / dy;
+        }
+      }
+    }
+  };
+
+  /**
+   * @brief Getter function for outx, the matrix containing velocities in
+   * x-Direction; returns a reference to outx
+   *
+   * */
+  const Eigen::MatrixX<T> &getOutx() const { return outx; }
+
+  /**
+   * @brief Getter function for outy, the matrix containing velocities in
+   * y-Direction; return a reference to outy
+   */
+  const Eigen::MatrixX<T> &getOuty() const { return outy; }
+
+  /**
+   * @brief Simulation of hydraulic charge either until convergence,
+   * or for a number of selected timesteps. Calculation of Darcy-velocities.
+   */
+  void run() {
+    // if iterations < 1 calculate hydraulic charge until steady state is
+    // reached
+    if (iterations < 1) {
+      // Calculate largest possible timestep, depending on K and gridsize
+      const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
+      const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow();
+      const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare);
+      setTimestep(minDeltaSquare / (4 * K));
+
+      Eigen::MatrixX<T> oldConcentrations;
+      do {
+        oldConcentrations = grid.getConcentrations();
+        if (csv_output >= CSV_OUTPUT_VERBOSE) {
+          printChargeCSV(filename1);
+        }
+        for (int i = 0; i < (grid.getRow() + grid.getCol() - 2); i++) {
+          hydraulic_charge();
+        }
+      } while (checkConvergance(oldConcentrations, grid.getConcentrations()) ==
+               false);
+    }
+    // if iterations >= 1 calculate hydraulice charge for a given number of
+    // iterations
+    else {
+      if (timestep == -1) {
+        throw_invalid_argument("Timestep is not set");
+      }
+      for (int i = 0; i < iterations * innerIterations; i++) {
+        hydraulic_charge();
+      }
+    }
+
+    velocities();
+
+    if (csv_output > CSV_OUTPUT_OFF) {
+      printChargeCSV(filename1);
+      printVelocitiesCSV(filename2, filename3);
+    }
+  };
+
+private:
+  int iterations{-1};
+  int innerIterations{1};
+  bool injhIsSet{false};
+  T timestep{-1};
+  T inj_h{1};
+  T K{1};
+  T epsilon{1E-5};
+  int numThreads{omp_get_num_procs()};
+  Grid<T> &grid;
+  Boundary<T> &bc;
+  CSV_OUTPUT csv_output{CSV_OUTPUT_OFF};
+  Eigen::MatrixX<T> outx;
+  Eigen::MatrixX<T> outy;
+  std::pair<int, int> center;
+  std::string filename1;
+  std::string filename2;
+  std::string filename3;
+};
+} // namespace tug