How can I check mass conservation when solving the advection equation using an upwind scheme?

Question

My question is how to keep track of the "mass" being advected out of a model domain, for the 1-D advection equation, and an upwind differencing scheme. Following is the background

Consider the continuity equation

$$\frac{\partial u}{\partial t} + \frac{\partial \Phi}{\partial x} = 0$$

Suppose the flux $\Phi$ is given by an advective term

$$\Phi = a u$$

Suppose $a = 0.01x$, and we are considering only $x>0$. Therefore, $a > 0$, so the wave is traveling in the positive $x$ direction. I can discretize the spatial derivative with a first order upwind scheme. $j$ is grid cell centers, and $j-1/2$ is the left edge of the grid cell.

$$\frac{\partial \Phi}{\partial x}|_j \approx \frac{\Phi_{j} - \Phi_{j-1/2}}{\frac{1}{2}\Delta x}$$

$$= \frac{a_{j} u_j - a_{j-1/2} u_{j-1/2}}{\frac{1}{2}\Delta x}$$

I replace $u_{j-1/2}$ with the averge of the grid cell neighbors ($\frac{u_{j-1} + u_{j}}{2}$)

$$= \frac{a_j u_j - a_{j-1/2} \frac{u_{j-1} + u_{j}}{2}}{\frac{1}{2}\Delta x}$$

So, using the method of lines, I've now turned the PDE into a system of ODEs,

$$\frac{\partial u_j}{\partial t} = - \frac{2 a_j u_j - a_{j-1/2} u_j - a_{j-1/2} u_{j-1}}{\Delta x} \tag{1}$$

Now we can consider the left boundary ($j = 1$)

$$\frac{\partial \Phi}{\partial x}|_1 = \frac{\Phi_{1} - \Phi_{1/2}}{\frac{1}{2}\Delta x}$$

Here, $\Phi_{1/2} = \Phi_\text{in}$ is the "flux" entering the left side of the model domain. The left boundary is then given by

$$\frac{\partial u_1}{\partial t} = \frac{-2 a_1 u_1}{\Delta x} + \frac{2 \Phi_\text{in}}{\Delta x} \tag{2}$$

Equations (1) and (2) describe a system of ODEs which can be evolved forward in time. I have done this for an initial condition of $u_j = 0$, and $\Phi = 0.1$ for the domain $x = [0,100]$. The model reaches the following steady state:

Mass is being advected out of the right side of the model domain. How do I figure out this mass flux? The mass flux out should be the same as $\Phi_\text{in} = 0.1$, for flux conservation. Any other relevant tips for understanding/checking mass conservation in these scenarios is appreciated.

I would have thought $\Phi_\text{out} = u_na_n$, where $n$ is the right most grid cell, but this is $\approx 0.105$. So, mass is not quite being conserved, or I'm not computing $\Phi_\text{out}$ properly.

Thanks!

Answer 1

The problem is that my finite difference in the original post does not conserve mass. Mass conserving finite difference schemes exist, but I've found it easier to use and think in terms of finite volumes instead.

For details, see the finite volume textbook by Leveque

Below, I re-do the discretization, using finite volumes.

I can discretize the flux in the following way, considering the fluxes at the edges of each finite volume.

$$\frac{\partial \Phi}{\partial x}|_j \approx \frac{\Phi_{j+1/2} - \Phi_{j-1/2}}{\Delta x}$$

For a first-order upwind scheme, we can approximate the fluxes at the edges of the finite volumes in the following way

$$\Phi_{j+1/2} = a_j u_j$$ $$\Phi_{j-1/2} = a_{j-1} u_{j-1}$$

so,

$$\frac{\partial u_j}{\partial t} = - \frac{a_j u_j - a_{j-1} u_{j-1}}{\Delta x}$$

On the left boundary ($j=1$)

$$\frac{\partial u_1}{\partial t} = - \frac{a_1 u_1}{\Delta x} + \frac{\Phi_{1/2}}{\Delta x}$$

$\Phi_{1/2}$ is the flux at the left boundary. On the right boundary ($j = n$, where $n$ is the total number of finite volumes)

$$\frac{\partial u_n}{\partial t} = - \frac{\Phi_{n+1/2}}{\Delta x} + \frac{a_{n-1} u_{n-1}}{\Delta x}$$

$\Phi_{n+1/2}$ is the flux at the right boundary.

This finite volume upwind scheme conserves mass to machine precision.