consider changing dtnuc_e to use average de/dt over timestep instead of instantaneous

[zingale]

Currently we use the instantaneous energy release from the network in determining the burning limited timestep. This could be problematic if we are in a region where we are trying to come to equilibrium, so the sign of the energy generation term changing rapidly to get to some equilibrium. In this case, it might be better to use the integral of the energy release over the timestep in the energy limiter.

[maxpkatz]

A note on the history here: when this functionality was first introduced we had several ways of estimating de/dt:

https://github.com/AMReX-Astro/Castro/blob/e278d91f1af5c067d2c98fbfa902d5f3d1c6f9ea/Source/driver/timestep.F90#L216-L248

The versions other than "instantaneously call the RHS" were ditched in #272. I don't remember the exact reasoning but I recall that it included (1) no one was using the other modes and there was no particular evidence that one was better than the others, (2) I think there was concern about how the method that depends on the burner data would mesh with SDC, and (3) I think I may have been motivated by wanting to eventually remove the Reactions_Type data.

[maxpkatz]

In order to test out different ideas, I think we would need to concretely identify a case where a zone burns to NSE in a single timestep but where it is not caught by either the start of timestep limiter (because the burning was still fairly quiescent at the start) or the end of timestep limiter (because the burning to NSE had mostly completed already). Such a thing is easy to imagine but having the actual test case would be helpful in understanding our options.

[zingale]

one thing that would be useful would be to make a plot of the average de/dt over a timestep and another with the instantaneous de/dt -- we can hack the storage in Reactions_Type to change this just for making a plotfile.

[harpolea]

This is instantaneous and average values of de/dt (i.e. the value from ydot(net_ienuc) vs the one from burn_state.e / dt), for the mean positive value and the max(abs(dedt)). This is for the rotating star problem, run with NSE and Strang for 15 time steps