species constraint violated for simplified-SDC + Detonation during advection

[zingale]

This is the same problem as in #1596, but run with the numerical Jacobian instead of the analytic Jacobian.

./Castro1d.gnu.SMPLSDC.ex inputs-det-x.simplified_sdc amrex.fpe_trap_invalid=1 integrator.jacobian=2

Now we run for 32 steps before getting:

amrex::Error::0::Invalid mass fraction in Castro::normalize_species() !!!

The backtrace shows that this is happening after we construct the hydro source and apply it to the state, so the advection itself seems to be leading to the species issues.

This might be because the reaction predictor used in the Riemann problem is doing too much now? Maybe we need to do characteristic tracing on it?

[zingale]

a very crude characeteristic tracing of the passive's source term gets use 2x further in simulation time (but eventually crashes for the same reason)

[zingale]

we make no attempt to ensure that the primitive variable source I_q sums to 0, so maybe we should explicitly normalize (and clip?) the q's that are input to the states?

We do normalize the species fluxes using the CMA procedure, but that won't prevent undershooting.

[maxpkatz]

we make no attempt to ensure that the primitive variable source I_q sums to 0

In the failing run, what do they sum to?

[zingale]

I don't know. I'll find out. In fact, I'll do a PR to make that abort more informative.

[zingale]

also, we might want to have the species non-conservation from advection trigger a retry

[zingale]

here's the problematic species:

(i, j, k) = 315 0 0 , X[16] = -0.01412742191

[zingale]

note: zeroing out Iq makes the problem go away. I don't think it is in the construction of Iq anymore, but now I suspect that it is in the source tracing. When we do a parabolic reconstruction of Iq for the tracing to the interface states, we break the sum{Iq} = 0 that we should have via construction.

So maybe the thing to do for this is just a simple dt/2 addition of the cell-center source to the interface (we can still upwind it, if desired)

[zingale]

note: PR #1601 changed the step where we hit this problem by a bit, but we still encounter it around step 80

[zingale]

now I'm wondering if the species source in the hydro predictor is making the interface q's go outside [0, 1], so maybe we need to clip the interface values?

[zingale]

explicitly clipping the interface values to lie between [0, 1] in trace_ppm doesn't fix this either. So something it's still something about the Iq predictor, but it is strange how it is interacting with the interface states...

[zingale]

actually, I only did the clipping in the normal states, but we might also need to do this after the transverse correction

[zingale]

actually, I am not sure if it makes sense to do this check on SDC state, since we are checking on S_new when it has seen the advection with a prediction of the reaction, but not the final reaction term yet, so the state is kinda incomplete. I don't think it makes sense, as it does with Strang, to check in this case. We should instead wait until the final update to do this check.

[zingale]

This was fixed by Microphysics PR 643, linked above. It ran as above for 2000 steps and reached a time 1.86e-4 s.