Merge branch 'spatial_fix' into 'main'

Fix spatial discretization on outer cells of inlet

See merge request mluebke/diffusion!13

.gitlab-ci.yml

@@ -19,7 +19,7 @@ build:
         expire_in: 100s
     script:
         - mkdir build && cd build
-        - cmake ..
+        - cmake -DCMAKE_BUILD_TYPE=Debug ..
         - make
 
 run_1D:
@@ -41,7 +41,7 @@ lint:
     script:
         - mkdir lint && cd lint
         - cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*,readability-*, modernize-*" ..
-        - make
+        - make diffusion
 
 memcheck_1D:
     stage: dynamic_analyze

app/CMakeLists.txt

@@ -6,3 +6,6 @@ target_link_libraries(2D PUBLIC diffusion)
 
 add_executable(Comp2D main_2D_mdl.cpp)
 target_link_libraries(Comp2D PUBLIC diffusion)
+
+add_executable(Const2D main_2D_const.cpp)
+target_link_libraries(Const2D PUBLIC diffusion)

app/main_1D.cpp

@@ -20,7 +20,7 @@ int main(int argc, char *argv[]) {
   // create input + diffusion coefficients for each grid cell
   std::vector<double> alpha(n, 1 * pow(10, -1));
   std::vector<double> field(n, 1 * std::pow(10, -6));
-  std::vector<boundary_condition> bc(n, {0,0});
+  std::vector<boundary_condition> bc(n + 2, {0, 0});
 
   // create instance of diffusion module
   BTCSDiffusion diffu(dim);

app/main_2D.cpp

@@ -1,7 +1,7 @@
-#include <diffusion/BTCSDiffusion.hpp>
-#include <diffusion/BoundaryCondition.hpp>
 #include <algorithm> // for copy, max
 #include <cmath>
+#include <diffusion/BTCSDiffusion.hpp>
+#include <diffusion/BoundaryCondition.hpp>
 #include <iomanip>
 #include <iostream> // for std
 #include <vector>   // for vector
@@ -20,7 +20,7 @@ int main(int argc, char *argv[]) {
   // create input + diffusion coefficients for each grid cell
   std::vector<double> alpha(n * m, 1 * pow(10, -1));
   std::vector<double> field(n * m, 1 * std::pow(10, -6));
-  std::vector<boundary_condition> bc(n*m, {0,0});
+  std::vector<boundary_condition> bc((n + 2) * (m + 2), {0, 0});
 
   // create instance of diffusion module
   BTCSDiffusion diffu(dim);

app/main_2D_const.cpp

@@ -0,0 +1,57 @@
+#include <algorithm> // for copy, max
+#include <cmath>
+#include <diffusion/BTCSDiffusion.hpp>
+#include <diffusion/BoundaryCondition.hpp>
+#include <iomanip>
+#include <iostream> // for std
+#include <vector>   // for vector
+using namespace std;
+
+using namespace Diffusion;
+
+int main(int argc, char *argv[]) {
+
+  // dimension of grid
+  int dim = 2;
+
+  int n = 5;
+  int m = 5;
+
+  // create input + diffusion coefficients for each grid cell
+  std::vector<double> alpha(n * m, 1 * pow(10, -1));
+  std::vector<double> field(n * m, 1 * std::pow(10, -6));
+  std::vector<boundary_condition> bc((n + 2) * (m + 2), {0, 0});
+
+  // create instance of diffusion module
+  BTCSDiffusion diffu(dim);
+
+  diffu.setXDimensions(1, n);
+  diffu.setYDimensions(1, m);
+
+  for (int i = 1; i <= n; i++) {
+    bc[(n + 2) * i] = {Diffusion::BC_CONSTANT, 5 * std::pow(10, -6)};
+  }
+  // set timestep for simulation to 1 second
+  diffu.setTimestep(1.);
+
+  cout << setprecision(12);
+
+  for (int t = 0; t < 10; t++) {
+    diffu.simulate(field.data(), alpha.data(), bc.data());
+
+    cout << "Iteration: " << t << "\n\n";
+
+    // iterate through rows
+    for (int i = 0; i < m; i++) {
+      // iterate through columns
+      for (int j = 0; j < n; j++) {
+        cout << left << std::setw(20) << field[i * n + j];
+      }
+      cout << "\n";
+    }
+
+    cout << "\n" << endl;
+  }
+
+  return 0;
+}

app/main_2D_mdl.cpp

@@ -1,7 +1,7 @@
-#include <diffusion/BTCSDiffusion.hpp>
-#include <diffusion/BoundaryCondition.hpp>
 #include <algorithm> // for copy, max
 #include <cmath>
+#include <diffusion/BTCSDiffusion.hpp>
+#include <diffusion/BoundaryCondition.hpp>
 #include <iomanip>
 #include <iostream> // for std
 #include <vector>   // for vector
@@ -20,7 +20,7 @@ int main(int argc, char *argv[]) {
   // create input + diffusion coefficients for each grid cell
   std::vector<double> alpha(n * m, 1 * pow(10, -1));
   std::vector<double> field(n * m, 0.);
-  std::vector<boundary_condition> bc(n*m, {0,0});
+  std::vector<boundary_condition> bc((n + 2) * (m + 2), {0, 0});
 
   field[125500] = 1;
 

doc/ADI_scheme.org

@@ -1,9 +1,314 @@
-#+TITLE: Adi 2D Scheme
+#+TITLE: Numerical solution of diffusion equation in 2D with ADI Scheme
 #+LaTeX_CLASS_OPTIONS: [a4paper,10pt]
 #+LATEX_HEADER: \usepackage{fullpage}
+#+LATEX_HEADER: \usepackage{amsmath, systeme}
 #+OPTIONS: toc:nil
 
-* Input
+
+* Diffusion in 1D
+
+** Finite differences with nodes as cells' centres
+
+The 1D diffusion equation is:
+
+#+NAME: eqn:1
+\begin{align}
+\frac{\partial C }{\partial t} & = \frac{\partial}{\partial x} \left(\alpha \frac{\partial C }{\partial x} \right) \nonumber \\
+   & = \alpha \frac{\partial^2 C}{\partial x^2}
+\end{align}
+
+We aim at numerically solving [[eqn:1]] on a spatial grid such as:
+
+[[./grid_pqc.pdf]]
+
+The left boundary is defined on $x=0$ while the center of the first
+cell - which are the points constituting the finite difference nodes -
+is in $x=dx/2$, with $dx=L/n$.
+
+
+** The explicit FTCS scheme (as in PHREEQC)
+
+We start by discretizing [[eqn:1]] following an explicit Euler scheme and
+specifically a Forward Time, Centered Space finite difference. 
+
+For each cell index $i \in 1, \dots, n-1$ and assuming constant
+$\alpha$, we can write:
+
+#+NAME: eqn:2
+\begin{equation}\displaystyle
+   \frac{C_i^{t+1} -C_i^{t}}{\Delta t} = \alpha\frac{\frac{C^t_{i+1}-C^t_{i}}{\Delta x}-\frac{C^t_{i}-C^t_{i-1}}{\Delta x}}{\Delta x}
+\end{equation}
+
+In practice, we evaluate the first derivatives of $C$ w.r.t. $x$ on
+the boundaries of each cell (i.e., $(C_{i+1}-C_i)/\Delta x$ on the
+right boundary of the i-th cell and $(C_{i}-C_{i-1})/\Delta x$ on its
+left cell boundary) and then repeat the differentiation to get the
+second derivative of $C$ on the the cell centre $i$.
+
+This discretization works for all internal cells, but not for the
+domain boundaries ($i=0$ and $i=n$). To properly treat them, we need
+to account for the discrepancy in the discretization.
+
+For the first (left) cell, whose center is at $x=dx/2$, we can
+evaluate the left gradient with the left boundary using such distance,
+calling $l$ the numerical value of a constant boundary condition:
+
+#+NAME: eqn:3
+\begin{equation}\displaystyle
+\frac{C_0^{t+1} -C_0^{t}}{\Delta t} = \alpha\frac{\frac{C^t_{1}-C^t_{0}}{\Delta x}-
+\frac{C^t_{0}-l}{\frac{\Delta x}{2}}}{\Delta x}
+\end{equation}
+
+This expression, once developed, yields:
+
+#+NAME: eqn:4
+\begin{align}\displaystyle
+C_0^{t+1} & =  C_0^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^t_{1}-C^t_{0}- 2 C^t_{0}+2l \right) \nonumber \\
+          & =  C_0^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^t_{1}- 3 C^t_{0} +2l \right)
+\end{align}
+
+
+In case of constant right boundary, the finite difference of point
+$C_n$ - calling $r$ the right boundary value - is:
+
+#+NAME: eqn:5
+\begin{equation}\displaystyle
+\frac{C_n^{t+1} -C_n^t}{\Delta t} = \alpha\frac{\frac{r - C^t_{n}}{\frac{\Delta x}{2}}-
+\frac{C^t_{n}-C^t_{n-1}}{\Delta x}}{\Delta x}
+\end{equation}
+
+Which, developed, gives
+#+NAME: eqn:6
+\begin{align}\displaystyle
+C_n^{t+1} & =  C_n^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^t_{n} -C^t_{n} + C^t_{n-1} \right) \nonumber \\
+          & =  C_n^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^t_{n} + C^t_{n-1} \right)
+\end{align}
+
+If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:5]]
+becomes zero and we are left with:
+
+
+#+NAME: eqn:7
+\begin{equation}\displaystyle
+C_n^{t+1} = C_n^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^t_{n-1} - C^t_n)
+\end{equation}
+
+
+
+A similar treatment can be applied to the BTCS implicit scheme.
+
+** Implicit BTCS scheme
+
+First, we define the Backward Time difference:
+
+\begin{equation}
+    \frac{\partial C^{t+1} }{\partial t} = \frac{C^{t+1}_i - C^{t}_i}{\Delta t}
+\end{equation}
+
+Second the spatial derivative approximation, evaluated at time level $t+1$:
+
+#+NAME: eqn:secondderivative
+\begin{equation}
+    \frac{\partial^2 C^{t+1} }{\partial x^2} = \frac{\frac{C^{t+1}_{i+1}-C^{t+1}_{i}}{\Delta x}-\frac{C^{t+1}_{i}-C^{t+1}_{i-1}}{\Delta x}}{\Delta x}
+\end{equation}
+
+Taking the 1D diffusion equation from [[eqn:1]] and substituting each term by the
+equations given above leads to the following equation:
+
+
+# \begin{equation}\displaystyle
+#    \frac{C_i^{j+1} -C_i^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{i+1}-C^{j+1}_{i}}{\Delta x}-\frac{C^{j+1}_{i}-C^{j+1}_{i-1}}{\Delta x}}{\Delta x}
+# \end{equation}
+
+# Since we are not able to solve this system w.r.t unknown values in $C^{j-1}$ we
+# are shifting each j by 1 to $j \to (j+1)$ and $(j-1) \to j$ which leads to:
+
+#+NAME: eqn:1DBTCS
+\begin{align}\displaystyle
+\frac{C_i^{t+1} - C_i^{t}}{\Delta t}    & = \alpha\frac{\frac{C^{t+1}_{i+1}-C^{t+1}_{i}}{\Delta x}-\frac{C^{t+1}_{i}-C^{t+1}_{i-1}}{\Delta x}}{\Delta x} \nonumber \\
+                                        & = \alpha\frac{C^{t+1}_{i-1} - 2C^{t+1}_{i} + C^{t+1}_{i+1}}{\Delta x^2}
+\end{align}
+
+This only applies to inlet cells with no ghost node as neighbor. For the left
+cell with its center at $\frac{dx}{2}$ and the constant concentration on the
+left ghost node called $l$ the equation goes as followed:
+
+\begin{equation}\displaystyle
+\frac{C_0^{t+1} -C_0^{t}}{\Delta t} = \alpha\frac{\frac{C^{t+1}_{1}-C^{t+1}_{0}}{\Delta x}-
+\frac{C^{t+1}_{0}-l}{\frac{\Delta x}{2}}}{\Delta x}
+\end{equation}
+
+This expression, once developed, yields:
+
+\begin{align}\displaystyle
+C_0^{t+1} & =  C_0^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{t+1}_{1}-C^{t+1}_{0}- 2 C^{t+1}_{0}+2l \right) \nonumber \\
+          & =  C_0^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{t+1}_{1}- 3 C^{t+1}_{0} +2l \right)
+\end{align}
+
+Now we define variable $s_x$ as followed:
+
+\begin{equation}
+    s_x = \frac{\alpha \cdot \Delta t}{\Delta x^2}
+\end{equation}
+
+Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model:
+
+\begin{equation}\displaystyle
+    -C^t_0 = (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{t+1}_0 + s_x \cdot C^{t+1}_1
+\end{equation}
+
+The right boundary follows the same scheme. We now want to show the equation for the rightmost inlet cell $C_n$ with right boundary value $r$:
+
+\begin{equation}\displaystyle
+\frac{C_n^{t+1} -C_n^{t}}{\Delta t} = \alpha\frac{\frac{r-C^{t+1}_{n}}{\frac{\Delta x}{2}}-
+\frac{C^{t+1}_{n}-C^{t+1}_{n-1}}{\Delta x}}{\Delta x}
+\end{equation}
+
+This expression, once developed, yields:
+
+\begin{align}\displaystyle
+C_n^{t+1} & =  C_n^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 2C^{t+1}_{n} - C^{t+1}_{n} + C^{t+1}_{n-1} \right) \nonumber \\
+          & =  C_0^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 3C^{t+1}_{n} + C^{t+1}_{n-1} \right)
+\end{align}
+
+Now rearrange terms and substituting with $s_x$ leads to:
+
+\begin{equation}\displaystyle
+    -C^t_n = s_x \cdot C^{t+1}_{n-1} + (-1 - 3s_x) \cdot C^{t+1}_n + (2s_x) \cdot r
+\end{equation}
+
+*TODO*
+- Tridiagonal matrix filling
+
+#+LATEX: \clearpage
+
+* Diffusion in 2D: the Alternating Direction Implicit scheme
+
+
+In 2D, the diffusion equation (in absence of source terms) and
+assuming homogeneous but anisotropic diffusion coefficient
+$\alpha=(\alpha_x,\alpha_y)$ becomes:
+
+#+NAME: eqn:2d
+\begin{equation}
+\displaystyle  \frac{\partial C}{\partial t} = \alpha_x \frac{\partial^2 C}{\partial x^2} + \alpha_y\frac{\partial^2 C}{\partial y^2}
+\end{equation}
+
+** 2D ADI using BTCS scheme
+
+The Alternating Direction Implicit method consists in splitting the
+integration of eq. [[eqn:2d]] in two half-steps, each of which represents
+implicitly the derivatives in one direction, and explicitly in the
+other. Therefore we make use of second derivative operator defined in
+equation [[eqn:secondderivative]] in both $x$ and $y$ direction for half
+the time step $\Delta t$.
+
+Denoting $i$ the grid cell index along $x$ direction, $j$ the index in
+$y$ direction and $t$ as the time level, the spatially centered second
+derivatives can be written as:
+
+\begin{align}\displaystyle
+\frac{\partial^2 C^t_{i,j}}{\partial x^2} &= \frac{C^{t}_{i-1,j} - 2C^{t}_{i,j} + C^{t}_{i+1,j}}{\Delta x^2} \\
+\frac{\partial^2 C^t_{i,j}}{\partial y^2} &= \frac{C^{t}_{i,j-1} - 2C^{t}_{i,j} + C^{t}_{i,j+1}}{\Delta y^2}
+\end{align}
+
+The ADI scheme is formally defined by the equations:
+
+#+NAME: eqn:genADI
+\begin{equation}
+\systeme{ 
+  \displaystyle  \frac{C^{t+1/2}_{i,j}-C^t_{i,j}}{\Delta t/2}     = \displaystyle \alpha_x \frac{\partial^2 C^{t+1/2}_{i,j}}{\partial x^2} + \alpha_y \frac{\partial^2 C^{t}_{i,j}}{\partial y^2},
+  \displaystyle  \frac{C^{t+1}_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = \displaystyle \alpha_x \frac{\partial^2 C^{t+1/2}_{i,j}}{\partial x^2} + \alpha_y \frac{\partial^2 C^{t+1}_{i,j}}{\partial y^2}
+}
+\end{equation}
+
+\noindent The first of equations [[eqn:genADI]], which writes implicitly
+the spatial derivatives in $x$ direction, after bringing the $\Delta t
+/ 2$ terms on the right hand side and substituting $s_x=\frac{\alpha_x
+\cdot \Delta t}{2 \Delta x^2}$ and $s_y=\frac{\alpha_y\cdot\Delta
+t}{2\Delta y^2}$ reads:
+
+\begin{equation}\displaystyle
+C^{t+1/2}_{i,j}-C^t_{i,j} = s_x (C^{t+1/2}_{i-1,j} - 2C^{t+1/2}_{i,j} + C^{t+1/2}_{i+1,j}) + s_y (C^{t}_{i,j-1} - 2C^{t}_{i,j} + C^{t}_{i,j+1})
+\end{equation}
+
+\noindent Separating the known terms (at time level $t$) on the left
+hand side and the implicit terms (at time level $t+1/2$) on the right
+hand side, we get:
+
+#+NAME: eqn:sweepX
+\begin{equation}\displaystyle
+-C^t_{i,j} - s_y (C^{t}_{i,j-1} - 2C^{t}_{i,j} + C^{t}_{i,j+1}) = - C^{t+1/2}_{i,j} + s_x (C^{t+1/2}_{i-1,j} - 2C^{t+1/2}_{i,j} + C^{t+1/2}_{i+1,j}) 
+\end{equation}
+
+\noindent Equation [[eqn:sweepX]] can be solved with a BTCS scheme since
+it corresponds to a tridiagonal system of equations, and resolves the
+$C^{t+1/2}$ at each inner grid cell.
+
+The second of equations [[eqn:genADI]] can be treated the same way and
+yields:
+
+#+NAME: eqn:sweepY
+\begin{equation}\displaystyle
+-C^{t + 1/2}_{i,j} - s_x (C^{t + 1/2}_{i-1,j} - 2C^{t + 1/2}_{i,j} + C^{t + 1/2}_{i+1,j}) = - C^{t+1}_{i,j} + s_y (C^{t+1}_{i,j-1} - 2C^{t+1}_{i,j} + C^{t+1}_{i,j+1}) 
+\end{equation}
+
+This scheme only applies to inner cells, or else $\forall i,j \in [1,
+n-1] \times [1, n-1]$. Following an analogous treatment as for the 1D
+case, and noting $l_x$ and $l_y$ the constant left boundary values and
+$r_x$ and $r_y$ the right ones for each direction $x$ and $y$, we can
+modify equations [[eqn:sweepX]] for $i=0, j \in [1, n-1]$
+
+#+NAME: eqn:boundXleft
+\begin{equation}\displaystyle
+-C^t_{0,j} - s_y (C^{t}_{0,j-1} - 2C^{t}_{0,j} + C^{t}_{0,j+1}) = - C^{t+1/2}_{0,j} + s_x (C^{t+1/2}_{1,j} - 3C^{t+1/2}_{0,j} + 2 l_x) 
+\end{equation}
+
+\noindent Similarly for  $i=n, j \in [1, n-1]$:
+#+NAME: eqn:boundXright
+\begin{equation}\displaystyle
+-C^t_{n,j} - s_y (C^{t}_{n,j-1} - 2C^{t}_{n,j} + C^{t}_{n,j+1}) = - C^{t+1/2}_{n,j} + s_x (C^{t+1/2}_{n-1,j} - 3C^{t+1/2}_{n,j} + 2 r_x) 
+\end{equation}
+
+\noindent For $i=j=0$:
+#+NAME: eqn:bound00
+\begin{equation}\displaystyle
+-C^t_{0,0} - s_y (C^{t}_{0,1} - 3C^{t}_{0,0} + 2l_y) = - C^{t+1/2}_{0,0} + s_x (C^{t+1/2}_{1,0} - 3C^{t+1/2}_{0,0} + 2 l_x) 
+\end{equation}
+
+Analogous expressions are readily derived for all possible
+combinations of $i,j \in 0\times n$. In practice, wherever an index
+$i$ or $j$ is $0$ or $n$, the centered spatial derivatives in $x$ or
+$y$ directions must be substituted in relevant parts of the sweeping
+equations \textbf{in both the implicit or the explicit sides} of
+equations [[eqn:sweepX]] and [[eqn:sweepY]] by a term
+
+#+NAME: eqn:bound00
+\begin{equation}\displaystyle
+ s(C_{forw} - 3C + 2 bc) 
+\end{equation}
+\noindent where $bc$ is the boundary condition in the given direction,
+$s$ is either $s_x$ or $s_y$, and $C_{forw}$ indicates the contiguous
+cell opposite to the boundary. Alternatively, noting the second
+derivative operator as $\partial_{dir}^2$, we can write in compact
+form:
+
+\begin{equation}
+\systeme{
+  \displaystyle  \partial_x^2 C_{0,j} = 2l_x - 3C_{0,j} + C_{1,j} ,
+  \displaystyle  \partial_x^2 C_{n,j} = 2r_x - 3C_{n,j} + C_{n-1,j} ,
+  \displaystyle  \partial_y^2 C_{i,0} = 2l_y - 3C_{i,0} + C_{i,1} ,
+  \displaystyle  \partial_y^2 C_{i,n} = 2r_y - 3C_{i,n} + C_{i,n-1}
+}
+\end{equation}
+
+
+
+#+LATEX: \clearpage
+
+* Old stuff
+
+** Input
 
 - =c= $\rightarrow c$
   - containing current concentrations at each grid cell for species
@@ -18,7 +323,7 @@
   - size: $N \times M$
   - row-major
 
-* Internals
+** Internals
 
 - =A_matrix= $\rightarrow A$
   - coefficient matrix for linear equation system implemented as sparse matrix
@@ -34,7 +339,7 @@
   - size: $(N+2) \cdot M$
   - column-major (not relevant)
 
-* Calculation for $\frac{1}{2}$ timestep
+** Calculation for $\frac{1}{2}$ timestep
 
 ** Symbolic addressing of grid cells
 [[./grid.png]]
@@ -77,9 +382,6 @@ s_x(i,j) & \text{else}
 - Adressing would look like this: $(i,j) = b(i \cdot (N+2) + j)$
   - $\rightarrow$ for simplicity we will write $b(i,j)$
 
-
-
-
 *** Rules
 
 **** Ghost nodes
@@ -96,11 +398,11 @@ c(i,N-1) & \text{else}
 
 *** Inlet
 
-$p(i,j) = \frac{\Delta t}{2}\alpha(i,j)\frac{c(i-1,j) - 2\cdot c(i,j) + c(i+1,j)}{\Delta x^2}$[fn:1]
+$p(i,j) = \frac{\Delta t}{2}\alpha(i,j)\frac{c(i-1,j) - 2\cdot c(i,j) + c(i+1,j)}{\Delta x^2}$
+
+\noindent $p$ is called =t0_c= inside code
 
 $b(i,j) = \begin{cases}
 bc(i,j).\text{value} & \text{if } bc(i,N-1) = \text{constant} \\
 -c(i,j)-p(i,j) & \text{else}
 \end{cases}$
-
-[fn:1] $p$ is called =t0_c= inside code

include/diffusion/BTCSDiffusion.hpp

@@ -135,16 +135,17 @@ private:
                  const BCMatrixRowMajor &bc, double time_step, double dx)
       -> DMatrixRowMajor;
 
-  void fillMatrixFromRow(Eigen::SparseMatrix<double> &A_matrix,
+  static void fillMatrixFromRow(Eigen::SparseMatrix<double> &A_matrix,
                                 const DVectorRowMajor &alpha,
                                 const BCVectorRowMajor &bc, int size, double dx,
                                 double time_step);
 
-  void fillVectorFromRow(Eigen::VectorXd &b_vector, const DVectorRowMajor &c,
+  static void fillVectorFromRow(Eigen::VectorXd &b_vector,
+                                const DVectorRowMajor &c,
                                 const DVectorRowMajor &alpha,
                                 const BCVectorRowMajor &bc,
-                         const DVectorRowMajor &t0_c, int size, double dx,
+                                const DVectorRowMajor &t0_c, int size,
-                         double time_step);
+                                double dx, double time_step);
   void simulate3D(std::vector<double> &c);
 
   inline static auto getBCFromFlux(Diffusion::boundary_condition bc,

src/BTCSDiffusion.cpp

@@ -28,8 +28,6 @@
 
 constexpr int BTCS_MAX_DEP_PER_CELL = 3;
 constexpr int BTCS_2D_DT_SIZE = 2;
-constexpr double center_eq(double sx) { return -1. - 2. * sx; }
-
 Diffusion::BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) {
 
   grid_cells.resize(dim, 1);
@@ -104,7 +102,7 @@ void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c,
   b_vector.resize(size + 2);
   x_vector.resize(size + 2);
 
-  fillMatrixFromRow(A_matrix, alpha.row(0), bc.row(0), size, dx, time_step);
+  fillMatrixFromRow(A_matrix, alpha, bc, size, dx, time_step);
   fillVectorFromRow(b_vector, c, alpha, bc, Eigen::VectorXd::Constant(size, 0),
                     size, dx, time_step);
 
@@ -143,17 +141,16 @@ void Diffusion::BTCSDiffusion::simulate2D(
   int n_rows = this->grid_cells[1];
   int n_cols = this->grid_cells[0];
   double dx = this->deltas[0];
-  DMatrixRowMajor t0_c;
 
   double local_dt = this->time_step / BTCS_2D_DT_SIZE;
 
-  t0_c = calc_t0_c(c, alpha, bc, local_dt, dx);
+  DMatrixRowMajor t0_c = calc_t0_c(c, alpha, bc, local_dt, dx);
 
 #pragma omp parallel for schedule(dynamic)
   for (int i = 0; i < n_rows; i++) {
     DVectorRowMajor input_field = c.row(i);
-    simulate_base(input_field, bc.row(i), alpha.row(i), dx, local_dt, n_cols,
+    simulate_base(input_field, bc.row(i + 1), alpha.row(i), dx, local_dt,
-                  t0_c.row(i));
+                  n_cols, t0_c.row(i));
     c.row(i) << input_field;
   }
 
@@ -165,8 +162,8 @@ void Diffusion::BTCSDiffusion::simulate2D(
 #pragma omp parallel for schedule(dynamic)
   for (int i = 0; i < n_cols; i++) {
     DVectorRowMajor input_field = c.col(i);
-    simulate_base(input_field, bc.col(i), alpha.col(i), dx, local_dt, n_rows,
+    simulate_base(input_field, bc.col(i + 1), alpha.col(i), dx, local_dt,
-                  t0_c.row(i));
+                  n_rows, t0_c.row(i));
     c.col(i) << input_field.transpose();
   }
 }
@@ -182,16 +179,22 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c,
 
   DMatrixRowMajor t0_c(n_rows, n_cols);
 
-  std::array<double, 3> y_values;
+  std::array<double, 3> y_values{};
 
   // first, iterate over first row
   for (int j = 0; j < n_cols; j++) {
-    y_values[0] = getBCFromFlux(bc(0, j), c(0, j), alpha(0, j));
+    boundary_condition tmp_bc = bc(0, j + 1);
+
+    if (tmp_bc.type == Diffusion::BC_CLOSED){
+      continue;
+    }
+
+    y_values[0] = getBCFromFlux(tmp_bc, c(0, j), alpha(0, j));
     y_values[1] = c(0, j);
     y_values[2] = c(1, j);
 
     t0_c(0, j) = time_step * alpha(0, j) *
-                 (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx);
+                 (2 * y_values[0] - 3 * y_values[1] + y_values[2]) / (dx * dx);
   }
 
 // then iterate over inlet
@@ -212,12 +215,19 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c,
 
   // and finally over last row
   for (int j = 0; j < n_cols; j++) {
+    boundary_condition tmp_bc = bc(end + 1, j + 1);
+
+    if (tmp_bc.type == Diffusion::BC_CLOSED) {
+      continue;
+    }
+
     y_values[0] = c(end - 1, j);
     y_values[1] = c(end, j);
-    y_values[2] = getBCFromFlux(bc(end, j), c(end, j), alpha(end, j));
+    y_values[2] = getBCFromFlux(tmp_bc, c(end, j), alpha(end, j));
 
     t0_c(end, j) = time_step * alpha(end, j) *
-                   (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx);
+                   (y_values[0] - 3 * y_values[1] + 2 * y_values[2]) /
+                   (dx * dx);
   }
 
   return t0_c;
@@ -227,8 +237,8 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow(
     Eigen::SparseMatrix<double> &A_matrix, const DVectorRowMajor &alpha,
     const BCVectorRowMajor &bc, int size, double dx, double time_step) {
 
-  Diffusion::boundary_condition left = bc[0];
+  Diffusion::boundary_condition left = bc[1];
-  Diffusion::boundary_condition right = bc[size - 1];
+  Diffusion::boundary_condition right = bc[size];
 
   bool left_constant = (left.type == Diffusion::BC_CONSTANT);
   bool right_constant = (right.type == Diffusion::BC_CONSTANT);
@@ -239,27 +249,36 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow(
 
   if (left_constant) {
     A_matrix.insert(1, 1) = 1;
+  } else {
+    double sx = (alpha[0] * time_step) / (dx * dx);
+    A_matrix.insert(1, 1) = -1. - 3. * sx;
+    A_matrix.insert(1, 0) = 2. * sx;
+    A_matrix.insert(1, 2) = sx;
   }
 
-  A_matrix.insert(A_size - 1, A_size - 1) = 1;
+  for (int j = 2, k = j - 1; k < size - 1; j++, k++) {
-
-  if (right_constant) {
-    A_matrix.insert(A_size - 2, A_size - 2) = 1;
-  }
-
-  for (int j = 1 + (int)left_constant, k = j - 1;
-       k < size - (int)right_constant; j++, k++) {
     double sx = (alpha[k] * time_step) / (dx * dx);
 
-    if (bc[k].type == Diffusion::BC_CONSTANT) {
+    if (bc[k + 1].type == Diffusion::BC_CONSTANT) {
       A_matrix.insert(j, j) = 1;
       continue;
     }
 
-    A_matrix.insert(j, j) = center_eq(sx);
+    A_matrix.insert(j, j) = -1. - 2. * sx;
     A_matrix.insert(j, (j - 1)) = sx;
     A_matrix.insert(j, (j + 1)) = sx;
   }
+
+  if (right_constant) {
+    A_matrix.insert(A_size - 2, A_size - 2) = 1;
+  } else {
+    double sx = (alpha[size - 1] * time_step) / (dx * dx);
+    A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx;
+    A_matrix.insert(A_size - 2, A_size - 3) = sx;
+    A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx;
+  }
+
+  A_matrix.insert(A_size - 1, A_size - 1) = 1;
 }
 
 void Diffusion::BTCSDiffusion::fillVectorFromRow(
@@ -268,7 +287,7 @@ void Diffusion::BTCSDiffusion::fillVectorFromRow(
     const DVectorRowMajor &t0_c, int size, double dx, double time_step) {
 
   Diffusion::boundary_condition left = bc[0];
-  Diffusion::boundary_condition right = bc[size - 1];
+  Diffusion::boundary_condition right = bc[size + 1];
 
   bool left_constant = (left.type == Diffusion::BC_CONSTANT);
   bool right_constant = (right.type == Diffusion::BC_CONSTANT);
@@ -276,7 +295,7 @@ void Diffusion::BTCSDiffusion::fillVectorFromRow(
   int b_size = b_vector.size();
 
   for (int j = 0; j < size; j++) {
-    boundary_condition tmp_bc = bc[j];
+    boundary_condition tmp_bc = bc[j + 1];
 
     if (tmp_bc.type == Diffusion::BC_CONSTANT) {
       b_vector[j + 1] = tmp_bc.value;
@@ -287,15 +306,14 @@ void Diffusion::BTCSDiffusion::fillVectorFromRow(
     b_vector[j + 1] = -c[j] - t0_c_j;
   }
 
-  if (!left_constant) {
   // this is not correct currently.We will fix this when we are able to define
   // FLUX boundary conditions
-    b_vector[0] = getBCFromFlux(left, c[0], alpha[0]);
+  b_vector[0] =
-  }
+      (left_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
 
-  if (!right_constant) {
+  b_vector[b_size - 1] =
-    b_vector[b_size - 1] = getBCFromFlux(right, c[size - 1], alpha[size - 1]);
+      (right_constant ? right.value
-  }
+                      : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
 }
 
 void Diffusion::BTCSDiffusion::setTimestep(double time_step) {
@@ -312,7 +330,7 @@ auto Diffusion::BTCSDiffusion::simulate(double *c, double *alpha,
   if (this->grid_dim == 1) {
     Eigen::Map<DVectorRowMajor> c_in(c, this->grid_cells[0]);
     Eigen::Map<const DVectorRowMajor> alpha_in(alpha, this->grid_cells[0]);
-    Eigen::Map<const BCVectorRowMajor> bc_in(bc, this->grid_cells[0]);
+    Eigen::Map<const BCVectorRowMajor> bc_in(bc, this->grid_cells[0] + 2);
 
     simulate1D(c_in, alpha_in, bc_in);
   }
@@ -323,8 +341,8 @@ auto Diffusion::BTCSDiffusion::simulate(double *c, double *alpha,
     Eigen::Map<const DMatrixRowMajor> alpha_in(alpha, this->grid_cells[1],
                                                this->grid_cells[0]);
 
-    Eigen::Map<const BCMatrixRowMajor> bc_in(bc, this->grid_cells[1],
+    Eigen::Map<const BCMatrixRowMajor> bc_in(bc, this->grid_cells[1] + 2,
-                                             this->grid_cells[0]);
+                                             this->grid_cells[0] + 2);
 
     simulate2D(c_in, alpha_in, bc_in);
   }