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}
\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}
+
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}
+
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}
+
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}
+
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}
\begin{array}{rl}
(1 + 2A_x + 2A_y)c^*&
- (A_x+B_x)c^*(x+\Delta x,y)
- (A_x-B_x)c^*(x-\Delta x,y) \\[1em]&
- (A_y+B_y)c^*(x,y+\Delta y)
- (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}{rl}
s = (1 - 2A_x - 2A_y)c&
+ (A_x+B_x)c(x+\Delta x,y)
+ (A_x-B_x)c(x-\Delta x,y) \\[1em]&
+ (A_y+B_y)c(x,y+\Delta y)
+ (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.