diff --git a/examples/Advection.cpp b/examples/Advection.cpp
index 7876243..c2ca66a 100644
--- a/examples/Advection.cpp
+++ b/examples/Advection.cpp
@@ -1,3 +1,6 @@
+#include "tug/Advection/Advection.hpp"
+#include "tug/Advection/Velocities.hpp"
+#include "tug/Core/Matrix.hpp"
 #include <Eigen/Eigen>
 #include <iostream>
 #include <tug/tug.hpp>
@@ -6,53 +9,53 @@ using namespace Eigen;
 using namespace tug;
 
 int main(int argc, char *argv[]) {
-  int row = 21;
-  int col = 21;
+  int row = 5;
+  int col = 5;
   // 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);
+  hydHeads(row / 2, col / 2) = 10;
+  concentrations(row / 2, col / 2) = 1;
 
-  Grid64 gridHeads(hydHeads);
-  Grid64 gridConc(concentrations);
-  gridHeads.setDomain(100, 100);
-  gridConc.setDomain(100, 100);
+  Velocities<double, tug::HYDRAULIC_MODE::STEADY_STATE,
+             tug::HYDRAULIC_RESOLVE::EXPLICIT>
+      velocities(hydHeads);
+  velocities.setDomain(100, 100);
+  velocities.setPermKX(RowMajMat<double>::Constant(row, col, 1E-8));
+  velocities.setPermKY(RowMajMat<double>::Constant(row, col, 1E-8));
+  velocities.setEpsilon(1E-8);
+
+  Advection advection = Advection(concentrations, velocities);
+
+  advection.setPorosity(RowMajMat<double>::Constant(row, col, 0.2));
+  advection.setIterations(1);
+  advection.setTimestep(10000);
 
   // create boundaries
-  Boundary bcH = Boundary(gridHeads);
-  bcH.setBoundarySideConstant(BC_SIDE_LEFT, 10);
+  Boundary<double> &bcH = velocities.getBoundaryConditions();
+  bcH.setBoundarySideConstant(BC_SIDE_LEFT, 1);
   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);
+  bcH.setBoundarySideConstant(BC_SIDE_TOP, 1);
+  bcH.setBoundarySideConstant(BC_SIDE_BOTTOM, 1);
+  // bcH.setBoundarySideClosed(BC_SIDE_TOP);
+  // bcH.setBoundarySideClosed(BC_SIDE_BOTTOM);
+  bcH.setInnerBoundary(row / 2, col / 2, 10);
+  //
+  Boundary<double> &bcC = advection.getBoundaryConditions();
+  // 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 << velocities.getConcentrationMatrix() << std::endl;
+  std::cout << velocities.getVelocitiesX() << std::endl;
+  std::cout << velocities.getVelocitiesY() << std::endl;
+
   std::cout << "Advection simulation finished." << std::endl;
 
   std::cout << concentrations << std::endl;
-}
\ No newline at end of file
+}
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 8924b9c..c6f2eb8 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,9 +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(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(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)
diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp
index 855e1b5..ff448c5 100644
--- a/include/tug/Advection/Advection.hpp
+++ b/include/tug/Advection/Advection.hpp
@@ -8,18 +8,16 @@
 #pragma once
 
 #include "tug/Core/Matrix.hpp"
+#include <Eigen/src/Core/Array.h>
+#include <Eigen/src/Core/Matrix.h>
 #include <algorithm>
-#include <filesystem>
-#include <fstream>
-#include <iostream>
+#include <cstddef>
 #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>
@@ -29,7 +27,15 @@
 
 using namespace Eigen;
 namespace tug {
-template <class T> class Advection {
+template <class T, HYDRAULIC_MODE hyd_mode, HYDRAULIC_RESOLVE hyd_resolve>
+class Advection : public BaseSimulationGrid<T> {
+private:
+  T timestep{-1};
+  int numThreads{omp_get_num_procs()};
+
+  Velocities<T, hyd_mode, hyd_resolve> &velocities;
+  RowMajMat<T> porosity;
+
 public:
   /**
    * @brief Construct a new Advection object, used to calculate material
@@ -46,53 +52,45 @@ public:
    * @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) {};
+  Advection(RowMajMat<T> &origin,
+            Velocities<T, hyd_mode, hyd_resolve> &_velocities)
+      : BaseSimulationGrid<T>(origin), velocities(_velocities) {
+    tug_assert(origin.rows() == velocities.getConcentrationMatrix().rows() &&
+                   origin.cols() == velocities.getConcentrationMatrix().cols(),
+               "Advection grid and Velocities must have the same dimensions");
+  };
+
+  Advection(T *data, std::size_t rows, std::size_t cols,
+            Velocities<T, hyd_mode, hyd_resolve> &_velocities)
+      : BaseSimulationGrid<T>(data, rows, cols), velocities(_velocities) {
+    tug_assert(rows == velocities.getConcentrationMatrix().rows() &&
+                   cols == velocities.getConcentrationMatrix().cols(),
+               "Advection grid and Velocities must have the same dimensions");
+  };
 
   /**
    * @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)");
-    }
+  void setPorosity(const RowMajMat<T> &porosity) {
+    tug_assert(porosity.rows() == this->rows() &&
+                   porosity.cols() == this->cols(),
+               "Porosity matrix must have the same dimensions as the grid");
+    tug_assert(porosity.minCoeff() >= 0 && porosity.maxCoeff() <= 1,
+               "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;
-    }
+    tug_assert(timestep > 0, "Timestep must be greater than zero");
+    this->timestep = timestep;
   }
 
   /**
@@ -118,113 +116,54 @@ public:
     }
   }
 
-  /**
-   * @brief Set the option to output the results to a CSV file. Off by default.
-   *
-   *
-   * @param csv_output Valid output option. The following options can be set
-   *                   here:
-   *                     - CSV_OUTPUT_OFF: do not produce csv output
-   *                     - CSV_OUTPUT_ON: produce csv output with last
-   *                       concentration matrix
-    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!");
-    }
+  void run() {
+    this->setDomain(velocities.domainX(), velocities.domainY());
 
-    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();
+    velocities.run();
+
+    for (int i = 0; i < this->getIterations(); i++) {
+      if constexpr (hyd_mode == HYDRAULIC_MODE::TRANSIENT) {
+        velocities.run();
+      }
+      adv();
+    }
   }
 
+private:
   /**
    * @brief function calculating material transport for one timestep
    */
   void adv() {
-    int rows = grid.getRow();
-    int cols = grid.getCol();
-    T volume = grid.getDeltaRow() * grid.getDeltaCol();
+    const std::size_t rows = this->rows();
+    const std::size_t cols = this->cols();
+    const T volume = this->deltaCol() * this->deltaRow();
 
-    RowMajMat<T> &newConcentrations = grid.getConcentrations();
+    RowMajMatMap<T> &newConcentrations = this->getConcentrationMatrix();
+
+    const RowMajMat<T> &outx = this->velocities.getVelocitiesX();
+    const RowMajMat<T> &outy = this->velocities.getVelocitiesY();
+    const Boundary<T> &bc = this->getBoundaryConditions();
 
     // 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);
+    const T maxFx = std::max(abs(outx.maxCoeff()), abs(outx.minCoeff()));
+    const T maxFy = std::max(abs(outy.maxCoeff()), abs(outy.minCoeff()));
+    const T maxF = std::max(maxFx, maxFy);
 
-    if (maxF == 0) {
-      throw std::runtime_error("Division by zero: maxF is zero.");
-    }
+    tug_assert(maxF != 0, "Division by zero: maxF is zero.");
 
-    T cvf = abs((volume * porosity) / maxF);
-    int innerSteps = (int)ceil(timestep / cvf);
-    T innerTimestep = timestep / innerSteps;
+    const RowMajMat<T> volumeTimesPorosity = volume * porosity;
+
+    const T cvf = (volumeTimesPorosity / maxF).maxCoeff();
+    const int innerSteps = (int)ceil(timestep / cvf);
+    const T innerTimestep = timestep / innerSteps;
+
+    const RowMajMat<T> multiplier = volumeTimesPorosity * (1 / innerTimestep);
 
     for (int k = 0; k < innerSteps; k++) {
-      const Eigen::MatrixX<T> oldConcentrations = newConcentrations;
+      const RowMajMat<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)
+#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) {
@@ -318,51 +257,15 @@ public:
           }
         }
       }
-      newConcentrations =
-          oldConcentrations +
-          newConcentrations * (innerTimestep / (porosity * volume));
+
+      for (std::size_t row_i = 0; row_i < rows; row_i++) {
+        for (std::size_t col_i = 0; col_i < cols; col_i++) {
+          newConcentrations(row_i, col_i) =
+              oldConcentrations(row_i, col_i) +
+              newConcentrations(row_i, col_i) * multiplier(row_i, col_i);
+        }
+      }
     }
   }
-
-  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
diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp
index bcb78ba..0c50671 100644
--- a/include/tug/Core/Numeric/FTCS.hpp
+++ b/include/tug/Core/Numeric/FTCS.hpp
@@ -82,14 +82,13 @@ template <class T> static void FTCS_1D(SimulationInput<T> &input) {
   const T &timestep = input.timestep;
 
   RowMajMatMap<T> &concentrations_grid = input.concentrations;
-  // matrix for concentrations at time t+1
-  RowMajMat<T> concentrations_t1 = concentrations_grid;
 
   const auto &alphaX = input.alphaX;
   const auto &bc = input.boundaries;
-
   checkAndSetConstantInnerCells(bc, concentrations_grid, input.rowMax,
                                 input.colMax);
+  // matrix for concentrations at time t+1
+  RowMajMat<T> concentrations_t1 = concentrations_grid;
 
   // only one row in 1D case -> row constant at index 0
   int row = 0;
@@ -99,6 +98,9 @@ template <class T> static void FTCS_1D(SimulationInput<T> &input) {
   // inner cells
   // independent of boundary condition type
   for (int col = 1; col < colMax - 1; col++) {
+    if (inner_bc[col].first) {
+      continue;
+    }
     const T &conc_c = concentrations_grid(row, col);
     const T &conc_left = concentrations_grid(row, col - 1);
     const T &conc_right = concentrations_grid(row, col + 1);
@@ -115,8 +117,8 @@ template <class T> static void FTCS_1D(SimulationInput<T> &input) {
   }
 
   // left boundary; hold column constant at index 0
-  {
-    int col = 0;
+  int col = 0;
+  if (!inner_bc[col].first) {
     const T &conc_c = concentrations_grid(row, col);
     const T &conc_right = concentrations_grid(row, col + 1);
     const T &alpha_c = alphaX(row, col);
@@ -132,8 +134,8 @@ template <class T> static void FTCS_1D(SimulationInput<T> &input) {
   }
 
   // right boundary; hold column constant at max index
-  {
-    int col = colMax - 1;
+  col = colMax - 1;
+  if (!inner_bc[col].first) {
     const T &conc_c = concentrations_grid(row, col);
     const T &conc_left = concentrations_grid(row, col - 1);
     const T &alpha_c = alphaX(row, col);
@@ -161,27 +163,28 @@ static void FTCS_2D(SimulationInput<T> &input, int numThreads) {
   const T &timestep = input.timestep;
 
   RowMajMatMap<T> &concentrations_grid = input.concentrations;
-
-  // matrix for concentrations at time t+1
-  RowMajMat<T> concentrations_t1 = concentrations_grid;
-
   const auto &alphaX = input.alphaX;
   const auto &alphaY = input.alphaY;
   const auto &bc = input.boundaries;
 
+  const T sx = timestep / (deltaCol * deltaCol);
+  const T sy = timestep / (deltaRow * deltaRow);
+
   checkAndSetConstantInnerCells(bc, concentrations_grid, input.rowMax,
                                 input.colMax);
 
-  const T sx = timestep / (deltaCol * deltaCol);
-  const T sy = timestep / (deltaRow * deltaRow);
+  // matrix for concentrations at time t+1
+  RowMajMat<T> concentrations_t1 = concentrations_grid;
 
 #pragma omp parallel for num_threads(numThreads)
   for (std::size_t row_i = 0; row_i < rowMax; row_i++) {
     for (std::size_t col_i = 0; col_i < colMax; col_i++) {
+      if (bc.getInnerBoundary(row_i, col_i).first) {
+        continue;
+      }
       // horizontal change
       T horizontal_change;
       {
-
         const T &conc_c = concentrations_grid(row_i, col_i);
         const T &alpha_c = alphaX(row_i, col_i);
 
diff --git a/include/tug/tug.hpp b/include/tug/tug.hpp
index bc797ae..6770de3 100644
--- a/include/tug/tug.hpp
+++ b/include/tug/tug.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
-// #include <tug/Advection/Advection.hpp>
+#include <tug/Advection/Advection.hpp>
+#include <tug/Advection/Velocities.hpp>
 #include <tug/Boundary.hpp>
 #include <tug/Core/Matrix.hpp>
 #include <tug/Diffusion/Diffusion.hpp>