tug/doc/ADI_scheme.org

Numerical solution of diffusion equation in 2D with ADI Scheme

Homogeneous diffusion in 1D

Finite differences with nodes as cells' centres

The 1D diffusion equation for spatially constant diffusion coefficients $\alpha$ 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/9815ebce9c7ed97d3c192a1f638aa54d5b5a5428/doc/images/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}

Heterogeneous diffusion

If the diffusion coefficient $\alpha$ is spatially variable, eqn:1 can be rewritten:

\begin{align} \frac{\partial C }{\partial t} & = \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x} \right) \end{align}

Discretization of the equation using chain rule

\noindent From the product rule for derivatives we obtain:

\begin{equation} \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right) = \frac{\partial \alpha}{\partial x} \cdot \frac{\partial C}{\partial x} + \alpha \frac{\partial^2 C }{\partial x^2} \end{equation}

\noindent Using a spatially centred second order finite difference approximation at $x=x_i$ for both $\alpha$ and $C$, we have

\begin{align} \frac{\partial \alpha}{\partial x} \cdot \frac{\partial C}{\partial x} + \alpha \frac{\partial^2 C }{\partial x^2} & \simeq \frac{\alpha_{i+1} - \alpha_{i-1}}{2\Delta x}\cdot\frac{C_{i+1} - C_{i-1}}{2\Delta x} + \alpha_i \frac{C_{i+1} - 2 C_{i} + C_{i-1}}{\Delta x^2} \\ \nonumber & = \frac{1}{\Delta x^2} \frac{\alpha_{i+1} - \alpha_{i-1}}{4}(C_{i+1} - C_{i-1}) + \frac{\alpha_{i}}{\Delta x^2}(C_{i+1}-2C_i+C_{i-1})\\ \nonumber & = \frac{1}{\Delta x^2} \left\{A C_{i+1} -2\alpha_i C_i + AC_{i-1})\right\} \end{align}

\noindent having set \[ A = \frac{\alpha_{i+1}-\alpha_{i-1}}{4} + \alpha_i \]

\noindent In 2D the ADI scheme eqn:genADI with heterogeneous diffusion coefficients can thus be written:

\begin{equation} \systeme{ \displaystyle \frac{C^{t+1/2}_{i,j}-C^{t }_{i,j}}{\Delta t/2} = \displaystyle \frac{\partial}{\partial x} \left( \alpha^x_{i,j} \frac{\partial C^{t+1/2}_{i,j}}{\partial x}\right) + \frac{\partial}{\partial y} \left( \alpha^y_{i,j} \frac{\partial C^{t }_{i,j}}{\partial y}\right), \displaystyle \frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = \displaystyle \frac{\partial}{\partial x} \left( \alpha^x_{i,j} \frac{\partial C^{t+1/2}_{i,j}}{\partial x}\right) + \frac{\partial}{\partial y} \left( \alpha^y_{i,j} \frac{\partial C^{t+1}_{i,j}}{\partial y}\right) } \end{equation}

\noindent We define for compactness $S_x=\frac{\Delta t}{2\Delta x^2}$ and $S_y=\frac{\Delta t}{2\Delta y^2}$ and

\begin{equation} \systeme{ \displaystyle A_{i,j} = \displaystyle \frac{\alpha^x_{i+1,j} -\alpha^x_{i-1,j }}{4} + \alpha^x_{i,j}, \displaystyle B_{i,j} = \displaystyle \frac{\alpha^y_{i, j+1}-\alpha^y_{i ,j-1}}{4} + \alpha^y_{i,j} } \end{equation}

\noindent Plugging eq. (eqn:hetdiff_fd) into the first of equations (eqn:genADI_het) - so called "sweep by x" - and putting all implicit terms (at time level $t+1/2$) on the left hand side we obtain:

\begin{equation}\displaystyle

\begin{split} -S_x A_{i,j} C^{t+1/2}_{i+1,j} + (1 + 2S_x\alpha^x_{i,j})C^{t+1/2}_{i,j} - S_x A_{i,j}C^{t+1/2}_{i-1,j} = \\ S_y B_{i,j} C^{t }_{i,j+1} + (1 - 2S_y\alpha^y_{i,j})C^{t }_{i,j} + S_y B_{i,j}C^{t }_{i,j-1} \end{split}

\end{equation}

\noindent In the same way for the second of eq. eqn:genADI_het we have:

\begin{equation}\displaystyle

\begin{split} -S_y B_{i,j} C^{t+1 }_{i,j+1} + (1 + 2S_y\alpha^y_{i,j})C^{t+1 }_{i,j} - S_y B_{i,j}C^{t+1 }_{i,j-1} = \\ S_x A_{i,j} C^{t+1/2}_{i+1,j} + (1 - 2S_x\alpha^x_{i,j})C^{t+1/2}_{i,j} + S_x A_{i,j}C^{t+1/2}_{i-1,j} \end{split}

\end{equation}

\noindent If the diffusion coefficients are constant, $A_{i,j}=B_{i,j}=\alpha$ and the scheme reverts to the homogeneous case. Problem with this discretization is that the terms in $A_{ij}$ and $B_{ij}$ can be negative depending on the derivative of the diffusion coefficient, resulting in unphysical values for the concentrations.

Direct discretization

As noted in literature (LeVeque and Numerical Recipes) a better way is to discretize directly the physical problem (eqn:hetdiff) at points halfway between grid points:

\begin{align*} \begin{cases} \displaystyle \alpha(x_{i+1/2}) \frac{\partial C }{\partial x}(x_{i+1/2}) & \displaystyle = \alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) \\ \displaystyle \alpha(x_{i-1/2}) \frac{\partial C }{\partial x}(x_{i-1/2}) & \displaystyle = \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) \end{cases} \end{align*}

\noindent A further differentiation gives us the centered approximation of $\frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right)$:

\begin{align*} \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right)(x_i) & \simeq \frac{1}{\Delta x}\left[\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) - \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) \right]\\ &\displaystyle =\frac{1}{\Delta x^2} \left[ \alpha_{i+1/2}C_{i+1} - (\alpha_{i+1/2}+\alpha_{i-1/2}) C_{i} + \alpha_{i-1/2}C_{i-1}\right] \end{align*}

\noindent The ADI scheme with this approach becomes:

\begin{equation} \left\{ \begin{aligned} \frac{C^{t+1/2}_{i,j}-C^{t }_{i,j}}{\Delta t/2} = & \frac{1}{\Delta x^2} \left[ \alpha_{i+1/2,j} C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C^{t+1/2}_{i,j} + \alpha_{i-1/2,j} C^{t+1/2}_{i-1,j}\right] + \\ & \frac{1}{\Delta y^2} \left[ \alpha_{i,j+1/2}C^{t}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^t_{i} + \alpha_{i,j-1/2}C^{t}_{i,j-1}\right]\\ \frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = & \frac{1}{\Delta y^2} \left[ \alpha_{i+1/2,j}C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C_{i} + \alpha_{i-1/2,j}C^{t+1/2}_{i-1,j}\right] + \\ & \frac{1}{\Delta x^2} \left[ \alpha_{i,j+1/2}C^{t+1}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^{t+1}_{i,j} + \alpha_{i,j-1/2}C^{t+1}_{i,j-1}\right] \end{aligned} \right. \end{equation}

\noindent Doing the usual algebra and separating implicit from explicit terms, the two sweeps become:

\begin{equation} \left\{ \begin{aligned} -S_x \alpha^x_{i+1/2,j} C^{t+1/2}_{i+1,j} + (1 + S_x(\alpha^x_{i+1/2,j}+ \alpha^x_{i-1/2,j}))C^{t+1/2}_{i,j} - S_x \alpha^x_{i-1/2,j} C^{t+1/2}_{i-1,j} = \quad & \\ S_y \alpha^y_{i,j+1/2} C^{t }_{i,j+1} + (1 - S_y(\alpha^y_{i,j+1/2}+ \alpha^y_{i,j-1/2}))C^{t }_{i,j} + S_y \alpha^y_{i,j-1/2} C^{t }_{i,j-1} & \\[1em] -S_y \alpha^y_{i,j+1/2} C^{t+1 }_{i,j+1} + (1 + S_y(\alpha^y_{i,j+1/2}+ \alpha^y_{i,j-1/2}))C^{t+1 }_{i,j} - S_y \alpha^y_{i,j-1/2} C^{t+1 }_{i,j-1} = \qquad & \\ S_x \alpha^x_{i+1/2,j} C^{t+1/2}_{i+1,j} + (1 - S_x(\alpha^x_{i+1/2,j}+ \alpha^x_{i-1/2,j}))C^{t+1/2}_{i,j} + S_x \alpha^x_{i-1/2,j} C^{t+1/2}_{i-1,j} \end{aligned} \right. \end{equation}