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?

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? Jul 21 '18 at 17:42

It turns out that the main problem here was that I thought just naively extending the 1D advection $$\rho^{n+1}_{i} = \rho^n_{i} + u\frac{\Delta t}{\Delta x}(\rho^n_{i-1/2} - \rho^n_{i+1/2})$$

to 2D just by adding the $$y$$ term like this:

$$\rho^{n+1}_{i,j} = \rho^n_{i,j} + u\frac{\Delta t}{\Delta x}(\rho^n_{i-1/2,j} - \rho^n_{i+1/2,j}) + v\frac{\Delta t}{\Delta x}(\rho^n_{i, j-1/2} - \rho^n_{i,j+1/2})$$

would be a good idea.

It isn't.

Suppose we have $$u = v = 1$$. Then what one would naturally expect is that the value $$\rho_{i-1, j-1}$$ advects to the position $$\rho_{i, j}$$ over one adequately chosen timestep. But looking at the formula above which I tried to implement, the value of $$\rho_{i-1, j-1}$$ never comes into play. Hence the stencil that I'm using is inadequate. Looking for example at the upper left corner of the square with higher density in the initial conditions, what happens is that we do get diffusion from the right and from below, but do not interact with the cell one place below and to the left, so more material will advect into the cell. Over many time steps, the result is that "stripes" perpendicular to the direction of advection will develop.

Instead, I should have been using an adequate multidimensional technique. Since I wanted to extend my 1D solvers to 2D, a proper extension is dimensional splitting, where you compute an intermediate state first:

$$\rho^{n+1/2}_{i} = \rho^n_{i} + u\frac{\Delta t}{\Delta x}(\rho^n_{i-1/2, j} - \rho^n_{i+1/2, j})$$

called the "$$x$$-sweep", and then we do the "$$y$$-sweep":

$$\rho^{n+1}_{i} = \rho^{n+1/2}_{i} + v\frac{\Delta t}{\Delta x}(\rho^{n+1/2}_{i,j-1/2} - \rho^{n+1/2}_{i,j+1/2})$$

This way, the values "on the corners" are taken care of.

Another really cool thing about this dimensional splitting is that one can show that if we keep the time steps $$\Delta t$$ constant and switch the order of the sweeps every step, the method will become second order accurate in time for basically no extra effort. This is called "Strang splitting".

To demonstrate: Here is what happens if no dimensional/Strang splitting is used: notice the stripes I was talking about.

And this is what happens if the only change is doing proper Strang splitting: • Just a quick comment out of curiosity: I'm not saying that your approach is incorrect but the solution of the pure advection equation that you are trying to solve looks pretty trivial to me. Cause you have: $$\frac{\partial \rho}{\partial t} + \mathbf{u} \cdot \nabla \rho = 0$$ Where $\mathbf{u}$ is the velocity that transport your scalar variable $\rho$. But, you know that total time-differentiation operator is defined as: $$\frac{D}{D t} = \frac{\partial }{\partial t} + \mathbf{u} \cdot \nabla$$ So, your advection equation is written as: $$\frac{D \rho}{D t} = 0$$ Feb 21 '20 at 15:19
• It means based on $\frac{D \rho}{D t} = 0$, scalar parameter of $\rho$ is just constant over each streamline defined as $\mathbf{r} + \mathbf{u} t$. So, I'm not sure what could be found here. And also it is pretty well known that when advection becomes stronger and which its extreme case is in your equation that you only have advection, all sorts of crazy things and instabilities might happen. So my main point here is that: Strang splitting is nice but it won't work for real case applications when you have diffusion or some sort of sources for $\rho$ cause it's too expensive. Feb 21 '20 at 15:24
• That is correct, however the point of the exercise was to get familiar with eulerian finite volume methods, hence the use of this approach. Of course there are better ways to do it, but the idea was to get familiar with the caveats and problems of this method without applying it to more complicated conservation laws. Feb 21 '20 at 15:25
• As far as I have seen, Strang splitting is a general way of solving multidimensional equations based on one dimensional solvers. The number of operations per dimension per time step doesn't increase if programmed properly. However the second order accuracy is only valid for pairwise equal time steps. If that isn't enforced, it falls back to being first order accurate. Feb 21 '20 at 15:28
• So it won't necessarily be a high order method, but at this point I was just happy with it not producing instabilities. If you want to go higher order, you should look into unsplit methods. Feb 21 '20 at 15:29

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.

• Spurious oscillations and instability are not the same thing.
– mfnx
Nov 22 '20 at 23:28

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.