Conservation of a physical quantity when using Neumann boundary conditions applied to the advection-diffusion equation

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?

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).

• Even with Neumann zero-flux boundary conditions $\frac{\partial \phi}{\partial x} = 0$ the mass can still leave the system, this is because the correct boundary conditions in this case are Robin boundary conditions in which the total flux is specified $j = D\frac{\partial \phi}{\partial x} + \boldsymbol{v}\phi = 0$. Moreover the Neunmann condition specifies that mass cannot leave the domain via diffusion, it says nothing about advection. In essence what we have hear are closed boundary conditions to diffusion and open boundary conditions to advection. For more information see the answer here, Implementation of gradient zero boundary conditon in advection-diffusion equation.

Would you agree?

• It seems that the boundary conditions are not implemented correctly. Could you show us how you've imposed the boundary conditions? – David Ketcheson Mar 5 '13 at 6:46
• OK I updated with the implementation and I think I spotted the error regarding applying $\beta$=0.5 at the boundary rows only. I have updated my "current understanding" at the bottom of the question. Do you have any comment? – boyfarrell Mar 5 '13 at 10:49
• So...what does the discretization on the boundaries look like in case of the Robin boundaries? You've shown it for the Neumann boundaries, but not the Robin boundaries. – user7513 Mar 3 '14 at 16:17

I think that one of your problems is that (as you observed in your comments) Neumann conditions are not the conditions you are looking for, in the sense that they do not imply the conservation of your quantity. To find the correct condition, rewrite your PDE as

$$\frac{\partial \phi}{\partial t} = \frac{\partial}{\partial x}\left( D\frac{\partial \phi}{\partial x} + v \phi \right) + S(x,t) .$$

Now, the term that appears in parentheses, $D\frac{\partial \phi}{\partial x} + v \phi = 0$ is the total flux and this is the quantity that you must put to zero on the boundaries to conserve $\phi$. (I have added $S(x,t)$ for the sake of generality and for your comments.) The boundary conditions that you have to impose are then (supposing your space domain is $(-10,10)$)

$$D\frac{\partial \phi}{\partial x}(-10) + v \phi(-10) = 0$$

for the left side and

$$D\frac{\partial \phi}{\partial x}(10) + v \phi(10) = 0$$

for the right side. These are the so-called Robin boundary condition (note that Wikipedia explicitly says these are the insulating conditions for advection-diffusion equations).

If you set up these boundary conditions, you get the conservation properties that you were looking for. Indeed, integrating over the space domain, we have

$$\int \frac{\partial \phi}{\partial t} dx = \int \frac{\partial}{\partial x} \left( D \frac{\partial \phi}{\partial x} + v \phi \right) dx + \int S(x,t) dx$$

Using integration by parts on the right hand side, we have

$$\int \frac{\partial \phi}{\partial t} dx = \left( D \frac{\partial \phi}{\partial x} + v \phi \right)(10) - \left( D \frac{\partial \phi}{\partial x} + v \phi \right)(-10) + \int S(x,t) dx$$

Now, the two central terms vanish thanks to the boundary conditions. Integrating in time, we obtain

$$\int_0^T \int \frac{\partial \phi}{\partial t} dx dt = \int_0^T \int S(x,t) dx dt$$

$$\int \phi(x,T) dx - \int \phi(x,0) dx = \int_0^T \int S(x,t) dx$$

This shows that the domain is insulated from the exterior. In particular, if $S=0$, we get the conservation of $\phi$.

• I realise now why the it only worked when $\boldsymbol{v}$=0; because this would imply conservation following your approach above. What would be the consequence of using this boundary condition on the above, would the wave reflect? I thought this wouldn't be possible because there is nothing in the equation that would give me negative velocity? – boyfarrell Mar 5 '13 at 12:52
• The best way to know is probably to try! But if this behaves correctly (and IMO it does), you should see a certain amount of $\phi$ that starts accumulating on the left side of the domain: the advection pushes $\phi$ in that direction but the boundary is closed. The accumulation stops when the diffusion is sufficiently large to balance it. So no, there should be no reflected wave. – Dr_Sam Mar 5 '13 at 13:08
• @DrSam Just a question regarding the implementation. I understand your point about how to make the quantity zero on the left hand side. But when you say "on the right just a boundary term" what does that mean? I thought that boundary conditions should be either Neumann or Dirichlet (or a mix of both)? – boyfarrell Mar 6 '13 at 9:31
• @boyfarrell The left/right in the answer was refering to a derivation of the correct boundary conditions, not the way it is implemented (edited for clarity). Robin conditions are classical conditions, although there are less known than Dirichlet and Neumann. – Dr_Sam Mar 6 '13 at 10:16
• So regarding the implementation, do you think I should derive Robin boundary conditions for both boundaries? Also, if the equation has a reaction term (e.g. $$\frac{\partial \phi}{\partial t} = \frac{\partial}{\partial x}\left( D\frac{\partial \phi}{\partial x} + v \phi \right) + S(x,t)$$ does the boundary condition need to also include this term? – boyfarrell Mar 6 '13 at 11:11