Test errors caused by oscillatory solutions

[conradtchan]

There are small oscillations to the solution near discontinuities, noticeable in the sod shock test problem. Although the amplitude of the oscillations is bound and does not affect the overall quality of the solution, they are extremely sensitive to small perturbations. For example, changing the courant factor by 1.e-15 (machine precision) completely changes the shape of these oscillations without affecting the amplitude.

Here, I compare courant=0.9d0 with courant=0.9d0 - 1.d-15:

Shock tube solution overlaid with error:

Zoomed in on oscillatory part of solution:

For testing purposes, this is problematic because when the oscillations are out of phase to the oscillations in the reference solution, errors of order 1.e-2 appear, which dominate any errors that may have been caused by code changes, hiding potential issues. I don't think we should simply adjust the test tolerances to so that the test passes, because 1.e-2 would let through many real issues.

In practice, these errors can be introduced when restarting from an output file that was written at a manually specified timestep, since a shorter timestep would have been taken to land at the target timestep exactly. These errors can also be introduced with MPI, where the order of operations may vary slightly compared to the serial code.

These oscillations are not present when no slope limiter is used (flux_limiter='flat'), and errors in the solution are only of order 1.e-15, as expected.

Is it worth considering:

• Trying different slope limiters?
• Testing against the analytic solution, when it is known?

[ryosuke-hirai]

This is a well-known problem of spatially higher-order reconstruction. The choice of slope limiter does not improve things as long as it is applied on arbitrary variables. At the moment, I do the reconstruction and apply the flux limiters on conserved quantities like density, momentum and total energy. The proper way to do reconstruction is to use the so-called "characteristic variables", which are the eigenfunctions to the linearized Euler equations. But it does involve some matrix inversion to go from characteristic to primitive variables, so it sacrifices quite a bit of computation time.

I know there is a cheap and simple solution that was presented in this paper. I do want to implement this eventually but haven't got around to do it yet. It will involve heavily rewriting interpolation.f90 and numflux.f90 so I'm not sure if it is feasible within the duration of the current MAP project.

So for now, I think comparing against the analytic solution will be the better option.

[conradtchan]

An alternative is we can run 2 versions of the test: with and without the slope limiter. I'm guessing you never use flat reconstruction in real simulations, but for the test it will allow us to set a very tight tolerance to pick up regressions in the code (except for regressions in the slope limiter subroutine 😄)

[ryosuke-hirai]

That would also work :) I agree it would be better to start off with flat reconstruction for MPI testing before we further complicate the problem with slope limiters!