Fix a potential bug in the characteristic projection

bcaddy commented 1 year ago

Prerequisite checklist

[x] My code follows the Athena++ Style Guide
[ ] My change requires a change to the documentation.
[ ] I have updated the documentation in the Wiki accordingly.
[ ] I have added tests to cover my changes.
[x] All new and existing tests passed. (I don't have the FFT library on my system so those tests failed but this change shouldn't impact the fft)

Description

I'm implementing MHD into Cholla and think I might have found a bug in the characteristic projection in Athena++. I'm not positive this is a bug but I'd like some feedback on if it is or not.

In the left and right EigenMatrixDotVector functions bet2 defaults to 0.0. In the equivalent Athena4.2 functions bet2 defaults to 1.0. This is only really an issue when By and Bz are exactly zero but Bx isn't and in that case it can yield incorrect values for the interface states. This case shows up when finding the Y direction interface states in the MHD Einfeldt strong rarefaction where the Vz interface state should be 2.0 but the projection reduces it to 0.0.

Testing and validation

I ran the tests according to the documentation in the Regression Testing page on my macbook, GCC 13 (I think, I don't know how Athena++ finds the right compiler so it might have gotten Apple Clang instead).

felker commented 1 year ago

As mentioned in our email correspondence, it seems likely that this is an undiscovered bug since few people to my knowledge are using characteristic projection, and it would only show up in pathological 1D problems like the Einfeldt strong rarefaction test.

For reference, the Athena (C) implementation is:

https://github.com/PrincetonUniversity/Athena-Cversion/blob/bb328f4ac6fd2c5a2f1881b12fe5e356ed60849b/src/reconstruction/esystem_prim.c#L236-L245

/* Compute beta(s) (eq A17) */

  bt  = sqrt(btsq);
  if (bt == 0.0) {
    bet2 = 1.0;
    bet3 = 0.0;
  } else {
    bet2 = b2/bt;
    bet3 = b3/bt;
  }

From the Athena method paper appendix, Stone et al (2008)

It does not mention the beta=0 conditional, so I am just trying to understand why we handle the case in this way.

felker commented 1 year ago

specifically, why would bet3 be handled differently than bet2?

c-white commented 1 year ago

FYI, this is on my list to look over closely before we make any change to the code.

c-white commented 1 year ago

I agree this is a bug. When $B_y^2 + B_z^2 = 0$, $(\beta_y, \beta_z)$ can be $(1, 0)$ like in Athena 4.2, or $(0, 1)$, or even $(1, 1)$. The eigenvectors are correct for any pair of values. However when they both vanish, the set of eigenvectors becomes degenerate. That is, $\mathbf{L}$ and $\mathbf{R}$ have rank 3 rather than 7, in reference to A12 and A18 of the 2008 Athena paper quoted above.

felker commented 1 year ago

Thanks @bcaddy! I will add a comment in the source code to help prevent this tricky error from being propagated in the future.

See the following references:

Jeffrey and Tanuiti (1964)
Brio and Wu (1988)
Roe and Balsara (1996)

Especially the last one, since Stone (2008) Appendix A3 follows the formalism closely. Brio and Wu basically pick arbitrary values in this case:

Roe (1996), citing BW, do the same thing:

The degenerate case discussed in the Athena 2008 paper appendix: is a subset of the case you encountered, when the transverse fields vanish and the Alfven speed is equal to the sound speed so 2x eigenvalues $u\pm v_A$ each have multiplicity 3--- the "triple umbilic" case V in Roe (1996)

Also, I was refreshing my knowledge of the phase velocities, I came across 2x plots in the Wikipedia article for MHD, e.g. : https://en.wikipedia.org/wiki/Magnetohydrodynamics#/media/File:MHD_wave_mode_1.svg

The legend is definitely wrong--- the green line should be $v{\text{slow}}$ and the orange line should be $v{\text{shear}}$ Alfven wave, no?

felker commented 1 year ago

@c-white I added some comments, but I am not completely satisfied with reconciling Athena's edge case handling with the formulation in Roe (1996). In particular, in the first (the triple umbilic) case, why in the code and in Stone (2008) Appendix A we choose $\alpha_f=1,\alpha_s=0$. Roe (1996) eq 15 provides expansions close to the singularity, which seems to imply that they both approach 1/sqrt(2), but I kinda lost the paper's argument after that.

Also, why are the 3x conditionals expressed as they are? E.g. case IV:

(asq - cssq) <= 0.0)

instead of

std::abs(asq - cssq)  <= epsilon

https://github.com/PrincetonUniversity/athena/blob/07f3bbb291a8cb8f64b8fe8a014bda93788c86eb/src/reconstruct/characteristic.cpp#L92-L110

PrincetonUniversity / athena