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.


  • $\begingroup$ It seems like your $\Phi_{in}$ should be zero, since $a(x=0)=0$. So I would say that the problem you've stated is mathematically inconsistent. $\endgroup$ Commented Jul 1, 2022 at 6:42
  • $\begingroup$ Changing to $a=0.01x +0.01$ does not really change the result much. Mass is still not conserved $\endgroup$ Commented Jul 1, 2022 at 12:52
  • 1
    $\begingroup$ I think the finite difference scheme I used is just not perfectly mass conserving. I need to be using a finite volume method, which will conserve mass perfectly. $\endgroup$ Commented Jul 1, 2022 at 14:11

1 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}$$


$$\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.

  • $\begingroup$ Non FV schemes still are conservative. You just need to use a conservative method. For example the edge-based discretization which is commonly referred to as node-centered finite volume is actually a finite difference scheme. It maintains conservation because it still uses a flux balance. There are conservative FD, FV, and FE discretizations. $\endgroup$
    – EMP
    Commented Jul 5, 2022 at 17:30
  • $\begingroup$ My understanding is that FV methods are the subset of FD methods that are conserving. $\endgroup$ Commented Jul 5, 2022 at 22:00
  • $\begingroup$ That is not correct. FD is solving the strong form, FV and FEM the weak form. $\endgroup$
    – EMP
    Commented Jul 6, 2022 at 0:30
  • 1
    $\begingroup$ Subset isn't right. But it seems like there is a Venn diagram between FV and FD. You can start with the strong form with FD, or the weak form and apply FV, and end up with the same approximation. $\endgroup$ Commented Jul 6, 2022 at 6:47
  • 1
    $\begingroup$ Philip Roe left a great answer to this subquestion in scicomp.stackexchange.com/a/30277/37438 $\endgroup$
    – IPribec
    Commented Jul 6, 2022 at 8:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.