I don't understand the different behaviour of the advection-diffusion equation when I apply different boundary conditions. My motivation is the simulation of a real physical quantity (particle density) under diffusion and advection. Particle density should be conserved in the interior unless it flows out from the edges. By this logic, if I enforce Neumann boundary conditions the ends of the system such as $\frac{\partial \phi}{\partial x}=0$ (on the left and the right sides) then the system should be "closed" i.e. if the flux at the boundary is zero then no particles can escape.

For all the simulations below, I have applied the Crank-Nicolson discretization to the advection-diffusion equation and all simulation have $\frac{\partial \phi}{\partial x}=0$ boundary conditions. However, for the first and last rows of the matrix (the boundary condition rows) I allow $\beta$ to be changed independently of the interior value. This allows the end points to be fully implicit.

Below I discuss 4 different configurations, only one of them is what I expected. At the end I discuss my implementation.

Diffusion only limit

Here the advection terms are turned off by setting the velocity to zero.

Diffusion only, with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at all points

The quantity is not conserved as can be seen by the pulse area reducing.

Diffusion only, with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at interior points, and $\boldsymbol{\beta}$=1 (full implicit) at the boundaries

By using fully implicit equation on the boundaries I achieve what I expect: no particles escape. You can see this by the area being conserved as the particle diffuse. Why should the choice of $\beta$ at the boundary points influence the physics of the situation? Is this a bug or expected?

Diffusion and advection

When the advection term is included, the value of $\beta$ at the boundaries does not seem to influence the solution. However, for all cases when the boundaries seem to be "open" i.e. particles can escape the boundaries. Why is this the case?

Advection and Diffusion with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at all points

Advection and Diffusion with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at interior points, and $\boldsymbol{\beta}$=1 (full implicit) at the boundaries

Implementation of the advection-diffusion equation

Starting with the advection-diffusion equation,

$ \frac{\partial \phi}{\partial t} = D\frac{\partial^2 \phi}{\partial x^2} + \boldsymbol{v}\frac{\partial \phi}{\partial x} $

Writing using Crank-Nicolson gives,

$ \frac{\phi_{j}^{n+1} - \phi_{j}^{n}}{\Delta t} = D \left[ \frac{1 - \beta}{(\Delta x)^2} \left( \phi_{j-1}^{n} - 2\phi_{j}^{n} + \phi_{j+1}^{n} \right) + \frac{\beta}{(\Delta x)^2} \left( \phi_{j-1}^{n+1} - 2\phi_{j}^{n+1} + \phi_{j+1}^{n+1} \right) \right] + \boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( \phi_{j+1}^{n} - \phi_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( \phi_{j+1}^{n+1} - \phi_{j-1}^{n+1} \right) \right] $

Note that $\beta$=0.5 for Crank-Nicolson, $\beta$=1 for fully implicit, and, $\beta$=0 for fully explicit.

To simplify the notation let's make the substitution,

$ s = D\frac{\Delta t}{(\Delta x)^2} \\ r = \boldsymbol{v}\frac{\Delta t}{2 \Delta x} $

and move the known value $\phi_{j}^{n}$ of the time derivative to the right-hand side,

$ \phi_{j}^{n+1} = \phi_{j}^{n} + s \left( 1-\beta \right) \left( \phi_{j-1}^{n} - 2\phi_{j}^{n} + \phi_{j+1}^{n} \right) + s \beta \left( \phi_{j-1}^{n+1} - 2\phi_{j}^{n+1} + \phi_{j+1}^{n+1} \right) + r \left( 1 - \beta \right) \left( \phi_{j+1}^{n} - \phi_{j-1}^{n} \right) + r \beta \left( \phi_{j+1}^{n+1} - \phi_{j-1}^{n+1} \right) $

Factoring the $\phi$ terms gives,

$ \underbrace{\beta(r - s)\phi_{j-1}^{n+1} + (1 + 2s\beta)\phi_{j}^{n+1} -\beta(s + r)\phi_{j+1}^{n+1}}_{\boldsymbol{A}\cdot\boldsymbol{\phi^{n+1}}} = \underbrace{ (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n}}_{\boldsymbol{M\cdot}\boldsymbol{\phi^n}} $

which we can write in matrix form as $\boldsymbol{A}\cdot\boldsymbol{\phi^{n+1}} = \boldsymbol{M}\cdot\boldsymbol{\phi^{n}}$ where,

$ \boldsymbol{A} = \left( \begin{matrix} 1+2s\beta & -\beta(s + r) & & 0 \\ \beta(r-s) & 1+2s\beta & -\beta (s + r) & \\ & \ddots & \ddots & \ddots \\ & \beta(r-s) & 1+2s\beta & -\beta (s + r) \\ 0 & & \beta(r-s) & 1+2s\beta \\ \end{matrix} \right) $

$ \boldsymbol{M} = \left( \begin{matrix} 1-2s(1-\beta) & (1 - \beta)(s + r) & & 0 \\ (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) & \\ & \ddots & \ddots & \ddots \\ & (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) \\ 0 & & (1 - \beta)(s - r) & 1-2s(1-\beta) \\ \end{matrix} \right) $

Applying Neumann boundary conditions

NB is working through the derivation again I think I have spotted the error. I assumed a fully implicit scheme ($\beta$=1) when writing the finite difference of the boundary condition. If you assume a Crank-Niscolson scheme here the complexity become too great and I could not solve the resulting equations to eliminate the nodes which are outside the domain. However, it would appear possible, there are two equation with two unknowns, but I couldn't manage it. This probably explains the difference between the first and second plots above. I think we can conclude that only the plots with $\beta$=0.5 at the boundary points are valid.

Assuming the flux at the left-hand side is known (assuming a fully implicit form),

$ \frac{\partial\phi_1^{n+1}}{\partial x} = \sigma_L $

Writing this as a centred-difference gives,

$ \frac{\partial\phi_1^{n+1}}{\partial x} \approx \frac{\phi_2^{n+1} - \phi_0^{n+1}}{2\Delta x} = \sigma_L $

therefore, $ \phi_0^{n+1} = \phi_{2}^{n+1} - 2 \Delta x\sigma_L $

Note that this introduces a node $\phi_0^{n+1}$ which is outside the domain of the problem. This node can be eliminated by using a second equation. We can write the $j=1$ node as,

$ \beta(r - s)\phi_0^{n+1} + (1+2s\beta)\phi_1^{n+1} - \beta(s+r)\phi_2^{n+1} = (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n} $

Substituting in the value of $\phi_0^{n+1}$ found from the boundary condition gives the following result for the $j$=1 row,

$ (1+2s\beta)\phi_1^{n+1} - 2s\beta\phi_2^{n+1} = (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n} + 2\beta(r-s)\Delta x\sigma_L $

Performing the same procedure for the final row (at $j$=$J$) yields,

$ -2s\beta\phi_{J-1}^{n+1} + (1+2s\beta)\phi_J^{n+1} = (1-\beta)(s - r)\phi_{J-1}^{n} + (1 - 2s(1-\beta))\phi_{J}^{n} + 2\beta(s+r)\Delta x\sigma_R $

Finally making the boundary rows implicit (setting $\beta$=1) gives,

$ (1+2s)\phi_1^{n+1} - 2s\phi_2^{n+1} = \phi_{j-1}^{n} + 1\phi_{j}^{n} + 2(r-s)\Delta x\sigma_L $

$ -2s\phi_{J-1}^{n+1} + (1+2s)\phi_J^{n+1} = \phi_{J}^{n} + 2(s+r)\Delta x\sigma_R $

Therefore with Neumann boundary conditions we can write the matrix equation, $\boldsymbol{A}\cdot\phi^{n+1} = \boldsymbol{M}\cdot\phi^{n} + \boldsymbol{b_N}$,

where,

$ \boldsymbol{A} = \left( \begin{matrix} 1+2s & -2s & & 0 \\ \beta(r-s) & 1+2s\beta & -\beta (s + r) & \\ & \ddots & \ddots & \ddots \\ & \beta(r-s) & 1+2s\beta & -\beta (s + r) \\ 0 & & -2s & 1+2s \\ \end{matrix} \right) $

$ \boldsymbol{M} = \left( \begin{matrix} 1 & 0 & & 0 \\ (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) & \\ & \ddots & \ddots & \ddots \\ & (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) \\ 0 & & 0 & 1 \\ \end{matrix} \right) $

$ \boldsymbol{b_N} = \left( \begin{matrix} 2 (r - s) \Delta x \sigma_L & 0 & \ldots & 0 & 2 (s + r) \Delta x \sigma_R \end{matrix} \right)^{T} $

My current understanding

• I think the difference between the first and second plots is explained by noting the error outlined above.

• Regarding the conservation of the physical quantity. I believe the cause is that, as pointed out here, the advection equation in the form I have written it doesn't allow propagation in the reverse direction so the wave just passes through even with zero-flux boundary conditions. My initial intuition regarding conservation only applied when advection term is zero (this is solution in plot number 2 where the area is conserved).

Would you agree?