Conservation of Mass in 1D Advection-Diffusion Equation

My long-term goal is to numerically solve the 1D advection-diffusion equation of the form:

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

To start I've set set $v(x)$ instead of $v(x,t)$ and I've been using an implicit finite difference method with various upwind approximations (I've tried first order and also second order linear). While the behaviour of the various solutions appears reasonable (and agrees with Matlab's pdepe solver) I've been unable to get conservation of mass, and in all cases the total mass eventually just reduces to 0 (this doesn't happen with Matlab's pdepe).

I'm confident it's not leaking out of my boundaries (I've tried both periodic BC's and Robin BC's), and I'm beginning to think this is a more fundamental problem with my approach.

My question is, is it reasonable to use an upwind-based FDM to solve such a conservation equation, or am I always going to run into conservation-of-mass problems?

Thanks!

-
If you give more information around your implementation details then somebody might be able to spot where things are going wrong. –  boyfarrell Dec 11 '13 at 23:04
what do you call "mass", an integral of $u$ or something? And why would it be conserved in the presence of diffusion? –  chris Dec 12 '13 at 7:45
en.m.wikipedia.org/wiki/Convection–diffusion_equation this is a conservation equation where a physical quantity (normally "mass") moves under the power of a velocity field (for example gravitational potential) or because of a gradient in concentration. For divergence zero field without and sources or sinks the there is no way to create or destroy the physical quantity so it must remain constant over the domain (boundary conditions permitting). –  boyfarrell Dec 12 '13 at 14:15
So you compute "total mass" as $m\equiv\int_x u dx$ and expect $m$ to be constant in time. And you assume that $v(x,t)$ is constant in space (divergence-free in 1D). –  chris Dec 12 '13 at 14:46
Hey, thanks for your help! I still couldn't manage to get the upwind method to conserve mass (though I got some really interesting results, for instance, setting the diffusivity to 0 and using the compressible flow-field $v(x)=sin(x)$ on $[0,2\pi]$ we reach a steady state where the (numerical) diffusion balances the inward compressible flow, but with both periodic and Robin BC's we lose exactly half the mass in getting there!) In the end I got this to work with a Lax-Friedrich scheme, and I used an implicit method for modelling the diffusion. –  tom Dec 16 '13 at 4:00

Use a first order upwind (for the convection component) and a second order central difference (for the diffusion component). So the end result would be equivalent to discretising the equation,

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

So using the $\theta$-method you will end up with,

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

where $u^n$ and $u^{n+1}$ terms are the present and future time step, respectively. So a fully implicit scheme would be recovered by setting $\theta=1$ and the Crank-Nicolson by setting $\theta=1/2$.

To solve this equation you need to write as a linear system so group in terms of the solution variable $u$ and move the unknowns and known to different sides and you should be able to solve.

For example, the end result would have the form,

$$\underbrace{\boldsymbol{A}}_{A}\cdot\underbrace{\boldsymbol{u^{n+1}}}_{x} = \underbrace{\boldsymbol{M}\cdot\boldsymbol{u^{n}}}_{d}$$

where,

$$\boldsymbol{A}=\begin{pmatrix} 1+2s\theta & -\theta(s + r) & & 0 \\ \theta(r-s) & 1+2s\theta & -\theta (s + r) & \\ & \ddots & \ddots & \ddots \\ & \theta(r-s) & 1+2s\theta & -\theta (s + r) \\ 0 & & \theta(r-s) & 1+2s\theta \\ \end{pmatrix}$$

and

$$\boldsymbol{M}=\begin{pmatrix} 1-2s(1-\theta) & (1 - \theta)(s + r) & & 0 \\ (1 - \theta)(s - r) & 1-2s(1-\theta) & (1 - \theta)(s + r) & \\ & \ddots & \ddots & \ddots \\ & (1 - \theta)(s - r) & 1-2s(1-\theta) & (1 - \theta)(s + r) \\ 0 & & (1 - \theta)(s - r) & 1-2s(1-\theta) \\ \end{pmatrix}$$

Finally $\boldsymbol{u^n}$ and $\boldsymbol{u^{n+1}}$ are just column vectors of the solution variable discretised over the space. And very finally, $s=D\frac{\Delta t}{(\Delta x)^2}$ and $r=\boldsymbol{v}\frac{\Delta t}{2 \Delta x}$.

A side note.

For 1D conservation problems you could also consider the finite volume method, it is just as easy to implement. It has some advantages in the fact that the conservative property is preserved to the discretisation level. You have to be a little be careful because discretisation will be come unstable for advection dominated problems so you need to introduce adaptive-upwinding or exponential fitting. For a step-by-step guide see my notes.

-
Thanks for that, I've actually done a very similar thing, but I'm losing mass. I'm in a very similar situation to yourself here: scicomp.stackexchange.com/questions/5434/… Did you manage to resolve your troubles by changing your BC's? BTW, your notes are fantastic, I'll definitely try FVM sometme –  tom Dec 11 '13 at 5:54
Actually I still don't understand why mass leaves the system in the final figure. The boundary conditions should prevent it, right? When I do the same simulation with the FVM I get conservation, see here scicomp.stackexchange.com/questions/7736/… Also note if you use trapezium rule to keep track of your mass you will see it continually change. With the FVM, where you use discretised the integral equations, you can so a simple sum to count the mass. Notice it stays exactly the same over time. –  boyfarrell Dec 11 '13 at 6:18
@tom figured it out! It's because a Neumann BC applied to the advection-diffusion equation only constrains the diffusion part of the flux; the boundary remains open to the advection component. See here, scicomp.stackexchange.com/questions/10407/… –  boyfarrell Jan 24 '14 at 0:12

A lot of this will depend on the magnitude of $D$ and mesh size $h$. For pure convection, upwinding works well; however, upwinding ends up too diffusive for $D \approx h$. If $h \ll D$, then upwinding should behave less diffusively, though only will $D=0$ will the upwind scheme show minimal numerical diffusivity. As $D$ increases relative to mesh size, an effective finite difference scheme should transition from upwinding to a central difference scheme (which is optimal for the case $v=0$).

Gresho talks a bit about how upwinded methods are generally over-diffusive, and this same idea is what is used to design stabilized schemes for finite elements (i.e. SUPG, etc). In both cases, respecting the balance between $h$, $D$, and $v$ is crucial to the design of the scheme.

-

What is "u" in your advection-diffusion equation? If it represents the mass-fraction of a species then the total mass of that species will likely vary over time. Also, in this case the advection-diffusion equation itself is the continuity equation of that species.

I also believe that the equation should have a negative sign in front of the first term in the parenthesis of the right hand side.

Conservation of total mass should be satisfied from the initial condition unless you plan on solving a pressure-velocity coupling as well.

-