diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp
index a20de13..8409279 100644
--- a/include/tug/Advection/Advection.hpp
+++ b/include/tug/Advection/Advection.hpp
@@ -119,6 +119,8 @@ public:
   void run() {
     this->setDomain(velocities.domainX(), velocities.domainY());
 
+    this->applyInnerBoundaries();
+
     if constexpr (hyd_mode == HYDRAULIC_MODE::STEADY_STATE) {
       velocities.run();
     }
diff --git a/include/tug/Advection/Velocities.hpp b/include/tug/Advection/Velocities.hpp
index f62f8d9..fd9580f 100644
--- a/include/tug/Advection/Velocities.hpp
+++ b/include/tug/Advection/Velocities.hpp
@@ -193,6 +193,8 @@ public:
     // if iterations < 1 calculate hydraulic charge until steady state is
     // reached
 
+    this->applyInnerBoundaries();
+
     SimulationInput<T> input = {.concentrations =
                                     this->getConcentrationMatrix(),
                                 .alphaX = this->getPermKX(),
diff --git a/include/tug/Core/BaseSimulation.hpp b/include/tug/Core/BaseSimulation.hpp
index 82a1480..3dc8268 100644
--- a/include/tug/Core/BaseSimulation.hpp
+++ b/include/tug/Core/BaseSimulation.hpp
@@ -64,6 +64,18 @@ private:
   T delta_col;
   T delta_row;
 
+protected:
+  void applyInnerBoundaries() {
+    const auto &inner_bc = boundaryConditions.getInnerBoundaries();
+    if (inner_bc.empty()) {
+      return;
+    }
+
+    for (const auto &[rowcol, value] : inner_bc) {
+      concentrationMatrix(rowcol.first, rowcol.second) = value;
+    }
+  }
+
 public:
   /**
    * @brief Constructs a BaseSimulationGrid from a given RowMajMat object.
diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp
index 28d2db6..4ccf1a4 100644
--- a/include/tug/Core/Numeric/FTCS.hpp
+++ b/include/tug/Core/Numeric/FTCS.hpp
@@ -57,23 +57,6 @@ constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center,
   return 0;
 }
 
-template <typename T>
-static inline void
-checkAndSetConstantInnerCells(const Boundary<T> &bc,
-                              RowMajMatMap<T> &concentrations, std::size_t rows,
-                              std::size_t cols) {
-  const auto &inner_boundaries = bc.getInnerBoundaries();
-  if (inner_boundaries.empty()) {
-    return;
-  }
-
-  for (const auto &[rowcol, value] : inner_boundaries) {
-    const auto &row = rowcol.first;
-    const auto &col = rowcol.second;
-    concentrations(row, col) = value;
-  }
-}
-
 // FTCS solution for 1D grid
 template <class T> static void FTCS_1D(SimulationInput<T> &input) {
   const std::size_t &colMax = input.colMax;
@@ -84,8 +67,6 @@ template <class T> static void FTCS_1D(SimulationInput<T> &input) {
 
   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;
 
@@ -169,9 +150,6 @@ static void FTCS_2D(SimulationInput<T> &input, int numThreads) {
   const T sx = timestep / (deltaCol * deltaCol);
   const T sy = timestep / (deltaRow * deltaRow);
 
-  checkAndSetConstantInnerCells(bc, concentrations_grid, input.rowMax,
-                                input.colMax);
-
   // matrix for concentrations at time t+1
   RowMajMat<T> concentrations_t1 = concentrations_grid;
 
diff --git a/include/tug/Diffusion/Diffusion.hpp b/include/tug/Diffusion/Diffusion.hpp
index e353550..f977514 100644
--- a/include/tug/Diffusion/Diffusion.hpp
+++ b/include/tug/Diffusion/Diffusion.hpp
@@ -347,6 +347,8 @@ public:
 
     auto begin = std::chrono::high_resolution_clock::now();
 
+    this->applyInnerBoundaries();
+
     SimulationInput<T> sim_input = {.concentrations =
                                         this->getConcentrationMatrix(),
                                     .alphaX = this->getAlphaX(),