Problems with setting up shock conditions in SR MHD (concerns dev and pub Athena++, problems run successfully in Athena 4.2)

[ppjanka]

Summary of issue

While trying to set up Rankine-Hugoniot shock conditions in SR MHD, I am encountering odd artifacts that, ultimately, crash the simulations, especially for the newest versions of Athena++.

Steps to reproduce

The general idea of the simulation is as follows:
-- The problem is 2D in special relativity with MHD (configuration settings are given in the athinput files I attach below), in cartesian coordinates.
-- Initial conditions are set up with Rankine-Hugoniot conditions for a vertical stationary shock at x1=10, with magnetic field parallel to the shock interface. For the version I attach below, I hard-set the pre- and post-shock conditions in the input files (with pre-shock Lorentz factor of 10) to limit code complexity. I would be happy to also send a full version of my problem generator, where the RH conditions are calculated on the fly for any given shock parameters.
-- Boundaries are set to be periodic perpendicular to the flow (x2 direction), the inflow (outer x1 bval) is set as a custom user boundary, and the outflow (inner x1) is set as "outflow".
-- The density at the inflow boundary can be modulated with a cosine function, as set in the input file (setting amplitude to 0 removes the modulation). This is intended to "corrugate" or bend the shock interface.

I have created stripped versions of my problem generators, I attach them as intsh2b-mini.cpp (to be used with Athena++) and IntSh2-mini.c (for Athena 4.2), with input files athinputPP.IntSh2b-mini and athinput4p2.IntSh2b-mini respectively, in the following folder:

https://drive.google.com/drive/folders/1dsfkrfF8QyXgrPNdrVM8MEH_Zhwgy3Fj?usp=sharing

The full set of results from my runs with .athdf and .vtk files is also available through my remote backup at:

https://www.dropbox.com/sh/vul99azm4ubyi5x/AADYoOQWJ2W4GS_nOkJalEECa?dl=0

Note that the latter link contains multiple GiB of data, so one may perhaps prefer to download single datasets from Athena*/bin/* folders rather than the entire folder. The "-clean" versions denote the most recent versions of the code at the time I was running, while "-myVer" is the Sept. 15 version I used previously (see commit tags below). Results are located in Athena*/bin/* folders with separate folders for runs with / without inflow density perturbations (pert / noPert), SMR, and Athena 4.2 delimiters copied over to Athena++ ("-Ath4p2delim"). I also include an intsh2.xml file with SR expressions useful with Visit, e.g. to calculate 3-velocities or the flow Lorentz factors (load with Controls -> Expressions -> Load).

Expected outcome

For expected outcome, see, e.g., the density plots from Athena 4.2 simulations (athena4p2-pert-rho.mpg) in the "movies" folder in:

https://drive.google.com/drive/folders/1dsfkrfF8QyXgrPNdrVM8MEH_Zhwgy3Fj?usp=sharing

or the raw vtk files in https://www.dropbox.com/sh/vul99azm4ubyi5x/AADYoOQWJ2W4GS_nOkJalEECa?dl=0 (as described above).

Run without perturbations, final density snaphsot:

Run with perturbations, final density snaphsot:

(Please note that the flow is moving from right to left.)

Actual outcome

What happens:
-- developer version, Sep 15th 2020, commit 714f7bc: I am able to run the initial pure RH conditions indefinitely without problems, but the code starts generating superluminal motions and, subsequently, NaNs once the shock interface starts corrugating (i.e., the density inflow perturbations reach the shock interface).
-- recent public (Jan 6th, 2021, commit 4351e638aa5e75110a4c24cd0e200d6fd7a52a41) and developer version (Jan 15th 2021, commit 0dfb14e): Even pure RH conditions (with perturbations off) will not run, producing instead odd noisy artifacts, as shown below. Note that the density plot below should correspond to the run without perturbations in the "Expected outcome" section.

-- Athena 4.2: Both RH conditions and perturbation / shock corrugation work without problems at low enough CFL (0.025).

For reference, I include the density movies for the current development version with no density perturbations (athenaPPdevCurrent-noPert-rho.mpg), the Sept 15th version with the same (athenaPPmyVer-noPert-rho.mpg), the Sept 15th version with perturbation (athenaPPmyVer-pert-rho.mpg), and Athena 4.2 version with perturbations (athena4p2-pert-rho.mpg), in the "movies" folder in:

https://drive.google.com/drive/folders/1dsfkrfF8QyXgrPNdrVM8MEH_Zhwgy3Fj?usp=sharing

Additional comments and/or proposed solution
I have tried copy-pasting Athena 4.2 delimiters to Athena++, but this did not seem to help with the issues.

Version info

• Athena++ version: developer version, Sep 15th 2020, commit 714f7bc; public version, Jan 6th, 2021, commit 4351e638aa5e75110a4c24cd0e200d6fd7a52a41; developer version, Jan 15th 2021, commit 0dfb14e
• Compiler and version: g++ (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0
• Libraries: libhdf5 1.10.4, openmpi 4.0.3
• Operating system: Ubuntu 20.04.1 LTS
• Hardware and cluster name (if applicable): Intel(R) Core(TM) i7-8750H CPU @ 2.20GHz; Razer Blade 15, 2018 version

[felker]

#316 could definitely be related to the differences in behavior between the September 2020 and current versions of Athena++. I have zero insight for how the SR solvers work in Athena 4.2, but my impression is that there are many possible algorithmic choices that could explain the discrepancies. Especially with flooring variables for this problem.

I confirmed that you are using VL2+(primitive) PLM in both your Athena 4.2 and Athena++ runs. Do you get the same behavior in Athena++ with PPM, RK3?

[c-white]

This will require some digging. That the two versions of Athena++ are different, and that both differ from Athena 4.2, is not surprising. In fact, it is by design. That Athena++ is worse is a little surprising, especially since there should be no nan's ever produced.

That the failures are not happening to all cells in a column is the most surprising. Am I understanding right that without perturbations the problem is independent of y?

[tomidakn]

This looks very serious indeed.

I am not sure this is relevant, but let me point out that this configuration is (at least in non-relativistic simulations) probably unstable against the Carbuncle instability. HLLE is stable against Carbuncle but HLLD and Roe are not. Athena had the H-correction option to suppress it by locally adding viscosity, but we have not implemented it on Athena++ yet. So these configurations may matter.

[ppjanka]

Thank you for your comments and suggestions!

That the failures are not happening to all cells in a column is the most surprising. Am I understanding right that without perturbations the problem is independent of y?

This is correct, yes. The no-perturbation version should be effectively 1-dimensional.

I confirmed that you are using VL2+(primitive) PLM in both your Athena 4.2 and Athena++ runs. Do you get the same behavior in Athena++ with PPM, RK3?

With PPM+RK3, the January versions fail very similarly (in no-perturbation runs), although now full columns of floor values appear before the random noise artifacts:

Adding PPM+RK3 to the Sept. version with perturbations does not seem to affect the result -- the code fails with NaNs as before, shortly after density perturbations reach the shock interface (as with PLM+VL).

I am not sure this is relevant, but let me point out that this configuration is (at least in non-relativistic simulations) probably unstable against the Carbuncle instability. HLLE is stable against Carbuncle but HLLD and Roe are not. Athena had the H-correction option to suppress it by locally adding viscosity, but we have not implemented it on Athena++ yet. So these configurations may matter.

I have now also tried using HLLE (with the original PLM+VL setup). In the January versions of Athena++, this removes the noisy artifacts, with the code now conserving 1-dimensional character of the problem without density perturbations. However, floor regions still appear (now as columns) and the Rankine-Hugoniot conditions do not remain stationary:

However, using HLLE does indeed stabilize the Sept. 15th version (provided that low enough CFL is used), which I can now run with perturbations indefinitely! (Thank you Kengo!)

Athena had the H-correction option to suppress it by locally adding viscosity, but we have not implemented it on Athena++ yet.

I have not deliberately turned on the H-correction option in Athena 4.2 (unless it is used by default?), and I have been using HLLD there.. However, I presume that the solver might be slightly more diffusive there, which perhaps prevented the instability from affecting the simulations?

[felker]

H-correction is off by default in Athena 4.2 and only relevant when using the Roe Riemann solver, I believe: https://princetonuniversity.github.io/Athena-Cversion/AthenaDocsUGH.html
Have you used it with HLLD in Athena 4.2 @tomidakn ?

#-------------------------------------------------------------------------------
# ALGORITHM FEATURE: turn on H-correction in multidimensional integrators
#   --enable-h-correction
# Note that H-correction only works with Roe flux.  Prints error message if
# H-correction is enabled with any other flux.

AC_SUBST(H_CORRECTION_MODE)
AC_ARG_ENABLE(h-correction,
        [--enable-h-correction  turn on H-correction],
        ok=$enableval, ok=no)
if test "$ok" = "yes"; then
if test "$with_flux" = "roe"; then
  H_CORRECTION_MODE="H_CORRECTION"
  H_CORRECTION_MODE_USER="ON"
else
  AC_MSG_ERROR([H-correction only works with Roe flux])
fi
else
  H_CORRECTION_MODE="NO_H_CORRECTION"
  H_CORRECTION_MODE_USER="OFF"
fi

Nothing has changed in the past ~2 years with Athena++ that warrants retesting this, does it @c-white ?

The only significant concern I still see from this issue is that the switch from hlle_mhd.cpp to hlld_mhd.cpp introduces asymmetries along y. NaNs only appeared in the September 2020 version, right @ppjanka ?

[c-white]

Nothing should have changed, @felker

[ppjanka]

The only significant concern I still see from this issue is that the switch from hlle_mhd.cpp to hlld_mhd.cpp introduces asymmetries along y. NaNs only appeared in the September 2020 version, right @ppjanka ?

As far as I remember yes, only the Sept 2020 version would fail with NaNs, while the Jan 2021 version would only generate floor artifacts. However, only the Sept 2020 version could be completely fixed by using HLLE. In the Jan 2021 version, while the assymetries along y would disappear when using HLLE, vertical stripes of floor values would still appear (see the second plot in my ).

Note: In the end, I have opted to use Athena 4.2 for my problem -- I was planning to use test particles, which were already implemented there (I only needed to ). So I have not looked into this problem since early 2021.

[tomidakn]

Sorry for not responding, I've been really busy these days.

@ppjanka , recently I implemented LHLLD which is a newly proposed modification of HLLD, and it is stable against Carbuncle. Internally, it uses slightly different method to detect a shock than the H-correction in Athena, but it should work more or less similarly. I'd appreciate if you could try it if you have time and are interested.

[tomidakn]

Also, the latest change (#457) in the HLLD solver regarding the degeneracy detection might be relevant to this issue.

[ppjanka]

I have just tested the most recent version of the Athena++ developer version (commit e34cd25; Oct 14th, 2022), with the setup identical to original (see my first post), and the following configurations:

• vl2 + plm + hlld,
• vl2 + plm + hlle,
• vl2 + ppm + hlld,
• vl2 + ppm + hlle,
• rk3 + ppm + hlld,
• rk3 + ppm + hlle.

[ @tomidakn : I have tried using LHLLD as well, but there does not seem to be a relativistic version available in the master branch -- should I use a specific branch / version? ]

Unfortunately, all of the configurations I tested continue to fail on the test case with no perturbations (pure 1D shock condition). As a remider, the expected behavior is for the initial condition to remain unchanged:

Instead, in the final frame of all the cases I run, I see vertical stripes of alternating density:

This is better than the previous behavior of HLLD -- the 1D symmetry is now conserved (no noise artifacts). However, the code still does not seem to be able to hold the expected solution...

[tomidakn]

Oh, I'm so sorry I forgot it's a relativistic problem. LHLLD is implemented only for non-relativistic.

[ppjanka]

No problem, @tomidakn. Just wanted to make sure that I wasn't missing something :)