improve conservation for "old" SDC

[zingale]

For the "old" SDC method, where we take a CTU advection term and integrate it along with reactions to get the updated state, the conservation depends on the number of substeps we take in the VODE integration, and can accumulate roundoff error. We should track this error and subtract it off and then explicit add in the advection over the timestep.

It might be as simple as accumulating the integration time, just by adding df/dt = 1 to the system, where f at the end will be the roundoff influence integration time delta t.

[maxpkatz]

I think I see what this is getting at. Would it be fair to say that this is only an issue if the two adjacent zones take a substantially different number of VODE steps? If they took the same integration path, I would expect there to still be conservation to machine precision.

[zingale]

yes, I suspect that is a big part of it. Remember in this case we are using a cell-centered advection term instead of fluxes through the interface.

[maxpkatz]

We're treating it as a cell-centered source term for the purpose of doing the burn, but it's still the case that the A term is composed of fluxes through interfaces, So to distill this issue down and ignore the reactions, we can think of this effect as being similar to what would happen if we used VODE to integrate the A term from t to t+dt, but by taking N steps of length dt / N. In each VODE substep, we're still adding (dt / N) * F_{i+1/2} to zone i+1 and subtracting (dt / N) * F_{i+1/2} from zone i. Any particular VODE update is then still conserved to machine precision, but we're accumulating roundoff error from that N times instead of one time, in the case of CTU.

[maxpkatz]

I think the correct way to capture this is to separately integrate A as we burn in VODE, and then subtract the difference between the integrated A and A * dt from the final update. If that's analogous to what you're saying, then we're in agreement.

[zingale]

it kinda is, but that would double the size of the system. So what is was suggesting is to add d(tau)/dt = 1 to the system. And then integrate this. At the end of the integration, correct as:

(rho X_k) += (dt - tau)*A(rho X_k)

where dt was the size of the step we asked for.

[maxpkatz]

That is definitely more elegant. I'm not certain it does what we want, but it seems plausible enough to be worth trying.

[zingale]

we might be seeing this in the flame problem when you run with inputs.1d.aprox13 with Strang and simplfiied-SDC and look at enuc

[zingale]

this might be related to #1597

also, we store "burn_weights" in the plotfile, so we can look at that to see how big an effect we expect roundoff to be in the region ahead of the burning

[zingale]

it looks like the "preheating" ahead of the flame is there even if I zero out the source predictor to the hydro interface states.

Now my confusion is -- why isn't there any burning ahead of the flame for Strang?

[maxpkatz]

With #1602 we concluded that the preheating ahead of the flame was an artifact of incorrectly including the advection term in "enuc," but this idea is still interesting on the merits, so we'll leave it open for now.

[zingale]

I don't think we want to do this anymore