tug/doc/ADI_scheme.org

#+TITLE: Adi 2D 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

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 discretize it following a Forward Time, Centered Space finite
difference scheme where the nodes correspond to the centers of a 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$.

We discretize [[eqn:1]] as following, for each index i in 1, \dots, n-1
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}
\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:

#+NAME: eqn:3
\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}
\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}_{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)
\end{align}

Which leads to the following term applicable to our approach:

#+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
\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:5]]
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.