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.
As this changes science runs, it should be extensively tested.
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.
I agree with this. After this gets accepted, perhaps it is worth having an additional flag for what @ngoldbaum suggests, which may provide speed benefits in certain circumstances.
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.
On testing: I've run the push suite, and only two tests don't pass.
Spherical Infall, which seems to only have a very minor answer change. This is expected, as it a.) has DualEnergy on b.) is AMR c.) has RadiativeCooling off.
RadiatingShockWave also gives minor answer changes. This is a little surprising to me, since it has RadiativeCooling on, so should not be affected by this PR.
I want to understand what went wrong with RadiatingShockWave, but other than that I think this is ready to go.
More on testing:
a.) Upon further study, it's very strange that RadiativeCooling is altered, since it's a unigrid run. Is it possible that this test is excessively touchy?
b.) The following is a list of the test that pass the three criteria above. It is now a little odd to me that they actually passed.
I'll look into this a little more and report back when I'm satisfied that I understand everything that it breaks (or doesn't)
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.
It changes the order of some algebra, so if it's really a nonlinear dynamics problem this might make a difference. Also in that case, I reassert that it's not a good test problem.
@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 did it on dffa33b10d8c+cc3f80ebe344, which was my change set plus the tip from just before today's flurry of work. I won't worry about it too much, then.
To update-- The tests all pass Jenkins, so I'm signing off on this. I was expecting more catastrophe, but that's certainly not a good reason to hold up the PR.
Hi all: as a followup to Dave's email, I'm going to merge this PR tomorrow (Friday) afternoon, unless there are strong objections.
I object to merging this without more people signing off on this. I'd also like to hear more discussion on the 0.5 geslice(i,j) comment below.
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.
Would anybody else like to look at this before we merge? Speak now or forever hold your peace...
Final chance to comment on this merge - I will accept tomorrow by noon EST unless somebody chimes in.