Add constant inner cell concentration with test cases

include/tug/Core/BTCS.hpp

@@ -16,6 +16,7 @@
 #include <tug/Boundary.hpp>
 #include <tug/Grid.hpp>
 #include <utility>
+#include <vector>
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -52,7 +53,8 @@ template <class T>
 static Eigen::SparseMatrix<T>
 createCoeffMatrix(const Eigen::MatrixX<T> &alpha,
                   const std::vector<BoundaryElement<T>> &bcLeft,
-                  const std::vector<BoundaryElement<T>> &bcRight, int numCols,
+                  const std::vector<BoundaryElement<T>> &bcRight,
+                  const std::vector<std::pair<bool, T>> &inner_bc, int numCols,
                   int rowIndex, T sx) {
 
   // square matrix of column^2 dimension for the coefficients
@@ -60,7 +62,9 @@ createCoeffMatrix(const Eigen::MatrixX<T> &alpha,
   cm.reserve(Eigen::VectorXi::Constant(numCols, 3));
 
   // left column
-
+  if (inner_bc[0].first) {
+    cm.insert(0, 0) = 1;
+  } else {
     switch (bcLeft[rowIndex].getType()) {
     case BC_TYPE_CONSTANT: {
       auto [centerCoeffTop, rightCoeffTop] =
@@ -81,10 +85,15 @@ createCoeffMatrix(const Eigen::MatrixX<T> &alpha,
           "Undefined Boundary Condition Type somewhere on Left or Top!");
     }
     }
+  }
 
   // inner columns
   int n = numCols - 1;
   for (int i = 1; i < n; i++) {
+    if (inner_bc[i].first) {
+      cm.insert(i, i) = 1;
+      continue;
+    }
     cm.insert(i, i - 1) =
         -sx * calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i));
     cm.insert(i, i) =
@@ -96,7 +105,9 @@ createCoeffMatrix(const Eigen::MatrixX<T> &alpha,
   }
 
   // right column
-
+  if (inner_bc[n].first) {
+    cm.insert(n, n) = 1;
+  } else {
     switch (bcRight[rowIndex].getType()) {
     case BC_TYPE_CONSTANT: {
       auto [centerCoeffBottom, leftCoeffBottom] = calcBoundaryCoeffConstant(
@@ -106,8 +117,8 @@ createCoeffMatrix(const Eigen::MatrixX<T> &alpha,
       break;
     }
     case BC_TYPE_CLOSED: {
-    auto [centerCoeffBottom, leftCoeffBottom] =
+      auto [centerCoeffBottom, leftCoeffBottom] = calcBoundaryCoeffClosed(
-        calcBoundaryCoeffClosed(alpha(rowIndex, n), alpha(rowIndex, n - 1), sx);
+          alpha(rowIndex, n), alpha(rowIndex, n - 1), sx);
       cm.insert(n, n - 1) = leftCoeffBottom;
       cm.insert(n, n) = centerCoeffBottom;
       break;
@@ -117,6 +128,7 @@ createCoeffMatrix(const Eigen::MatrixX<T> &alpha,
           "Undefined Boundary Condition Type somewhere on Right or Bottom!");
     }
     }
+  }
 
   cm.makeCompressed(); // important for Eigen solver
 
@@ -155,6 +167,7 @@ createSolutionVector(const Eigen::MatrixX<T> &concentrations,
                      const std::vector<BoundaryElement<T>> &bcRight,
                      const std::vector<BoundaryElement<T>> &bcTop,
                      const std::vector<BoundaryElement<T>> &bcBottom,
+                     const std::vector<std::pair<bool, T>> &inner_bc,
                      int length, int rowIndex, T sx, T sy) {
 
   Eigen::VectorX<T> sv(length);
@@ -163,6 +176,10 @@ createSolutionVector(const Eigen::MatrixX<T> &concentrations,
   // inner rows
   if (rowIndex > 0 && rowIndex < numRows - 1) {
     for (int i = 0; i < length; i++) {
+      if (inner_bc[i].first) {
+        sv(i) = inner_bc[i].second;
+        continue;
+      }
       sv(i) =
           sy *
               calcAlphaIntercell(alphaY(rowIndex, i), alphaY(rowIndex + 1, i)) *
@@ -181,6 +198,10 @@ createSolutionVector(const Eigen::MatrixX<T> &concentrations,
   // first row
   else if (rowIndex == 0) {
     for (int i = 0; i < length; i++) {
+      if (inner_bc[i].first) {
+        sv(i) = inner_bc[i].second;
+        continue;
+      }
       switch (bcTop[i].getType()) {
       case BC_TYPE_CONSTANT: {
         sv(i) = calcExplicitConcentrationsBoundaryConstant(
@@ -204,6 +225,10 @@ createSolutionVector(const Eigen::MatrixX<T> &concentrations,
   // last row
   else if (rowIndex == numRows - 1) {
     for (int i = 0; i < length; i++) {
+      if (inner_bc[i].first) {
+        sv(i) = inner_bc[i].second;
+        continue;
+      }
       switch (bcBottom[i].getType()) {
       case BC_TYPE_CONSTANT: {
         sv(i) = calcExplicitConcentrationsBoundaryConstant(
@@ -226,13 +251,14 @@ createSolutionVector(const Eigen::MatrixX<T> &concentrations,
 
   // first column -> additional fixed concentration change from perpendicular
   // dimension in constant bc case
-  if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) {
+  if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT && !inner_bc[0].first) {
     sv(0) += 2 * sx * alphaX(rowIndex, 0) * bcLeft[rowIndex].getValue();
   }
 
   // last column -> additional fixed concentration change from perpendicular
   // dimension in constant bc case
-  if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) {
+  if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT &&
+      !inner_bc[length - 1].first) {
     sv(length - 1) +=
         2 * sx * alphaX(rowIndex, length - 1) * bcRight[rowIndex].getValue();
   }
@@ -341,17 +367,21 @@ static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
   const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
   const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
 
+  const auto inner_bc = bc.getInnerBoundaryRow(0);
+
   Eigen::MatrixX<T> concentrations = grid.getConcentrations();
   int rowIndex = 0;
-  A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex,
+  A = createCoeffMatrix(alpha, bcLeft, bcRight, inner_bc, length, rowIndex,
                         sx); // this is exactly same as in 2D
   for (int i = 0; i < length; i++) {
     b(i) = concentrations(0, i);
   }
-  if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) {
+  if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT &&
+      !inner_bc[0].first) {
     b(0) += 2 * sx * alpha(0, 0) * bcLeft[0].getValue();
   }
-  if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) {
+  if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT &&
+      !inner_bc[length - 1].first) {
     b(length - 1) += 2 * sx * alpha(0, length - 1) * bcRight[0].getValue();
   }
 
@@ -394,10 +424,11 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
 
 #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
   for (int i = 0; i < rowMax; i++) {
+    auto inner_bc = bc.getInnerBoundaryRow(i);
 
-    A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx);
+    A = createCoeffMatrix(alphaX, bcLeft, bcRight, inner_bc, colMax, i, sx);
     b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
-                             bcTop, bcBottom, colMax, i, sx, sy);
+                             bcTop, bcBottom, inner_bc, colMax, i, sx, sy);
 
     row_t1 = solverFunc(A, b);
 
@@ -411,10 +442,11 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
 
 #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
   for (int i = 0; i < colMax; i++) {
+    auto inner_bc = bc.getInnerBoundaryCol(i);
     // swap alphas, boundary conditions and sx/sy for column-wise calculation
-    A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy);
+    A = createCoeffMatrix(alphaY, bcTop, bcBottom, inner_bc, rowMax, i, sy);
     b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom,
-                             bcLeft, bcRight, rowMax, i, sy, sx);
+                             bcLeft, bcRight, inner_bc, rowMax, i, sy, sx);
 
     row_t1 = solverFunc(A, b);
 

test/testSimulation.cpp

@@ -151,3 +151,36 @@ TEST_CASE("Closed Boundaries - no change expected") {
   CHECK(checkSimilarityV2(input_values, grid.getConcentrations(), 1E-12) ==
         true);
 }
+
+TEST_CASE("Constant inner cell - 'absorbing' concentration") {
+  constexpr std::uint32_t nrows = 5;
+  constexpr std::uint32_t ncols = 5;
+
+  tug::Grid64 grid(nrows, ncols);
+
+  auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0);
+  auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
+  auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
+
+  grid.setConcentrations(concentrations);
+  grid.setAlpha(alphax, alphay);
+
+  tug::Boundary bc(grid);
+  // inner
+  bc.setInnerBoundary(2, 2, 0);
+
+  tug::Simulation<double> sim(grid, bc);
+  sim.setTimestep(1);
+  sim.setIterations(1);
+
+  MatrixXd input_values(concentrations);
+  sim.run();
+
+  CHECK(grid.getConcentrations().sum() < input_values.sum());
+
+  const bool greater_one = (grid.getConcentrations().array() > 1.0).any();
+  CHECK(greater_one == false);
+
+  const bool less_zero = (grid.getConcentrations().array() < 0.0).any();
+  CHECK(less_zero == false);
+}