Polar Papercuts

[bprather]

Pushing through later INCITE runs and chasing down people's issues, I've learned that some issues I assumed were unique to high resolution and adiabatic index, were in fact universal to the WENO5/face-centered/transmitting scheme we had put together. This PR extends the mitigations I was using at high res back down to all simulations, and fixes some little issues with them.

This PR also fixes an issue with the implicit solver never "un-flagging" previously flagged zones, leading to comically large numbers of failed zones causing dissipation problems. Finally, it backports #133 for the stable branch, which could help dramatically for anyone using excised polar boundary conditions or cell-centered/Flux-CT transport. There's a small chance it fixes even more than that, see the PR for details.

Reconstruction

WENO5 is particularly unstable at the pole, and this led to large, rapid B field increases which would cause the code to stop due to large divB. In particular, if there is a coherent $B_{\phi}$ around the pole (as there often is) but the sign oscillates in X1 (imagine stacked, anti-aligned field loops around the pole), then WENO5 will cause positive feedback, reconstructing fluxes into the already positive loops and out of the already negative loops. Note this also disappears if we use reconstructed versions of all three magnetic field components (i.e., of X1, rather than using the value already present on the face). But reconstructing B leads to other problems (#79).

I don't believe the WENO feedback loop is the source of all polar issues (more on that below) but it's certainly a rake we're leaving lying around to step on. One can attempt to just prevent the prerequisite state from happening (just, never accumulate coherent layers of B in oscillatory directions) but I don't think there's a good heuristic or fix to always prevent these WENO runaways.

Instead, we should use a different reconstruction, and indeed the problem disappears when I do (despite persisting with different floors, reflecting or transmitting poles, FOFC on/off, lower CFL, etc). In particular I've tested the linear options, PPM, and "WENO5-AO at home," the reconstruction used in Phoebus. "WENO5-AO-AH" is usually WENO5Z, but interpolates to a linear reconstruction as the WENO5 smoothness parameters get bad, which surely happens whenever we encounter runaway positive feedback -- overall it's much like an "adaptive order" scheme, but less formal. It behaves pretty similarly to WENO5 even under duress (I tested it as a part of #79, and it only adds some wobbles). So, I've made WENO5-AO-AH the default in the sample mad.par and sane.par, and I recommend using it.

There was a recent issue noticed with this scheme in Phoebus, where the interpolation factor could slightly escape the range [0,1] -- I don't know if KHARMA suffers the same, but I've "fixed" it by just clipping the value to always lie in that range.

Reconnection

The "instant-death" WENO feedback loop is not the only way that the poles can cause problems, however. Smaller explosions seem to occur even under WENO5-AO-AH or other reconstructions, which seem to consist of high temperature and magnetic field buildup, pushing out a wave through the disk. These appear to be due to high magnetic field at the pole, as well, though this is less certain than the above case since explosions can take thousands of steps to develop over broad regions of X1/X3.

What we do know is that regardless of polar boundary conditions, the azimuthal components of velocity and magnetic field, U3 & B3, still build up at polar boundaries -- less in the transmitting case, but still enough to cause problems. Fixing the averaging to fully allow rotation under transmitting boundaries (#126) can't have hurt, but the behavior seems to persist after that fix as well. Physically, this buildup represents tight loops of field and quickly circulating material, which would reconnect or blow away if not for the strained connectivity of the logically Cartesian grid. Thus the idea of forcible reconnection (see #90).

Reconnecting B3 was off by default because it's a new algorithm no one else has tried -- indeed there was a bug (#131) around conserving energy with it enabled (which is fixed in this PR). Now at least a few more pairs of eyes have passed over the code, and it has more or less proven itself necessary if one adopts the other algorithmic changes (i.e., less dissipative face-centered B field, reliable Kastaun inversion which does not inject dissipation through failures). Thus I've also enabled it by default in mad.par and sane.par, and (somewhat reluctantly) recommend it, especially for MADs at adiabatic index 5./3 or if your simulation fails without it.

[vedantdhruv96]

Thanks @bprather! Will give this a try with my modified gravity runs.