tug/doc/ADI_scheme.org

Numerical solution of diffusion equation in 2D with ADI Scheme

Diffusion in 1D

Finite differences with nodes as cells' centres

The 1D diffusion equation is:

\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:

/max/tug/src/commit/8f3b19b63afdb6cd0b78b35bb4642969296466e6/doc/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:

\begin{equation}\displaystyle \frac{C_i^t+1 -C_i^t}{Δ t} = α\frac{\frac{C^t_i+1-C^t_i}{Δ x}-\frac{C^t_i-C^t_i-1}{Δ x}}{Δ 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:

\begin{equation}\displaystyle \frac{C_0^t+1 -C_0^t}{Δ t} = α\frac{\frac{C^t₁-C^t₀}{Δ x}- \frac{C^t₀-l}{\frac{Δ x}{2}}}{Δ x}

\end{equation}

This expression, once developed, yields:

\begin{align}\displaystyle C_0^t+1 & = C_0^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( C^t₁-C^t₀- 2 C^t₀+2l \right) \nonumber
& = C_0^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( C^t₁- 3 C^t₀ +2l \right)

\end{align}

In case of constant right boundary, the finite difference of point $C_n$ - calling $r$ the right boundary value - is:

\begin{equation}\displaystyle \frac{C_n^t+1 -C_n^t}{Δ t} = α\frac{\frac{r - C^t_n}{\frac{Δ x}{2}}- \frac{C^t_n-C^t_n-1}{Δ x}}{Δ x}

\end{equation}

Which, developed, gives

\begin{align}\displaystyle C_n^t+1 & = C_n^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( 2 r - 2 C^t_n -C^t_n + C^t_n-1 \right) \nonumber
& = C_n^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( 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:

\begin{equation}\displaystyle C_n^t+1 = C_n^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ (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$:

\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{align}\displaystyle \frac{C_i^t+1 - C_i^t}{Δ t} & = α\frac{\frac{C^t+1_i+1-C^t+1_i}{Δ x}-\frac{C^t+1_i-C^t+1_i-1}{Δ x}}{Δ x} \nonumber
& = α\frac{C^t+1_i-1 - 2C^t+1_i + C^t+1_i+1}{Δ 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}{Δ t} = α\frac{\frac{C^t+1₁-C^t+1₀}{Δ x}- \frac{C^t+1₀-l}{\frac{Δ x}{2}}}{Δ x}

\end{equation}

This expression, once developed, yields:

\begin{align}\displaystyle C_0^t+1 & = C_0^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( C^t+1₁-C^t+1₀- 2 C^t+1₀+2l \right) \nonumber
& = C_0^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( C^t+1₁- 3 C^t+1₀ +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) ⋅ l + (-1 - 3s_x) ⋅ C^t+1_0 + s_x ⋅ 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}{Δ t} = α\frac{\frac{r-C^t+1_n}{\frac{Δ x}{2}}- \frac{C^t+1_n-C^t+1_n-1}{Δ x}}{Δ x}

\end{equation}

This expression, once developed, yields:

\begin{align}\displaystyle C_n^t+1 & = C_n^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( 2r - 2C^t+1_n - C^t+1_n + C^t+1_n-1 \right) \nonumber
& = C_0^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( 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 ⋅ C^t+1_n-1 + (-1 - 3s_x) ⋅ C^t+1_n + (2s_x) ⋅ r

\end{equation}

TODO

• Tridiagonal matrix filling

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:

\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{∂^2 C^t_i,j}{∂ x^2} &= \frac{C^t_i-1,j - 2C^t_i,j + C^t_i+1,j}{Δ x^2}
\frac{∂^2 C^t_i,j}{∂ y^2} &= \frac{C^t_i,j-1 - 2C^t_i,j + C^t_i,j+1}{Δ y^2}

\end{align}

The ADI scheme is formally defined by the equations:

\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:

\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:

\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]$

\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]$:

\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$:

\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

\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}

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

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}$