The implosion test: Noise and the Illusion of formal correctness.

by Julius Donnert

This is a fun one ! The implosion test sets up an over density in a corner of a periodic box and let the shocks run through the computational domain. With our new WENO solver, it looks like this at 200^2 zones: 

This test is interesting, because the flow is highly non-linear. That means tiny differences in the algorithm and its implementation become apparent extremely quickly. In particular, the test exposes noise in the computation, which are tiny fluctuations in the fluid quantities around the "correct" solution. This is likely why Sijacki et al. 2012 chose it to compare a traditional SPH code with a moving mesh code. SPH performs pretty badly at the test, because of its inherent sub grid noise. Here is their result at t=0.7: 

Implosion test at t = 0.7 from Sijacki et al 2012.

AREPO performs quite well for a second order code. There is as much structure as in our WENO result, but also more noise. Here is the result from a second order grid code (TVD, courtesy of T. Jones): 

2D implosion test at t=0.7 with a (rather diffusive) second order TVD code (T. Jones, unpublished)

2D implosion test at t=0.7 with a (rather diffusive) second order TVD code (T. Jones, unpublished)

AREPO clearly has more structure at the same resolution (all tests are 200^2 zones). 

So which result is "correct"?  WENO just looks like a sharper version of the TVD result, which is a good sign for the grid codes. Between AREPO and TVD/WENO, there are clear differences  in the structure of the "blob" on the bottom left, e.g. in its arms.

We should re-compute the test at very high resolution and compare again. We expect that all codes should converge to the same "correct" solution, or do we ? 

Consequently,  here is WENO at 1024^2 (excuse the rotation) : 

Implosion test at 1024^2 with WENO.

Implosion test at 1024^2 with WENO.

The high resolution solution of the test shows features from both codes, other features are completely different ! AREPO seems to capture the blobs in the diagonal jet better than WENO, while with interfaces it seems to be the other way around. However, the central RT instability is not captured correctly by both codes. It grows faster at high resolution, that is the curse of non-linear solutions.

To make matters worse, if I had a 1024^2 solution from AREPO to show, it would again look different from the 1024^2 WENO solution above. So there is really no straight forward way to establish formal correctness in this problem ! 

In scientific applications, we never fully resolve the complex problems we are interested in. That means we will always miss some features in the flow. On the other hand, more small structure in the solution does not mean it is more "correct". The fixed-mesh TVD result is a good representation of the high-resolution solution, even though it show less structure than the moving mesh code.

Finally this test was very helpful for us to find an issue in the eigenvector calculation that injected noise at the level of 1e-8 in the density. That is the real meaning of the test. Here is what we were getting with one of the codes: 

WENO5 with numerical noise in the density at the 1e-8 level at t = 0.7.

WENO5 with numerical noise in the density at the 1e-8 level at t = 0.7.

The solution is clearly not correct. The problem here was the compiler injecting a value of 1e-16  inside a square root computation, where it should be 0. Peter Mendygral might write another blog post about it, because this is really about HPC technology. 

I am taking away from this week long exercise that finite volume schemes are incredible sensitive to changes in the implementation - all the seemingly boring code comparison papers are well justified and we should really be careful to not over-interpret our results.