# Finite volume piecewise linear 2D advection develops instability

I'm developing a finite volume solver for the simple twodimensional advection equation with constant velocities $u, v$ and constant mesh spaces $\Delta x$: $$\frac{\partial \rho}{\partial t} + u \frac{\partial \rho}{\partial x} + v \frac{\partial \rho}{\partial y} = 0$$

where $u,v > 0$.

According to this lecture, the piecewise linear update scheme in 1D is $$\rho_i^{n+1} = \rho_i^n - \frac{u \Delta t}{\Delta x} ( \rho_i ^n - \rho_{i-1}^n) - \frac{u \Delta t}{\Delta x} \frac{1}{2} (s_i^n - s_{i-1}^n)(\Delta x - u\Delta t)$$

where $s_i^n$ is a slope, that may be chosen as

Centered slope: $s_i^n = \frac{\rho_{i+1} - \rho_{i-1}}{2\Delta x}$ (Fromm's method)

Upwind slope: $s_i^n = \frac{\rho_{i} - \rho_{i-1}}{\Delta x}$ (Beam-Warming method)

Downwind slope: $s_i^n = \frac{\rho_{i+1} - \rho_{i}}{\Delta x}$ (Lax-Wendroff method)

and any of these choices results in second-order accurate methods.

I implemented these methods in 1D, and they work as expected. In the plot below is the 1D advection for $u=1$ with periodic boundaries for a step function at various time steps with 1000 cells. An integer time step means that the curve has gone that many times through the whole domain.

For the 2D method, I evaluate the new value simultaneously in both directions, according to the 2d advection equation as above, while approximating the divergence terms by the piecewise linear scheme as stated above.

If I only let the density advect in one direction, i.e. $u = 0, v=1$ or $u=1, v=0$, I get the expected results, corresponding to the 1D solution. If I choose the timestep to be $\Delta t = [u/\Delta x + v/\Delta y]^{-1}$ as wikipedia suggests (as opposed to $t< [u/\Delta x + v/\Delta y]^{-1}$), I even get perfect advection without diffusion, as it is predicted by the 1d case.

But if I give both $u$ and $v$ nonzero values, e.g. $u = v = \sqrt{2}/2$, at some point, apparently at some point in time, the solution becomes unstable, and I do not know why. This happens for any choice of the slope. Or does it only look unstable to the inexperienced me, but is in fact just the interference of the oscillations produced by the piecewise linear scheme without any slope limiters? (The situation gets a bit better with a minmod or Van Leer slope limiter).

Below you can see the results of a 2d simulation at various time steps for $nx = ny = 200$ for a 2D step function, using the downwind slope. (The timesteps are 0, 0.2, 0.4, 0.6, 0.8, 1.0, 2.0, 5.0, 10.0)

Can anyone tell me how/why this apparent instability occurs, and what to do against it?

Each numerical scheme can have a different stability condition, so you can not use the single one as suggested by the page in Wikipedia. Simply try a smaller time step and if your code is correct you should obtain a stable solution for enough small time step. I am on phone so i can not check it and write it in a mathematical mode, but in your case, I think, your time step must be smaller than the minimal time step of two related one dimensional stability restriction.

By the way if you are showing us a typical example you want to solve, you should check the so called Corner Transport Upwind scheme, e.g. in the book of LeVeque on hyperbolic problems, that has the stability condition you are quoting.

A numerical scheme can be unconditionally unstable, conditionally stable and unconditionally stable.

You would want the later two.

I would recommend that you learn to do Von-Neumann stability analysis for numerical schemes.

You learn it once. And you can find yourself the stability status of any arbitrary governing differential equation for any arbitrary numerical scheme.

Oscillations is a natural result of higher-order approximations near discontinuities/shocks for hyperbolic conservation laws. Recall that the finite-difference approximations you have listed are generally derived using truncated Taylor expansions, which requires a degree of smoothness which is not present in your model problem.

As you have observed, the single-point upwinding does not suffer from oscillations, but has low order. One solution is to use a limiter, which uses the low order approximation where the solution is not sufficiently smooth. Many possible choices are found in the literature, some of which are detailed on wikipedia: https://en.wikipedia.org/wiki/Flux_limiter

• I understand they develop oscillations, but in 1D, this didn't lead to instabilities, while in 2D, it did. My question: How come? – lemdan Jul 21 '18 at 17:42