From f42af515fb327117cdabe75a2df1612e52566508 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 10:06:38 +0200 Subject: [PATCH 01/13] Update applicable equation --- doc/ADI_scheme.org | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 8a3fd92..557fe34 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -163,11 +163,17 @@ C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) \end{align} -Which leads to the following term applicable to our approach: +Now we define variable $s_x$ as following: + +\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: #+NAME: eqn:5 \begin{equation}\displaystyle - -C^j_0} = \left[\frac{\alpha \cdot \Delta t}{\Delta x^2}\right] \cdot \left(C^{j+1}_1 + 2l \right) + \left[ -1 - 3 \cdot \frac{\alpha \cdot \Delta t}{\Delta x^2} \right] \cdot C^{j+1}_0 + -C^j_0} = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} In case of constant right boundary, the finite difference of point From f31937f4e0cff661611cb318b829e8d99f3a08b6 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 10:09:43 +0200 Subject: [PATCH 02/13] Fix broken reference --- doc/ADI_scheme.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 557fe34..3cabbf1 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -192,7 +192,7 @@ C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{n-1} \right) \end{align} -If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:5]] +If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:6]] becomes zero and we are left with: From 29d3d81ff4e41e649fe05ec81402d5a848326c36 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 11:27:30 +0200 Subject: [PATCH 03/13] Apply bc to ghost node --- app/main_1D.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/app/main_1D.cpp b/app/main_1D.cpp index e908039..2fb7c21 100644 --- a/app/main_1D.cpp +++ b/app/main_1D.cpp @@ -28,7 +28,7 @@ int main(int argc, char *argv[]) { diffu.setXDimensions(1, n); // set the boundary condition for the left ghost cell to dirichlet - bc[1] = {Diffusion::BC_CONSTANT, 5 * std::pow(10, -6)}; + bc[0] = {Diffusion::BC_CONSTANT, 5 * std::pow(10, -6)}; // diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT, // 5. * std::pow(10, -6)); From 08da0b47f44159f208319b81d9b995e79eed25e9 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 11:30:59 +0200 Subject: [PATCH 04/13] Revert to 'age' state and append with implicit BTCS --- doc/ADI_scheme.org | 96 ++++++++++++++++++++++++++++++---------------- 1 file changed, 64 insertions(+), 32 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 3cabbf1..09158ed 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -132,7 +132,7 @@ and assuming constant $\alpha$: #+NAME: eqn:2 \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} + \frac{C_i^{j+1} -C_i^{j}}{\Delta t} = \alpha\frac{\frac{C^j_{i+1}-C^j_{i}}{\Delta x}-\frac{C^j_{i}-C^j_{i-1}}{\Delta x}}{\Delta x} \end{equation} In practice, we evaluate the first derivatives of $C$ w.r.t. $x$ on @@ -150,6 +150,69 @@ 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^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^j_{1}-C^j_{0}}{\Delta x}- +\frac{C^j_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} +\end{equation} + +This expression, once developed, yields: + +#+NAME: eqn:4 +\begin{align}\displaystyle +C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^j_{1}-C^j_{0}- 2 C^j_{0}+2l \right) \nonumber \\ + & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^j_{1}- 3 C^j_{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^{j+1} -C_n^j}{\Delta t} = \alpha\frac{\frac{r - C^j_{n}}{\frac{\Delta x}{2}}- +\frac{C^j_{n}-C^j_{n-1}}{\Delta x}}{\Delta x} +\end{equation} + +Which, developed, gives +#+NAME: eqn:6 +\begin{align}\displaystyle +C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^j_{n} -C^j_{n} + C^j_{n-1} \right) \nonumber \\ + & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{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^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} - C^j_n) +\end{equation} + + + +A similar treatment can be applied to the BTCS implicit scheme. + +*** implicit BTCS + +\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} + +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 +boundaries. 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: + \begin{equation}\displaystyle \frac{C_0^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{1}-C^{j+1}_{0}}{\Delta x}- \frac{C^{j+1}_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} @@ -157,7 +220,6 @@ calling $l$ the numerical value of a constant boundary condition: This expression, once developed, yields: -#+NAME: eqn:4 \begin{align}\displaystyle C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}-C^{j+1}_{0}- 2 C^{j+1}_{0}+2l \right) \nonumber \\ & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) @@ -171,36 +233,6 @@ Now we define variable $s_x$ as following: Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model: -#+NAME: eqn:5 \begin{equation}\displaystyle -C^j_0} = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} - -In case of constant right boundary, the finite difference of point -$C_n$ - calling $r$ the right boundary value - is: - -#+NAME: eqn:6 -\begin{equation}\displaystyle -\frac{C_n^{j+1} -C_n^j}{\Delta t} = \alpha\frac{\frac{r - C^j_{n}}{\frac{\Delta x}{2}}- -\frac{C^j_{n}-C^j_{n-1}}{\Delta x}}{\Delta x} -\end{equation} - -Which, developed, gives -#+NAME: eqn:7 -\begin{align}\displaystyle -C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^j_{n} -C^j_{n} + C^j_{n-1} \right) \nonumber \\ - & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{n-1} \right) -\end{align} - -If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:6]] -becomes zero and we are left with: - - -#+NAME: eqn:8 -\begin{equation}\displaystyle -C_n^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} - C^j_n) -\end{equation} - - - -A similar treatment can be applied to the BTCS implicit scheme. From 0f646d478c63a780817eaded175a4d2a361677c9 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Thu, 21 Apr 2022 12:51:34 +0200 Subject: [PATCH 05/13] Fix unnecessary bracket --- doc/ADI_scheme.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 09158ed..d47f3cb 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -234,5 +234,5 @@ Now we define variable $s_x$ as following: Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model: \begin{equation}\displaystyle - -C^j_0} = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + -C^j_0 = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} From d54fe25cace10879e7846f67cedc20b3d8955db1 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Fri, 22 Apr 2022 09:49:53 +0200 Subject: [PATCH 06/13] Update documentation of implicit BTCS --- doc/ADI_scheme.org | 41 +++++++++++++++++++++++++++-------------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index d47f3cb..7a0b0c7 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -195,23 +195,36 @@ A similar treatment can be applied to the BTCS implicit scheme. *** implicit BTCS -\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} +First, we define the Backward time difference: + +\begin{equation} + \frac{\partial C }{\partial t} = \frac{C^j_i - C^{j-1}_i}{\Delta t} \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$. +Second the spatial derivative approximation: -This discretization works for all internal cells, but not for the -boundaries. To properly treat them, we need to account for the -discrepancy in the discretization. +\begin{equation} + \frac{\partial^2 C }{\partial t} = \frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{i-1}}{\Delta x}}{\Delta x} +\end{equation} -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: +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} -C_i^{j-1}}{\Delta t} = \alpha\frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{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: + +\begin{align}\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} \nonumber \\ + & = \alpha\frac{C^{j+1}_{i-1} - 2C^{j+1}_{i} + C^{j+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^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{1}-C^{j+1}_{0}}{\Delta x}- @@ -225,7 +238,7 @@ C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) \end{align} -Now we define variable $s_x$ as following: +Now we define variable $s_x$ as followed: \begin{equation} s_x = \frac{\alpha \cdot \Delta t}{\Delta x^2} From 8b4f1aae46164c707f52c1cc7e2c6732467c78bd Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Fri, 22 Apr 2022 13:55:50 +0200 Subject: [PATCH 07/13] MDL: some restructuring in ADI_scheme.org --- doc/ADI_scheme.org | 259 ++++++++++++++++++++++++--------------------- 1 file changed, 136 insertions(+), 123 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 7a0b0c7..a3497aa 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -1,114 +1,11 @@ -#+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} #+OPTIONS: toc:nil -* Input -- =c= $\rightarrow c$ - - containing current concentrations at each grid cell for species - - size: $N \times M$ - - row-major -- =alpha= $\rightarrow \alpha$ - - diffusion coefficient for both directions (x and y) - - size: $N \times M$ - - row-major -- =boundary_condition= $\rightarrow bc$ - - Defines closed or constant boundary condition for each grid cell - - size: $N \times M$ - - row-major - -* Internals - -- =A_matrix= $\rightarrow A$ - - coefficient matrix for linear equation system implemented as sparse matrix - - size: $((N+2)\cdot M) \times ((N+2)\cdot M)$ (including ghost zones in x direction) - - column-major (not relevant) - -- =b_vector= $\rightarrow b$ - - right hand side of the linear equation system - - size: $(N+2) \cdot M$ - - column-major (not relevant) -- =x_vector= $\rightarrow x$ - - solutions of the linear equation system - - size: $(N+2) \cdot M$ - - column-major (not relevant) - -* Calculation for $\frac{1}{2}$ timestep - -** Symbolic addressing of grid cells -[[./grid.png]] - -** Filling of matrix $A$ - -- row-wise iterating with $i$ over =c= and =\alpha= matrix respectively -- addressing each element of a row with $j$ -- matrix $A$ also containing $+2$ ghost nodes for each row of input matrix $\alpha$ - - $\rightarrow offset = N+2$ - - addressing each object $(i,j)$ in matrix $A$ with $(offset \cdot i + j, offset \cdot i + j)$ - -*** Rules - -$s_x(i,j) = \frac{\alpha(i,j)*\frac{t}{2}}{\Delta x^2}$ where $x$ defining the domain size in x direction. - -For the sake of simplicity we assume that each row of the $A$ matrix is addressed correctly with the given offset. - -**** Ghost nodes - -$A(i,-1) = 1$ - -$A(i,N) = 1$ - -**** Inlet - -$A(i,j) = \begin{cases} -1 & \text{if } bc(i,j) = \text{constant} \\ --1-2*s_x(i,j) & \text{else} -\end{cases}$ - -$A(i,j\pm 1) = \begin{cases} -0 & \text{if } bc(i,j) = \text{constant} \\ -s_x(i,j) & \text{else} -\end{cases}$ - -** Filling of vector $b$ - -- each elements assign a concrete value to the according value of the row of matrix $A$ -- 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 - -$b(i,-1) = \begin{cases} -0 & \text{if } bc(i,0) = \text{constant} \\ -c(i,0) & \text{else} -\end{cases}$ - -$b(i,N) = \begin{cases} -0 & \text{if } bc(i,N-1) = \text{constant} \\ -c(i,N-1) & \text{else} -\end{cases}$ - -*** 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] - -$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 - -** Finite differences with nodes as cells' centres - -*** The explicit FTCS scheme as in PHREEQC +* Finite differences with nodes as cells' centres The 1D diffusion equation is: @@ -118,17 +15,22 @@ The 1D diffusion equation is: & = \alpha \frac{\partial^2 C}{\partial x^2} \end{align} -We discretize it following a Forward Time, Centered Space finite -difference scheme where the nodes correspond to the centers of a grid -such as: +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 is in $x=dx/2$, with $dx=L/n$. +cell - which are the points constituting the finite difference nodes - +is in $x=dx/2$, with $dx=L/n$. -We discretize [[eqn:1]] as following, for each index i in 1, \dots, n-1 -and assuming constant $\alpha$: + +** 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 @@ -142,8 +44,8 @@ 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 -boundaries. To properly treat them, we need to account for the -discrepancy in the discretization. +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, @@ -193,29 +95,30 @@ C_n^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} A similar treatment can be applied to the BTCS implicit scheme. -*** implicit BTCS +** Implicit BTCS scheme -First, we define the Backward time difference: +First, we define the Backward Time difference: \begin{equation} - \frac{\partial C }{\partial t} = \frac{C^j_i - C^{j-1}_i}{\Delta t} + \frac{\partial C^{j+1} }{\partial t} = \frac{C^{j+1}_i - C^{j}_i}{\Delta t} \end{equation} -Second the spatial derivative approximation: +Second the spatial derivative approximation, evaluated at time level $j+1$: \begin{equation} - \frac{\partial^2 C }{\partial t} = \frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{i-1}}{\Delta x}}{\Delta x} + \frac{\partial^2 C^{j+1} }{\partial x^2} = \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} 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} -C_i^{j-1}}{\Delta t} = \alpha\frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{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: +# \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: \begin{align}\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} \nonumber \\ @@ -249,3 +152,113 @@ Substituting with the new variable $s_x$ and reordering of terms leads to the eq \begin{equation}\displaystyle -C^j_0 = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} + +*TODO* +- Right boundary +- Tridiagonal matrix filling + + + + +#+LATEX: \clearpage + +* Old stuff + +** Input + +- =c= $\rightarrow c$ + - containing current concentrations at each grid cell for species + - size: $N \times M$ + - row-major +- =alpha= $\rightarrow \alpha$ + - diffusion coefficient for both directions (x and y) + - size: $N \times M$ + - row-major +- =boundary_condition= $\rightarrow bc$ + - Defines closed or constant boundary condition for each grid cell + - size: $N \times M$ + - row-major + +** Internals + +- =A_matrix= $\rightarrow A$ + - coefficient matrix for linear equation system implemented as sparse matrix + - size: $((N+2)\cdot M) \times ((N+2)\cdot M)$ (including ghost zones in x direction) + - column-major (not relevant) + +- =b_vector= $\rightarrow b$ + - right hand side of the linear equation system + - size: $(N+2) \cdot M$ + - column-major (not relevant) +- =x_vector= $\rightarrow x$ + - solutions of the linear equation system + - size: $(N+2) \cdot M$ + - column-major (not relevant) + +** Calculation for $\frac{1}{2}$ timestep + +** Symbolic addressing of grid cells +[[./grid.png]] + +** Filling of matrix $A$ + +- row-wise iterating with $i$ over =c= and =\alpha= matrix respectively +- addressing each element of a row with $j$ +- matrix $A$ also containing $+2$ ghost nodes for each row of input matrix $\alpha$ + - $\rightarrow offset = N+2$ + - addressing each object $(i,j)$ in matrix $A$ with $(offset \cdot i + j, offset \cdot i + j)$ + +*** Rules + +$s_x(i,j) = \frac{\alpha(i,j)*\frac{t}{2}}{\Delta x^2}$ where $x$ defining the domain size in x direction. + +For the sake of simplicity we assume that each row of the $A$ matrix is addressed correctly with the given offset. + +**** Ghost nodes + +$A(i,-1) = 1$ + +$A(i,N) = 1$ + +**** Inlet + +$A(i,j) = \begin{cases} +1 & \text{if } bc(i,j) = \text{constant} \\ +-1-2*s_x(i,j) & \text{else} +\end{cases}$ + +$A(i,j\pm 1) = \begin{cases} +0 & \text{if } bc(i,j) = \text{constant} \\ +s_x(i,j) & \text{else} +\end{cases}$ + +** Filling of vector $b$ + +- each elements assign a concrete value to the according value of the row of matrix $A$ +- 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 + +$b(i,-1) = \begin{cases} +0 & \text{if } bc(i,0) = \text{constant} \\ +c(i,0) & \text{else} +\end{cases}$ + +$b(i,N) = \begin{cases} +0 & \text{if } bc(i,N-1) = \text{constant} \\ +c(i,N-1) & \text{else} +\end{cases}$ + +*** 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}$ + +\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}$ From aa91824eefe82eff854240d584df6b58e38c8216 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 25 Apr 2022 10:52:23 +0200 Subject: [PATCH 08/13] Reorder equation for better readability --- doc/ADI_scheme.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index a3497aa..d363194 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -150,7 +150,7 @@ Now we define variable $s_x$ as followed: Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model: \begin{equation}\displaystyle - -C^j_0 = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + -C^j_0 = (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + s_x \cdot C^{j+1}_1 \end{equation} *TODO* From e8b1c1daae4395dc29b36181cc8f1d8257a22ad7 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 25 Apr 2022 10:50:06 +0200 Subject: [PATCH 09/13] Added right boundary --- doc/ADI_scheme.org | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index d363194..391d3a6 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -153,8 +153,27 @@ Substituting with the new variable $s_x$ and reordering of terms leads to the eq -C^j_0 = (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + s_x \cdot C^{j+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^{j+1} -C_n^{j}}{\Delta t} = \alpha\frac{\frac{r-C^{j+1}_{n}}{\frac{\Delta x}{2}}- +\frac{C^{j+1}_{n}-C^{j+1}_{n-1}}{\Delta x}}{\Delta x} +\end{equation} + +This expression, once developed, yields: + +\begin{align}\displaystyle +C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 2C^{j+1}_{n} - C^{j+1}_{n} + C^{j+1}_{n-1} \right) \nonumber \\ + & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 3C^{j+1}_{n} + C^{j+1}_{n-1} \right) +\end{align} + +Now rearrange terms and substituting with $s_x$ leads to: + +\begin{equation}\displaystyle + -C^j_n = s_x \cdot C^{j+1}_{n-1} + (-1 - 3s_x) \cdot C^{j+1}_n + (2s_x) \cdot r +\end{equation} + *TODO* -- Right boundary - Tridiagonal matrix filling From ffefb6ab50339dde222e27ae4a44b47cbcc9f0e8 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 25 Apr 2022 12:31:27 +0200 Subject: [PATCH 10/13] Added 2D ADI scheme --- doc/ADI_scheme.org | 48 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 391d3a6..3f3e303 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -105,6 +105,7 @@ First, we define the Backward Time difference: Second the spatial derivative approximation, evaluated at time level $j+1$: +#+NAME: eqn:secondderivative \begin{equation} \frac{\partial^2 C^{j+1} }{\partial x^2} = \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} @@ -120,6 +121,7 @@ equations given above leads to the following 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^{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} \nonumber \\ & = \alpha\frac{C^{j+1}_{i-1} - 2C^{j+1}_{i} + C^{j+1}_{i+1}}{\Delta x^2} @@ -176,8 +178,54 @@ Now rearrange terms and substituting with $s_x$ leads to: *TODO* - Tridiagonal matrix filling +** 2D ADI using BTCS scheme +In the previous sections we described the usage of FTCS and BTCS on 1D grids. To +make usage of the BTCS scheme we are following the alternating-direction +implicit method. Therefore we make use of second difference operator defined in +equation [[eqn:secondderivative]] in both $x$ and $y$ direction for half the time +step $\Delta t$. We are denoting the numerator of equation [[eqn:1DBTCS]] as +$\delta^2_x C^n_{ij}$ and $\delta^2_y C^n_{ij}$ respectively with $i$ the +position in $x$, $j$ the position in $y$ direction and $n$ as the current time +step: +\begin{align}\displaystyle +\delta^2_x C^n_{ij} &= C^{n}_{i-1,j} - 2C^{n}_{i,j} + C^{n}_{i+1,j} \nonumber \\ +\delta^2_y C^n_{ij} &= C^{n}_{i,j-1} - 2C^{n}_{i,j} + C^{n}_{i,j+1} +\end{align} + +Assuming a constant $\alpha_x$ and $\alpha_y$ in each direction the equation can +be defined: + +#+NAME: eqn:genADI +\begin{align}\displaystyle +\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{n+1}_{ij}-C^{n+1/2}_{ij}}{\frac{\Delta t}{2}} &= \alpha_y \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n+1}_{ij}\right)}{\Delta y^2} +\end{align} + +Now we will define $s_x$ and $s_y$ respectively as followed: + +\begin{align}\displaystyle +s_x &= \frac{\alpha_x \cdot \frac{\Delta t}{2}}{\Delta x^2} \nonumber \\ +s_y &= \frac{\alpha_y \cdot \frac{\Delta t}{2}}{\Delta y^2} +\end{align} + +Equation [[eqn:genADI]], once developed in the $x$ direction, yields: + +\begin{align}\displaystyle +\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} \right)}{\Delta x^2} + \alpha_x \frac{\left(\delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ +C^{n+1/2}_{ij}-C^n_{ij} &= s_x \delta^2_x C^{n+1/2}_{ij} + s_x \delta^2_y C^{n}_{ij} \nonumber \\ +-C^n_{ij} - s_x \delta^2_y C^{n}_{ij} &= C^{n+1/2}_{ij} + s_x \delta^2_x C^{n+1/2}_{ij} \nonumber \\ +-C^n_{ij} - s_x \left(C^n_{i,j-1}-2C^n_{i,j}+C^n_{i,j+1} \right)&= C^{n+1/2}_{ij} + s_x \left( C^{n+1/2}_{i-1,j} - 2C^{n+1/2}_{i,j} + C^{n+1/2}_{i+1,j} \right) +\end{align} + +All values on the left side of the equation are known and we are solving for +$C^{1/2}$. The calculation in $y$ direction for another $\frac{\Delta t}{2}$ is +similar and can also be easily derived. Therefore we do not show it here. + +*TODO*: +- ADI at boundaries #+LATEX: \clearpage From 57cfa682c576311a3404b1bb3e127811cc7bc722 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 26 Apr 2022 10:07:51 +0200 Subject: [PATCH 11/13] Added 2D ADI scheme for boundaries --- doc/ADI_scheme.org | 39 +++++++++++++++++++++++++++------------ 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 3f3e303..431e0e5 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -186,12 +186,12 @@ implicit method. Therefore we make use of second difference operator defined in equation [[eqn:secondderivative]] in both $x$ and $y$ direction for half the time step $\Delta t$. We are denoting the numerator of equation [[eqn:1DBTCS]] as $\delta^2_x C^n_{ij}$ and $\delta^2_y C^n_{ij}$ respectively with $i$ the -position in $x$, $j$ the position in $y$ direction and $n$ as the current time +position in $x$, $j$ the position in $y$ direction and $T$ as the current time step: \begin{align}\displaystyle -\delta^2_x C^n_{ij} &= C^{n}_{i-1,j} - 2C^{n}_{i,j} + C^{n}_{i+1,j} \nonumber \\ -\delta^2_y C^n_{ij} &= C^{n}_{i,j-1} - 2C^{n}_{i,j} + C^{n}_{i,j+1} +\delta^2_x C^T_{ij} &= C^{T}_{i-1,j} - 2C^{T}_{i,j} + C^{T}_{i+1,j} \nonumber \\ +\delta^2_y C^T_{ij} &= C^{T}_{i,j-1} - 2C^{T}_{i,j} + C^{T}_{i,j+1} \end{align} Assuming a constant $\alpha_x$ and $\alpha_y$ in each direction the equation can @@ -199,8 +199,8 @@ be defined: #+NAME: eqn:genADI \begin{align}\displaystyle -\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ -\frac{C^{n+1}_{ij}-C^{n+1/2}_{ij}}{\frac{\Delta t}{2}} &= \alpha_y \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n+1}_{ij}\right)}{\Delta y^2} +\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{T+1}_{ij}-C^{T+1/2}_{ij}}{\frac{\Delta t}{2}} &= \alpha_y \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T+1}_{ij}\right)}{\Delta y^2} \end{align} Now we will define $s_x$ and $s_y$ respectively as followed: @@ -213,19 +213,34 @@ s_y &= \frac{\alpha_y \cdot \frac{\Delta t}{2}}{\Delta y^2} Equation [[eqn:genADI]], once developed in the $x$ direction, yields: \begin{align}\displaystyle -\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ -\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} \right)}{\Delta x^2} + \alpha_x \frac{\left(\delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ -C^{n+1/2}_{ij}-C^n_{ij} &= s_x \delta^2_x C^{n+1/2}_{ij} + s_x \delta^2_y C^{n}_{ij} \nonumber \\ --C^n_{ij} - s_x \delta^2_y C^{n}_{ij} &= C^{n+1/2}_{ij} + s_x \delta^2_x C^{n+1/2}_{ij} \nonumber \\ --C^n_{ij} - s_x \left(C^n_{i,j-1}-2C^n_{i,j}+C^n_{i,j+1} \right)&= C^{n+1/2}_{ij} + s_x \left( C^{n+1/2}_{i-1,j} - 2C^{n+1/2}_{i,j} + C^{n+1/2}_{i+1,j} \right) +\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} \right)}{\Delta x^2} + \alpha_x \frac{\left(\delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ +C^{T+1/2}_{ij}-C^T_{ij} &= s_x \delta^2_x C^{T+1/2}_{ij} + s_x \delta^2_y C^{T}_{ij} \nonumber \\ +-C^T_{ij} - s_x \delta^2_y C^{T}_{ij} &= C^{T+1/2}_{ij} + s_x \delta^2_x C^{T+1/2}_{ij} \end{align} All values on the left side of the equation are known and we are solving for $C^{1/2}$. The calculation in $y$ direction for another $\frac{\Delta t}{2}$ is similar and can also be easily derived. Therefore we do not show it here. +Substituting $\delta$ in the equation leads to following epxression: -*TODO*: -- ADI at boundaries +\begin{equation} + -C^T_{ij} - s_x \left(C^T_{i,j-1}-2C^T_{i,j}+C^T_{i,j+1} \right) = C^{T+1/2}_{ij} + s_x \left( C^{T+1/2}_{i-1,j} - 2C^{T+1/2}_{i,j} + C^{T+1/2}_{i+1,j} \right) +\end{equation} + +This scheme also only applies to inlet cells without a relation to boundaries. +Fortunately we already derived both cases of outer left and right inlet cell +respectively. Hence we are able to redefine each $\delta^2$ case in x and y +direction, assuming $l_x$ and $l_y$ the be the left boundary value and $r_x$ and +$r_y$ the right one for each direction $x$ and $y$. The equations are exemplary +for timestep $T+1/2$: + +\begin{align}\displaystyle +\delta^2_d C^{T+1/2}_{0,j} &= 2l_x - 3C^{T+1/2}_{0,j} + C^{T+1/2}_{1,j} \nonumber \\ +\delta^2_d C^{T+1/2}_{n,j} &= 2r_x - 3C^{T+1/2}_{n,j} + C^{T+1/2}_{n-1,j} \nonumber \\ +\delta^2_d C^{T+1/2}_{i,0} &= 2l_y - 3C^{T+1/2}_{i,0} + C^{T+1/2}_{i,1} \nonumber \\ +\delta^2_d C^{T+1/2}_{i,n} &= 2r_y - 3C^{T+1/2}_{i,n} + C^{T+1/2}_{i,n-1} +\end{align} #+LATEX: \clearpage From e9a1d067849a2b5dcb49512458c3c4bc005121d6 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 09:55:38 +0200 Subject: [PATCH 12/13] Apply new scheme to model (only 1D) --- src/BTCSDiffusion.cpp | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 1271156..606d6eb 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -239,8 +239,8 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow( A_matrix.insert(1, 1) = 1; } else { double sx = (alpha[0] * time_step) / (dx * dx); - A_matrix.insert(1, 1) = -1. - 2. * sx; - A_matrix.insert(1, 0) = sx; + A_matrix.insert(1, 1) = -1. - 3. * sx; + A_matrix.insert(1, 0) = 2. * sx; A_matrix.insert(1, 2) = sx; } @@ -261,9 +261,9 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow( 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. - 2. * sx; + 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) = sx; + A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx; } A_matrix.insert(A_size - 1, A_size - 1) = 1; @@ -294,15 +294,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]); - } + // this is not correct currently.We will fix this when we are able to define + // FLUX boundary conditions + b_vector[0] = + (left_constant ? left.value : getBCFromFlux(left, c[0], alpha[0])); - if (!right_constant) { - b_vector[b_size - 1] = getBCFromFlux(right, c[size - 1], alpha[size - 1]); - } + b_vector[b_size - 1] = + (right_constant ? right.value + : getBCFromFlux(right, c[size - 1], alpha[size - 1])); } void Diffusion::BTCSDiffusion::setTimestep(double time_step) { From d35a27f54a4adfa6b8b1ca5e59302f24b51ff9d1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 26 Apr 2022 13:29:12 +0200 Subject: [PATCH 13/13] Apply 2D scheme to model --- src/BTCSDiffusion.cpp | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 606d6eb..4f7b109 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -141,11 +141,10 @@ 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++) { @@ -184,12 +183,17 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c, // 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 @@ -210,12 +214,17 @@ 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;