As far as literature searches go in the direction of cataloguing finite difference methods explicitly, there appears to be no descriptive representation of reaction–diffusion equations in their full three spatial dimensions. It is therefore our task (albeit tedious) to assemble the matrices and vectors for the model.

We extend the two-dimensional operator formulation to the three-dimensional case by introducing an additional spatial index n \in \{0,1,\dots,L\} with step size \Delta z. The flux operators (\epsilon^{+}_x + \epsilon^{+}_y) are augmented by the z-direction term \epsilon^{+}_z, so that the discrete diffusion operator acting on the substrate S is given by the 3D seven-point stencil.

3D Stencil Formulation

(1)   \[ \begin{aligned} (\epsilon^{+}_x + \epsilon^{+}_y + \epsilon^{+}_z) S^k_{l,m,n} &= a^{-}_{l,m,n;x}\, S^{k+1}_{l-1,m,n} + a^{+}_{l,m,n;x}\, S^{k+1}_{l+1,m,n} \\ &\quad + a^{-}_{l,m,n;y}\, S^{k+1}_{l,m-1,n} + a^{+}_{l,m,n;y}\, S^{k+1}_{l,m+1,n} \\ &\quad + a^{-}_{l,m,n;z}\, S^{k+1}_{l,m,n-1} + a^{+}_{l,m,n;z}\, S^{k+1}_{l,m,n+1} \\ &\quad + b_{l,m,n}\, S^{k+1}_{l,m,n}. \end{aligned} \]

Coefficient Definitions

    \[ \begin{aligned} a^{-}_{l,m,n;x} &= \frac{\Delta t}{(\Delta x)^2}\, D\!\left(\frac{M^k_{l-1,m,n} + M^k_{l,m,n}}{2}\right),\\ a^{+}_{l,m,n;x} &= \frac{\Delta t}{(\Delta x)^2}\, D\!\left(\frac{M^k_{l,m,n} + M^k_{l+1,m,n}}{2}\right),\\[6pt] a^{-}_{l,m,n;y} &= \frac{\Delta t}{(\Delta y)^2}\, D\!\left(\frac{M^k_{l,m-1,n} + M^k_{l,m,n}}{2}\right),\\ a^{+}_{l,m,n;y} &= \frac{\Delta t}{(\Delta y)^2}\, D\!\left(\frac{M^k_{l,m,n} + M^k_{l,m+1,n}}{2}\right),\\[6pt] a^{-}_{l,m,n;z} &= \frac{\Delta t}{(\Delta z)^2}\, D\!\left(\frac{M^k_{l,m,n-1} + M^k_{l,m,n}}{2}\right),\\ a^{+}_{l,m,n;z} &= \frac{\Delta t}{(\Delta z)^2}\, D\!\left(\frac{M^k_{l,m,n} + M^k_{l,m,n+1}}{2}\right). \end{aligned} \]

(2)   \[ b_{l,m,n} = 1 - \left( a^{-}_{l,m,n;x} + a^{+}_{l,m,n;x} + a^{-}_{l,m,n;y} + a^{+}_{l,m,n;y} + a^{-}_{l,m,n;z} + a^{+}_{l,m,n;z} \right). \]

3D Block-Matrix Assembly

(3)   \[ \mathcal{A}^k_Y = \begin{bmatrix} A^k_Y & H^{+}_1 & 0 & \cdots & 0 \\ H^{-}_2 & A^k_Y & H^{+}_2 & \cdots & 0 \\ 0 & H^{-}_3 & A^k_Y & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & H^{+}_{L-1} \\ 0 & \cdots & 0 & H^{-}_L & A^k_Y \end{bmatrix}. \]

Global Coupled System

(4)   \[ \mathcal{M}^k = \begin{bmatrix} \mathcal{A}^k_S & \mathcal{Z}^k_1 & 0 & \mathcal{H}^k \\ \mathcal{G}^k_1 & \mathcal{A}^k_X & 0 & 0 \\ 0 & \mathcal{Z}^k_2 & \mathcal{A}^k_I & 0 \\ \mathcal{G}^k_2 & 0 & 0 & \mathcal{A}^k_E \end{bmatrix}. \]

Kronecker Representation

(5)   \[ \mathcal{A}^k_Y = I_z \otimes A^k_Y + T_z \otimes I_{xy}, \]

where I_z is the identity in the z-direction, I_{xy} is the identity in the (x,y)-plane, and T_z contains the coefficients a^{\pm}_{l,m,n;z}. This form preserves sparsity and allows direct use of sparse BLAS.

Leave a Reply

Your email address will not be published. Required fields are marked *