The general convection-diffusion equation is

\frac{\partial c}{\partial t} = \nabla\cdot (D \nabla c) - \nabla\cdot (\vec v\:c) + R

where

• c\: is concentration.
• D is diffusivity.
• \vec v\: is the velocity field.
• R\: describes "sources" or "sinks" of concentration, e.g. chemical reactions.

I want to model this in two dimensions with constant D and \vec v, and R = 0

\frac{\partial c}{\partial t} = \left(D\frac{\partial^2 c}{\partial x^2} + D\frac{\partial^2 c}{\partial y^2}\right) - \left(v_x\frac{\partial c}{\partial x} + v_y\frac{\partial c}{\partial y}\right)

And in one dimension in the general case of D, v, and R being dependent on time and location

\frac{\partial c}{\partial t} = \frac{\partial (D\frac{\partial c}{\partial x})}{\partial x} - \frac{\partial (vc)}{\partial x} + R.

## Two Dimensions

Note that for this section, for the functions of time and location defined below, the following shorthands are equivalent

f = f(x, y) = f(t, x, y)
f^*(t, x, y) = f(t + \Delta t, x, y)

Using the Crank-Nicholson finite difference method, the 2-D equation becomes

\begin{array}{lrrl} \frac{c^* - c}{\Delta t} &=\frac12\Bigg( \Bigg[& \bigg(& D\frac{ c^*(x + \Delta x, y) - 2c^* + c^*(x - \Delta x, y)}{(\Delta x)^2} \\[1em] &&+& D\frac{ c^*(x, y + \Delta y) - 2c^* + c^*(x, y - \Delta y)}{(\Delta y)^2} \bigg) \\[1em] &&- \bigg(& v_x\frac{ c^*(x + \Delta x, y) - c^*(x - \Delta x, y)}{2\Delta x} \\[1em] &&+& v_y\frac{ c^*(x, y + \Delta y) - c^*(x, y - \Delta y)}{2\Delta y} \bigg) \Bigg] \\[1em] &+ \Bigg[& \bigg(& D\frac{ c(x + \Delta x, y) - 2c + c(x - \Delta x, y)}{(\Delta x)^2} \\[1em] &&+& D\frac{ c(x, y + \Delta y) - 2c + c(x, y - \Delta y)}{(\Delta y)^2} \bigg) \\[1em] &&- \bigg(& v_x\frac{ c(x + \Delta x, y) - c(x - \Delta x, y)}{2\Delta x} \\[1em] &&+& v_y\frac{ c(x, y + \Delta y) - c(x, y - \Delta y)}{2\Delta y} \bigg) \Bigg] \Bigg). \end{array}

This can be rearranged to

\begin{array}{rl} (1 + 2A_x + 2A_y)c^*& - (A_x+B_x)c^*(x+\Delta x,y) \\[1em]& - (A_x-B_x)c^*(x-\Delta x,y) \\[1em]& - (A_y+B_y)c^*(x,y+\Delta y) \\[1em]& - (A_y-B_y)c^*(x,y-\Delta y) \\[1em] = & s \end{array}

where

\begin{array}{rl} s = (1 - 2A_x - 2A_y)c& + (A_x+B_x)c(x+\Delta x,y) \\[1em]& + (A_x-B_x)c(x-\Delta x,y) \\[1em]& + (A_y+B_y)c(x,y+\Delta y) \\[1em]& + (A_y-B_y)c(x,y-\Delta y) \end{array}
\begin{array}{lclclcl} A_x &=& \frac{v_x\Delta t}{4\Delta x} && B_x &=& \frac{D\Delta t}{2(\Delta x)^2} \\[1em] A_y &=& \frac{v_y\Delta t}{4\Delta y} && B_y &=& \frac{D\Delta t}{2(\Delta y)^2} \\[1em] \end{array}

This is actually a two-dimensional system of equations at each timestep, which, for a system simulating W\times H "active" (non boundary condition) nodes, can be represented by the tensor equation

\mathbf T\:\mathbf C^* = \mathbf G

where

• \mathbf T is a constant (time-independent) four-dimensional tensor with shape {\langle W+2},{H+2},{W+2},{H+2\rangle}. We will get back to this.

• \mathbf C^* is a (W+2)\times(H+2) matrix representing the state at the next timestep (the state being calculated).

\mathbf C^* = \begin{bmatrix} c^*(0, 0) & \cdots & c^*((W+1)\Delta x,0) \\ \vdots & \ddots & \vdots \\ c^*(0,(H+1)\Delta y) & \cdots & c^*((W+1)\Delta x,(H+1)\Delta y) \end{bmatrix}
• \mathbf G is a (W+2)\times(H+2) matrix representing the (processed) state at the current timestep (the known state).

\mathbf G = \begin{bmatrix} g(0, 0) & \cdots & g((W+1)\Delta x,0) \\ \vdots & \ddots & \vdots \\ g(0,(H+1)\Delta y) & \cdots & g((W+1)\Delta x,(H+1)\Delta y) \end{bmatrix}
g = \begin{cases} s &\text{if}\:(1\le x\le W\Delta x) \land (1\le y\le H\Delta y) \\ b^* &\text{if}\:(1\le x\le W\Delta x) \oplus (1\le y\le H\Delta y) \\ 0 &\text{otherwise} \end{cases}

g and s are functions of time and location defined above, and b is some function of time and location representing the boundary conditions.

c^* = g\;\text{if}\;\lnot((1\le x\le W\Delta x) \land (1\le y\le H\Delta y))

\mathbf T really did my head in, so I thought it would be useful to write everything out for a simple case first. For W = 2 and H = 2, with dots representing zeros, \mathbf T = {}

\begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix} \hspace{1.8em} \begin{bmatrix} \cdot & 1 & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix} \hspace{2.7em} \begin{bmatrix} \cdot & \cdot & 1 & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix} \hspace{1.8em} \begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix}
\begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix} \begin{bmatrix} \cdot & {\small k_{-y}} & \cdot & \cdot \\ {\small k_{-x}} & {\small k_\varnothing} & {\small k_{+x}} & \cdot \\ \cdot & {\small k_{+y}} & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix} \begin{bmatrix} \cdot & \cdot & {\small k_{-y}} & \cdot \\ \cdot & {\small k_{-x}} & {\small k_\varnothing} & {\small k_{+x}} \\ \cdot & \cdot & {\small k_{+y}} & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix} \begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & 1 \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix}
\begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix} \begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & {\small k_{-y}} & \cdot & \cdot \\ {\small k_{-x}} & {\small k_\varnothing} & {\small k_{+x}} & \cdot \\ \cdot & {\small k_{+y}} & \cdot & \cdot \end{bmatrix} \begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & {\small k_{-y}} & \cdot \\ \cdot & {\small k_{-x}} & {\small k_\varnothing} & {\small k_{+x}} \\ \cdot & \cdot & {\small k_{+y}} & \cdot \end{bmatrix} \begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & 1 \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix}
\begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix} \hspace{1.8em} \begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & 1 & \cdot & \cdot \end{bmatrix} \hspace{2.7em} \begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & 1 & \cdot \end{bmatrix} \hspace{1.8em} \begin{bmatrix} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{bmatrix}

where

k_\varnothing = (1 + 2A_x + 2A_y)
\begin{array}{lclclcl} k_{+x} &=& -(A_x + B_x) && k_{-x} &=& -(A_x - B_x) \\[1em] k_{+y} &=& -(A_y + B_y) && k_{-y} &=& -(A_y - B_y) \\[1em] \end{array}

Okay, now this is a bit easier to visualize and understand. The general formula for an element of \mathbf T is

\mathbf T_{i,j,m,n} = \begin{cases} k_\varnothing &\text{if}\:a \land i = m \land j = n \\ k_{+x} &\text{if}\:a \land i + 1 = m \land j = n \\ k_{-x} &\text{if}\:a \land i - 1 = m \land j = n \\ k_{+y} &\text{if}\:a \land i = m \land j + 1 = n \\ k_{-y} &\text{if}\:a \land i = m \land j - 1 = n \\ 1 &\text{if}\:((1\le i\le W) \oplus (1\le j\le H)) \land i = m \land j = n \\ 0 &\text{otherwise} \end{cases}

where a = ((1\le i\le W) \land (1\le j\le H)).

And there we go! We have a full description for the finite difference method for the two-dimensional convection-diffusion equation with constant D and \vec v, and R = 0. I can now write a program to generate \mathbf T and \mathbf G, and solve for \mathbf C^* with something like numpy's linalg.tensorsolve. But I'll leave that for another day.