tug/doc/Implementation.org

Implementation of 2D ADI BTCS

Introduction

As described in /max/tug/src/commit/9d663d8140b17e0df7095cf6fd13a891092dce72/doc/ADI_scheme.org the 1D-BTCS diffusion and the 2D-BTCS diffusion using the ADI scheme are currently implemented. In the following document the implementation is described.

To access columns and rows of input vectors and matrices with ease and to use build-in solver algorithms we're using the Eigen library inside the library. However, this not imply to use Eigen as the on-top data structure of the application. The API expects pointer to a continuous memory area, where concentrations, alpha and boundary conditions are stored. Those memory areas are mapped to the according Eigen data structure then.

Solving a linear equation system requires an $A$ matrix and $b$ and $x$ vectors. All of them are also implemented using Eigen data structures, whereby the $A$ matrix is implemented as a sparse matrix.

In the following all indices will start at 0.

Internal data structures

During instantiating of an object the dimension size is required. According to that parameter all other parameters are stored inside the instance of the diffusion module and are set with the according getter and setters (see code documentation for more details). Stored values are:

• grid dimension
• time step $\Delta t$
• number of grid cells $N$ in each direction
• domain size $x$ or $y$ resp. in each direction
• spatial discretization in each direction $\Delta x$ or $\Delta y$ resp.

Whenever a value is changed through according setter, the private member function updateInterals() is called, updating all other dependent values.

Diffusion in 1D

As described above to call simulate() properly, pointers to continuous memory containing the field with its concentrations, alphas and boundary conditions are required. The boundary condition data type is defined in the BoundayCondition.hpp header file. There are three assignable types:

• BC_CONSTANT - constant/Dirichlet boundary condition with value set to a fixed value.
• BC_CLOSED - closed/Neumann boundary condition with value set to a gradient of 0
• BC_FLUX - flux/Chauchy boundary with value set to user defined gradient

\noindent See the code documentation for more details on that. Defining the number of grid cells as $N$, the input values should have the following dimensions:

• vector: field c $\to N$
• vector: alpha alpha $\to N$
• vector: boundary condition bc $\to N+2$, where bc[0] and bc[N+1] describing the actual boundary condition on the left/right hand side. All other values can be used to define constant inlet cells.

All input values are mapped to Eigen vectors and passed to the simulate1D() function. With the only row (since we are just having a 1D vector) we will start the simulation by calling simulate_base(). There, a SparseMatrix A with a dimension of $N+1 \times N+1$ and a VectorXd b with a dimension of $N+1$ are created with, filled with according values and used to solve the LES by the built-in SparseLU solver of Eigen. The results are written to the VectorXd x and copied back to the memory area of c.

Filling of sparse matrix A

To fill the A matrix various information about the field are needed. This includes:

• alpha for each cell of the field
• bc for each cell of the field and boundary conditions
• size of the field
• dx as the $\Delta x$
• time_step as the current time step size $\Delta t$

We are using indices $j$ from $0 \to N+1$ to address the correct field of the boundary condition and $k$ from $0 \to N-1$ to address the according field of alpha

First of all the left and right boundary condition are extracted from bc evaluated. Then we will set the left boundary condition to the linear equation system, which is the left boundary condition, to 1 ($A(0,0) = 1$). Now solving for the leftmost inlet is required. Each $sx = (\alpha_k \cdot \Delta t)/ \Delta x^2$ is stored in the same named variable for each index pointing to the current cell of the input/output field.

• $A(1,0) = 2 \cdot sx$
• $A(1,1) = -1 -3 \cdot sx$
• $A(1,2) = sx$

After that we will iterate until $k = N-2$ or $j = N-1$ is reached respectively and setting A with alpha[k] and bc[j] to the following values:

• $A(j,j-1) = sx$
• $A(j,j) = -1 -2 \cdot sx$
• $A(j,j+1) = sx$

Now the equation for the rightmost inlet cell is applied:

• $A(N,N-1) = sx$
• $A(N,N) = -1 -3 \cdot sx$
• $A(N,N+1) = 2 \cdot sx$

The right boundary condition is also set to $A(N+1,N+1) = 1$.

A special case occurs if one cell of the input/output field is set to a constant value by setting the according bc field to BC_CONSTANT. Then, just the according element of the matrix is set to 1. For example cell $k=4$ is marked as constant will lead to $j = 5$ and $A(j,j) = 1$.

Filling of vector b

The right hand side of the equation system requires the following input values:

• c as the input field
• alpha for each cell of the grid
• bc for each cell of the field and boundary conditions
• d_ortho as the explicit parameter of left hand side of the equation for each cell of the field (only interesting if using 2D ADI, otherwise this parameter is set to 0 for all cells)
• size of the field
• dx as the $\Delta x$
• time_step as the current time step size $\Delta t$

Now we are able to iterating through the indices of $j$ from $0 \to N-1$ and filling up the b vector by applying the following rule:

\begin{equation}\displaystyle b_j+1 = \begin{cases} bc_j+1.\text{value} & bc_j+1.\text{type} = \texttt{BC\_CONSTANT}
-\texttt{c}_j - \texttt{d\_ortho}_j & \text{else}

\end{cases}

\end{equation}

\noindent With

\begin{equation}\displaystyle \texttt{d\_ortho}_j = Δ t ⋅ α_j ⋅ (\frac{\bot d_j}{Δ x^2})

\end{equation}

Both boundary conditions on the left and right side are assigned by either using a constant value or a derived value from a given gradient, depending on the assigned boundary condition type. Until now, only gradients of 0 (closed boundary condition) are supported. This leads to:

\begin{equation}\displaystyle b₀ = \begin{cases} bc_0.\text{value} & bc_0.\text{type} = \texttt{BC\_CONSTANT}
c_0 & \text{else}

\end{cases}

\end{equation}

\begin{equation}\displaystyle b_N+1 = \begin{cases} bc_N+1.\text{value} & bc_N+1.\text{type} = \texttt{BC\_CONSTANT}
c_N-1 & \text{else}

\end{cases}

\end{equation}

Diffusion in 2D

Simulating a 2D diffusion works in overall just like the 1D diffusion. Input parameters are should now be matrices or more precisly arrays, which can be interpreted as matrices, no matter if stored column or row wise. Defining the number of grid cells now as $N times M$, the input parameters should have the following dimensions:

• vector: field c $\to N \times M$
• vector: alpha alpha $\to N \times M$
• vector: boundary condition bc $\to N+2 \times M+2$

Those input pointers will now be interpreted as Eigen matrices for further use. In simulate2D() the matrix d_ortho is calculated in y-direction. Then, for each row of matrix c simulate_base() is called. The according row of alpha, bc and d_ortho are passed too and the time step is divided by two. Now the 2D problem is divided into an 1D problem. After each row the new resulting field is written back into the original memory area and serves as the input of the second simulation in the y-direction.

For this, d_ortho will be recalculated, now in x-direction. Now, for every column of c, the correspondending column alpha, bc and d_ortho is passed to simulate_base() with also half the time step.

Finally, the result of the simulation of each column get written back to c and the 2D simulation is finished.

Calculation of d_ortho

The function cal_d_ortho returns a matrix with the size of $N \times M$ with each cell holding the $\bot d$ value of a cell from c at the same index. This is done by first iterating over the first row from $0 \to N-1$. There, we apply the following rule:

\begin{equation}\displaystyle \bot d_0,j = \begin{cases} Δ t ⋅ α_0,j ⋅ ≤ft(\frac{2bc_0,j+1.\text{value} - 3 ⋅ c_0,j + c_1,j}{Δ x^2}\right) & bc_0,j+1 = \texttt{BC\_CONSTANT}
Δ t ⋅ α_0,j ⋅ ≤ft(\frac{c_0,j - 2 ⋅ c_0,j + c_1,j}{Δ x^2}\right) & \text{else}

\end{cases}

\end{equation}

\noindent Then we succeeding to the inner cells with $i = 1,\dots,M-2$ and $j = 0,\dots,N-1$:

\begin{equation}\displaystyle \bot d_i,j = Δ t ⋅ α_i,j ⋅ ≤ft(\frac{c_i-1,j - 2 ⋅ c_i,j + c_i+1,j}{Δ x^2}\right)

\end{equation}

\noindent And lastly calculating the last row by applying the following rule for each $j = 0,\dots,N-1$:

\begin{equation}\displaystyle \bot d_M-1,j = \begin{cases} Δ t ⋅ α_M-1,j ⋅ ≤ft(c_M-2,j -3 ⋅ c_M-1,j + 2\frac{bc_M+1,j+1.\text{value}}{Δ x^2}\right) & bc_M+1,j+1 = \texttt{BC\_CONSTANT}
Δ t ⋅ α_M-1,j ⋅ ≤ft(\frac{c_M-2,j - 2 ⋅ c_M-1,j + c_M-1,j}{Δ x^2}\right) & \text{else}

\end{cases}

\end{equation}