Improve the handling of unphysical results in Hydro_HancockPredict()

[hsinhaoHHuang]

Motivations

• Some numerical instabilities were seen in isolated galaxy simulations with feedback and in hydrodynamical cosmological simulations. The unphysical results in the half-step prediction in the MHM scheme may be the source leading to those artificial explosions.
• The original negative pressure check inside Hydro_HancockPredict() was performed on fcPri[f][4], which is the input value before the update. The negative pressure after the update may not be detected, and a floored value is returned at the end.

Goal

• To improve the accuracy and stability of the half-step Hancock prediction in the MHM scheme.

Changes

• Bring back the check of negative pressure after the half-step update.
• Add a velocity check against dh/dt in the unphysical checks.
• Add a three-stage rescue iteration locally when the unphysical results are found.
  • First, recalculate the half-step update with more substeps and with the midpoint integration method.
    • Ideally, it should help and could maintain the desired slope. In practice, however, this is not helpful to rescue the unphysical results most of the time. Therefore, in the default setup, MHM_REPREDICT_ITER_NUM = 2, this will NOT be used.
  • Second, reduce the slope for data reconstruction and then recalculate the half-step update
    • Although it is more diffuse, this usually helps and is not as extreme as the zero-slope method in the third stage.
  • Third, reset to the cell-centered variables (i.e., zero slope for the half-step update).
    • This is the original method to deal with the unphysical results, and it is the safest.

Verification Tests

1. Hydro/Riemann with --flu_scheme=MHM
   • Riemann_Prob = [0, 1, 2, 3, 4, 5, 9] was run by ALWAYS triggering the rescuing methods (separately) regardless of whether the original update is unphysical to test the correctness of the repredict solver. With the default parameters for the rescue methods, the results can still match the analytical solutions, and there is not much difference from the original update.
     • The results by running with more substeps look almost identical to the original update; while running by the reduced slope, small differences can be seen at some discontinuities.
   • Only Riemann_Prob = 9 can produce negative pressure in the default half-step update. The 3-stage rescue methods are triggered for that, and the final results can match the analytical solution as usual.
     • Running by more substeps actually looks slightly worse, while running by reduced slope looks better than the default.
2. Isolated galaxy with SN feedback
   The previous explosions are not found anymore with the new method.
3. Cosmological hydrodynamical simulations
   The previous explosions are gone after the negative pressure is properly checked and a lower slope is applied.

Future Work

• Support SRHD
• Turn the hard-coded macros MHM_REPREDICT_* into runtime options

Notes

• The default parameters of safety factors and the number of substeps are just rule-of-thumb estimations. They are not theoretically based or empirically optimized.
• The rescue methods are kind of experimental. Any suggestions on improving the accuracy, stability, and performance (especially for GPUs) for this PR will be greatly appreciated.