bprather
commented
6 months ago FAQ (not actually frequent):
- "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).
- "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.
- "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.
- "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...
- "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.
- "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.
- "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!
- "When are FOFC and the Kastaun solver hitting
dev?" Soon^TM - "How does Flux-CT do with WENO5 reconstruction & FOFC behind it?": Pretty well accuracy-wise, though it doesn't look terribly stable:
tl;dr people using
face_ctin KHARMA might consider setting the upcoming optionflux/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_bto 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).