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

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 :(

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.

• Thank you for explaining. And your intuition is correct; I am solving a system of equations where one of the equation is the velocity (a PDE of the other variables). The system of equations is 1D only, I am planning to use an adaptive 1st order upwind method (can flip between 2nd order central and upwind) maybe with exponential fitting. I am not using flux limiters, but the system is non-linear. Do I need to be "more careful" in this situation? – boyfarrell May 30 '13 at 23:59
• That all depends if you expect shock waves and the like to form, if there's a possibility that the velocity will go below zero in some regions, or if the velocity might become high enough that you will run afoul of the Courant-Friedrichs-Lewy condition at some point. That said, I would try the simple approach first to see if it works, which it may well do. If it's going to fail, it'll do so spectacularly and unambiguously, so I don't think you need to worry about having something wrong slip under your radar. – Daniel Shapero May 31 '13 at 4:47
• Yes I expect the velocity to be only non-zero only in the center of my domain and then rapidly approach zero as one moves away from the centre. I am choosing the time step so that the CFL condition is satisfied (using the max. velocity), the mesh is fixed. What is the criteria for a shock wave? I am not anticipating seeing that (but you never know). – boyfarrell May 31 '13 at 5:56

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.

• Regarding consistency. For a cell centered finite-volume approach. If I have chosen to use (for example) a linear interpolation to estimate the vertex values $u_{j+\frac{1}{2}}$ of the unknown, should I also use the same interpolation method to estimate the value of the coefficient $a(x_{j+\frac{1}{2}})$? Is it OK to "mix" estimates? For example, would it be valid to assume $a(x_{j+\frac{1}{2}}) \rightarrow a(x_{j+1})$ but use interpolation to find $u_{j+\frac{1}{2}}$? – boyfarrell Jun 10 '13 at 9:13
• @boyfarrell It would be ok in the sense that the method remains to be conservative. It does, however, affect the accuracy of the solution. Often times, e.g. in ENO schemes, one approximates the whole flux function and not velocity and solution separately. – mmirzadeh Jun 10 '13 at 20:49

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