Consistent face fields less accurate & stable in KHARMA

[bprather]

tl;dr people using face_ct in KHARMA might consider setting the upcoming option flux/consistent_face_b=false, if they are seeing instabilities peculiar to the magnetic field -- with the understanding that the resulting scheme is not ideal. I hope to better determine how and why this works, and better characterize the potential consequences of this choice, in the future.

Most face-centered transport schemes have a clause which looks something like this (from Stone & Gardiner 2009):

That is, they set e.g. B1 at face F1 to be the actual face-centered split value, rather than one reconstructed from cell centers (which values themselves were reconstructed from faces). This is logical, as it circumvents a whole game of telephone from faces to centers to faces again (though, it shouldn't have an effect on formal convergence).

However, disabling this feature in KHARMA appears to be more stable and more accurate. To be clear, I'm not claiming this is true in general, I'm claiming there's an issue in KHARMA making it true for us. Speculation on that issue later.

Re: stability, here's a blast wave using donor cell reconstruction, LLF fluxes, Kastaun inverter, Gardiner & Stone '07 simple upwinding, and CFL=0.3 of stability limit, with "consistent faces" i.e. true values. I'm plotting Lorentz factor, because it's the hardest thing to get right & good-looking.

And here's the same but with reconstructed face values:

The "wings" of very high velocity disappear! These "wings" are nasty pieces to integrate and bear every mark of the traditional failure modes of HARM & KHARMA: high velocity and internal energy, and checkerboarding at higher Courant numbers (meaning more than half the stability limit).

But wait, that clean blast wave still only hits Gamma=3! Doesn't that mean reconstructing these values is less accurate? Well, the lower velocity appears to be more an issue of the constrained transport scheme & reconstruction than of consistent faces. Here's the same problem with Balsara & Spicer reconstruction, and with the true centered upwinded fluxes from Gardiner & Stone '05. Both problems generate some nasty high velocities when the density and internal energy get low, but they get the point across that the B field can be accurate even when reconstructed.
Balsara & Spicer:

Gardiner & Stone '05 upwinded (not sure of my implementation here, grain of salt):

Note Balsara & Spicer looks exactly the same as using cell-centered Flux-CT, as we might hope:

Okay, what about with reconstruction schemes people actually use? Here's a blast wave with KHARMA's tried & true WENO5 with first-order fallback (FOFC) to Donor Cell, both original and fallback using consistent faces:

And here with reconstructed faces:

And the reference figure from Komissarov (1999) (not quite comparable time, but it gets the point across):

So, is any of this perfect? Clearly not. But barring a full solution to the underlying issue, there seems to be a significant stability and accuracy gain in KHARMA from using reconstructed B field values at faces, rather than replacing them with the "true" values, which is why I wanted to open this combo bug/PSA. I'll be adding an option flux/consistent_face_b to dev, which can be set to false to use the reconstructed values.

The validity of the resulting scheme is hard to gauge, since the values really only differ outside the regime where it's easy to test convergence. I'd be happy to know a) if the reconstructed-faces scheme is clearly bad from a theoretical perspective, or if not, how to verify it better, b) any speculation as to why replacement doesn't work well in KHARMA, when clearly it does in AthenaK, despite using a supposedly identical algorithm for this test (same inversion, B field upwinding, reconstruction, integrator, and floor values, at least).

[bprather]

FAQ (not actually frequent):

1. "Is it an off-by-one when replacing reconstructed face values?": No. I've checked about a million times. You're welcome to check too. Mostly, replacing face values doesn't modify them, with only a few crucial exceptions during the first dozen steps of this problem. The replacements always seem to be between the adjoining cell-centered/reconstructed values (recall we're using donor-cell).
2. "Are you just drawing from old memory?/Is it ordering during a step?": Everything I've shown behaves identically using an RK1 integrator, so likely not memory. The replacement/consistent values are drawn from the same face-centered conserved B field we just used to calculate the cell-centered values (I can add B_CT::UtoP just before calling the kernel and it changes nothing) -- then the replaced values are run through the same flux calculation as if reconstructed. CT is applied just after the flux calculation, using its results to replace any existing updates to the conserved face-centered values, so I'm not sure there's any ordering error in that procedure which would be subtle.
3. "What does AthenaK do?": Always replacement, both in original fluxes and first-order corrections. As with most codes which don't need to support both cell- and face-centered magnetic field options, the reconstructed version is never calculated at all, no code exists to do it. While one can coax vague "wings" out of AthenaK, they are far less pronounced than in KHARMA, and fade with smaller CFL number. Also, the high-order/FOFC versions contain much less trash.
4. "Is it just that reconstructed values are more dissipative?": I don't think so? I can use less-dissipative options for much of this without making the "wings" worse. What does make them worse is using Balsara & Spicer fluxes, or using a higher CFL number where they cause checkerboarding. The worst was using the Noble+ 1Dw inverter with averaging turned on, which just smeared the checkerboarding around to cover the whole domain...
5. "Is it floors?": Also don't think so. I can enable or disable the "normal" floors (as opposed to the new ones applied during Kastaun inversion) and get similar results, even turn on a max Lorentz factor, etc etc. It's irrelevant whether/how each code applies floors after reconstruction, since donor cell is obviously TVD.
6. "Does this mean KHARMA also gets XYZ wrong?": Probably not. I've dived as far as the sound speed calculations and Riemann solver looking for this thing, but the truth is there's not much that would pass modes and run tori but still somehow get these specific zones of this specific problem wrong. Is it causing instability? Yes. Can I see how the old algorithm would also be unstable like this? Yeah, take a look at the lorentz factors that flux-ct version can produce. But incorrect in any meaningful way? I can't see it.
7. "Is AthenaK faster than KHARMA?" Actually they're almost exactly the same on my laptop when running the same algorithm. That suggests that if we see differences with more complex problems, they're likely to be KHARMA implementations or Parthenon peripheral features, not inherent overhead of packing or Parthenon's generality, which is encouraging!
8. "When are FOFC and the Kastaun solver hitting dev?" Soon^TM
9. "How does Flux-CT do with WENO5 reconstruction & FOFC behind it?": Pretty well accuracy-wise, though it doesn't look terribly stable:
   

[bprather]

Well this looks familiar...

Oh, no.

What about with PPM/FOFC?

Oh, no.

Turns out my problem setup for AthenaK was a little easier than the KHARMA one. When I standardize them, I get the same results! Above is the same donor cell/LLF/low CFL blast wave as KHARMA ran, and below I simply switch to PPM/HLLE and enable FOFC. Note these are "velx" ~ U1, but share a scale so you can tell the top image is not reaching full velocity.

So, half of this mystery is solved: using the same algorithm (and problem!) KHARMA and AthenaK agree... and are both wrong, though in wonderfully similar ways. Bug-for-bug compatible! Exciting.

This means I can stop searching for gremlins that make face-B replacement not work "for me." I can also be pretty confident in my implementations of FOFC, Kastaun inversion, and the S&G '07 and G&S '05 "0" CT schemes, which is nice. So, expect those to come down the line pretty quick here when I get Kastaun working for curved spacetimes.

What this doesn't solve is the unreasonable effectiveness of reconstructed B fields: I don't get it, and I don't like it. Frankly, until and even when all the Kastaun/FOFC changes land, my advice is going to be that if you're having trouble with face-centered B, try reconstructing it. I don't know why, I don't know how, but it works.

[bprather]

Looking back at this, you could be forgiven for still seeing a major difference: AthenaK's full reconstruction/FOFC blast wave looks a lot cleaner than KHARMA's version with consistent faces (which recall is the only mode for AthenaK). Instead it looks more like the reconstructed-faces/cleaner version. What's up with that?

The discrepancy is reconstruction. Here's KHARMA using the same PPM reconstruction as used in the AthenaK test, with consistent faces. As far as I can tell this is pretty much identical in the ways that matter (but different CFL number, so don't stare at the noise). Same shape of blast, same artificial wobbles behind it from using PPM.

[bprather]

The saga continues!

Surprising no one, there are some downsides to reconstructed B fields (surely there were reasons no one does this? Here they are). When floors are low, using reconstructed B appears, for lack of a better word... wobbly.
Consistent fields (MAD high-spin):

Reconstructed (same setup & sim time):

Note that when the floors are higher (this instance is particularly restrictive, I don't know if that's necessary) reconstructed fields seem to work out fine. But of course, under these floors almost none of the domain is difficult to integrate, so we may just be seeing the power of WENO.

I think it's safe to say that for tori, reconstructed fields have a steep accuracy penalty wherever there's difficult integration -- which is exactly where they would have been useful, so compared to the other stability mitigations reconstructed fields are a Bad Idea^TM.

But, what about using reconstructed B only when first-order correcting (i.e., at faces using donor-cell reconstruction, where there's a definite stability bonus and potentially still an accuracy bonus)? Stay tuned!

[bprather]

Closing this since we now understand:

1. KHARMA & AthenaK match when using consistent B at faces. Both are prone to instability in the Komissarov blast, relative to reconstructed fields.
2. However, using reconstructed fields across the whole domain is very unstable/inaccurate in tori, and afaict shouldn't be done for serious simulations in GR.
3. Using reconstructed field in specifically FOFC-corrected zones (3-7% of all zones) might be acceptable? Safer to default to not using it, unless necessary for stability.