How should non-constant coefficients be treated with finite-volume first order upwind scheme?

Question

Starting with the advection equation in conservation form.

$$ u_t = (a(x)u)_x $$

where $a(x)$ is a velocity which depend on space, and $u$ is a concentration of a species which is conserved.

Discretising the flux (where the flux $f=a(x)u$, is defined on the edges of the cells between mesh points) gives, $$ u_t = \frac{1}{h}\left( f_{j-{\frac{1}{2}}} - f_{j+{\frac{1}{2}}} \right) $$

Using a first order upwind we approximate the fluxes as,

$$ f_{j-{\frac{1}{2}}} = a(x_{j-\frac{1}{2}})u_{j-1} \\ f_{j+{\frac{1}{2}}} = a(x_{j+\frac{1}{2}})u_{j} $$ Which gives, $$ u_t = \frac{1}{h} \left( a(x_{j-\frac{1}{2}})u_{j-1} - a(x_{j+\frac{1}{2}})u_{j} \right) $$

If $a(x)$ was constant then this will reduce to the familiar upwind scheme i.e., $ u_t = \frac{a}{h} \left( u_{j-1} - u_{j} \right) $.

My question is, how can we treat the non-constant coefficients of the advection equation? The velocity is defined at the cell centers, so a simple approach would be the following,

$$ a(x_{j-\frac{1}{2}}) \rightarrow a(x_{j-1}) \\ a(x_{j+\frac{1}{2}}) \rightarrow a(x_{j}) $$

This is my preferred approach because it is very simple to implement.

However, we could also use (I am guessing) an averaging scheme to define the velocity at the cell edges, $$ a(x_{j-\frac{1}{2}}) \rightarrow \frac{1}{2}a(x_{j-1} ) + \frac{1}{2}a(x_{j} ) \\ a(x_{j+\frac{1}{2}}) \rightarrow \frac{1}{2}a(x_{j} ) + \frac{1}{2}a(x_{j+1} ) \\ $$

In LeVeque's book he says,

So far we have assumed that the variable velocity $a(x)$ is specified by a constant value $a_j$ within the j-th grid cell. In some cases it is more natural to instead assume that a velocity $a_{j−{\frac{1}{2}}}$ is specified at each cell interface.

But he doesn't really elaborate too much after that. What is a common approach?

I am solving a conservation problem (I am using the advection equation as a continuity equation) so I want to make sure that after applying discretisation that the conservation property is preserved. I would like to avoid any hidden surprises regarding these variable coefficients! Does anybody have some general comments and guidance?

---

Update There are two really good answers below and I could only pick one :(

Answer 1

Depending on what kind of system you're looking at, it may be more convenient to consider the velocity $a$ as piecewise-constant within each cell, or that it's defined at the cell interfaces. For example, in meteorology, staggered grids are quite common, where pressure might be defined inside cells and velocity at cell interfaces. You could just as easily think of the velocity as defined within the cells. All told: the choice of representation should not affect the convergence of your method*, provided your discretization is stable and consistent.

What is most important (and you've already touched on this in your question) is that the discretized system is still conservative. Provided that your scheme can be written in the form

$ \frac{\partial u_j}{\partial t} = F_{j-\frac{1}{2}}(u_{j-1},u_j) - F_{j+\frac{1}{2}}(u_j,u_{j+1}) $

then it should be conservative, since

$ \frac{\partial}{\partial t}\int u dx = \sum_j\frac{\partial u_j}{\partial t}\delta x = \sum_j( F_{j-\frac{1}{2}}- F_{j+\frac{1}{2}} )\delta x = (F_{-\frac{1}{2}}-F_{N+\frac{1}{2}})\delta x $.

Your simple approach should work fine, as will averaging the velocity between cells to define it on cell interfaces, provided the velocity is always positive. Moreover, I don't think that averaging will net you any higher accuracy, so you are right to prefer the simple way.

If you are also solving for the velocity and you have a system of equations, you may well need to be more careful. Likewise, if you're solving a nonlinear hyperbolic PDE and using flux limiters, you have to be yet more cautious.

*However, for a system of hyperbolic PDEs, using staggered grids can substantially ameliorate artificial dispersion/diffusion. If you want to know more, look up Arakawa C-grids or check out chapter 4 of this book.

Answer 2

Others have said it all, but I just wanted to add a simple, yet sometimes subtle, point. Your upwind discretization remains conservative as long as you use a consistent interpolation of $a(x)$ on the cell boundaries.

What I mean by consistent is that the only condition that the interpolation needs to satisfy is

$$ a_{i+1/2}^+ = a_{i+1/2}^- $$

In other words, as long as your interpolation method is continuous across cell boundaries, your discretization is guaranteed to remain conservative.

This may not seem like a big issue here in 1D (and it should not) but can cause issues at the coarse-fine interfaces on multi-level AMR grids.

Answer 3

You can use any type of interpolation to determine $a(x_{j-\frac{1}{2}})$, and the method will remain conservative.

To see why this is so, consider that the analytic definition of conservative is that

$$\frac{\partial}{\partial t}\int_D u(x)\, dx = \int_{\partial D} a(x)u(x) dS,$$

where $D$ is the problem domain. This says that the change in the conserved quantity is equal to the flux at the boundary, and follows immediately from integrating the conservation law.

If our discretization is of the form

$$u_t(x_j) = \frac{1}{h}\left(a(x_{j-\frac{1}{2}})u_{j-\frac{1}{2}} - a(x_{j+\frac{1}{2}})u_{j+\frac{1}{2}}\right)$$

where $x_1,\ldots, x_n$ are our grid points, $D = [c,d]$, and $c = x_{\frac{1}{2}}$, $d = x_{n+\frac{1}{2}}$, then the equivalent discrete statement of conservation is

$$\frac{1}{h}\sum_{j=1}^n\left(a(x_{j-\frac{1}{2}})u_{j-\frac{1}{2}} - a(x_{j+\frac{1}{2}})u_{j+\frac{1}{2}}\right) = a(x_{\frac{1}{2}})u_{\frac{1}{2}} - a(x_{n+\frac{1}{2}})u_{n+\frac{1}{2}},$$

and this can easily be observed to hold by expanding the sum on the left-hand side. Note that, in your case with upwinding, $u_{j-\frac{1}{2}} = u_{j-1}$ and $u_{j+\frac{1}{2}} = u_j$, though it should be noted that this scheme is only an appropriate upwind scheme if $a(x)u$ is always positive.

For higher-order methods, provided that $a(x)$ is smooth, one can simply fit a polynomial to the points $a(x_{j-r}), \ldots, a(x_{j+s})$, and evaluate the polynomial at $a(x_{j-\frac{1}{2}})$.