This fix moves the 1/dx from the flux difference in euler.F to the definition of the flux. This only affects the internal energy field, other fields had this done some time ago. This is important for flux correction.
This will change the answer on tests that use dual energy but do not use radiative cooling.
I have a test galaxy simulation running with this change along with a patch that turns on flux correction for the energy fields.
The run has proceeded about ~75% of the way through the simulation - way past the part where other runs have crashed in the past - so it looks like this does indeed fix some stability issues, even in simulations that include cooling.
That said, the timestep is also significantly more variable than in a run where I prevent flux correction from changing the gas temperature. Due to these slowdowns, The CPU time cost of the simulation that includes this PR is about 5 times the simulation that does not flux correct energy.
I think it might make sense to allow both types of simulations to run. For example, if my setup was focused on a situation where the cooling time is much longer than a dynamical time, I might want to flux correct energy since the internal energy is changing slowly, so the flux corrections on the energy field will be more noticeable against the changes due to cooling. For my problem, the cooling is much shorter than the dynamical time, so the flux corrections are much smaller than the corrections due to cooling.
@dcollins4096, what do you think about adding some logic that will prevent the energy fields from being flux corrected (i.e. what the code does now) and in addition prevents the energy fields from being made conservative during flux correction? This behavior could be toggled on and off with a runtime parameter. That would fix the current implementation (do not flux correct, but incorrectly convert to a conservative internal energy despite the fact that flux corrections never happens) and also make it possible to turn on flux correction at runtime.
ngoldbaum, I think that's good going forward. The current PR should always be on, since this is an actual error. I think an additional PR for machinery like you discuss would be a good move, at least as we sort this temperature bug out.
It's not really a matter of speed - more a matter of doing the correct thing. Right now energy isn't flux corrected but still converted to conservative form and then back to specific form during flux correction. This is incorrect.
Followup on this: I've re-run several simulations using this fix, and it seems to work just fine. It has eliminated my "temperature goes to infinity and sim grinds to a halt" problems (in adiabatic cosmological sims lots of grids). This definitely fixes some stability issues. My one current problem is that setting CorrectParentBoundaryFlux to 1 causes temperature spikes, but that is true in sims with and without the fix in this PR. So, I am going to hit 'approve', since it's critical we get this into the codebase.
I have found RadiativeCooling to be a bit touchy (since it tests radiative cooling instabilities and all, differences tend to get exponentially magnified), but if you used the same machine+compilation options, I don't know why it would be any different at all.
@dcollins4096 are you seeing these small differences between the current tip (or close) and with this PR applied? The change to the buffer zones + valgrind cleanup might be a cause if this was comparing with an old-ish gold standard. The test bot didn't find any problems with this PR on the push suite: http://srs.slac.stanford.edu/hudson/job/enzo-dev/38/console
I have one more data point. I recently had an adiabatic simulation, experiencing the dreaded temperature spike (horn). With this bugfix, it no longer happens. Since the 0.5*geslice() is being pushed to its separate issue, I'm going to approve.