update to exact integration scheme for radiative cooling

[cmprussell]

Update to exact integration scheme (EIS) for radiative cooling, specifically eq. A7 in Townsend09

Type of PR: modification to existing code, bug fix

Description: One-sentence summary -- The index 'k' in Eq. A7 is not necessarily the same index 'k' in Eq. A5, so this pull request determines the new 'k' for A7, thereby correcting the implementation of exact cooling.

More details -- I believe that EIS was implemented in Phantom based on an EIS subroutine I had passed along years ago from our separate code. There was an issue with the implementation that I passed along, which subsequently made its way into the Phantom implementation. This pull request fixes that issue.

Using Townsend09 as the reference source, Eqs. A5 and A7 both use the index 'k' to refer to specific segments of the piecewise-power-law cooling curve. The prior implementation of EIS in Phantom assumed that the 'k' in each instance was the same, but this is not necessarily true, as explained here.

The 'k' in Eq. A5 depends on the temperature grid, which leads to the expression for A5 on lines 262/264, "y = yk + QrefTgrid(k)/...". Eq. 26 then adds a component to this newly found y, which is done in lines 267/268, "dy = -Qrefdt*T_on_u/Tref, y = y + dy". If dy is big enough (i.e., if the cooling is appreciable), then the old y (line 262/264) from Eq. A5 and the new y from "y = y + dy" (line 268) are in separate portions of the temporal evolution function, so they should have different 'k' indexes. The newly added code in this pull request computes the new index 'k' in the same fashion as it was calculated earlier, though now the requirement for ending the calculation is that Y_k<=ynew (i.e. y = y + dy) < Y{k+1}. Once the correct index 'k' and value yk are found for Eq. A7, the calculation proceeds as before.

Without finding the new 'k' for Eq A7, the algorithm results in constant power-law cooling (according to the initial temperature bin) across the whole temperature range. With this fix, the cooling changes according to variations in the cooling curve as intended.

The one additional change I made is a "Temp = max(...,T_floor)" after "if (Yinv > 0.) then" to ensure that the cooling does not go below the floor temperature. This can happen if Tgrid extends below the floor temperature.

Testing: I made some changes to read in a cooling table (analogous to what was done with cooltable.dat long ago before that was removed) and then ran Phantom with a rho=const, vel=0, T_init=const lattice for a single timestep. By varying the length of that timestep, as well as the initial temperature, I recreated the cooling plots in Townsend09, namely the solid lines in Fig. 1, 2, and 4. The comparison of the old method (k_A5 = k_A7) clearly shows constant power-law cooling, while the newly corrected method shows the variations in cooling as expected.

I'm happy to share these plots, just let me know who all I should share them with.

Did you run the bots? yes I ran "make test" and everything passed. I ran "./testbot.sh" and there were no changes related to the file that I changed -- cooling_solver.f90 -- so I did not make any changes. There were many results from "./testbot.sh" that had nothing to do with the file that I changed, so I did not incorporate these.

Did you update relevant documentation in the docs directory? I don't think this is necessary. If there is a place where this level of detail should go, or if we should mention something like "EIS implementation fixed on ", then I'll do that.

[danieljprice]

This is extremely helpful - thanks Chris !! Please can you post the files to this P-R?

[danieljprice]

what we really need is a unit test that checks and reproduces the plots in Townsend (2009) as part of the test suite. I am working on something similar with Luis Bermudez-Bustamante where we digitised some plots from a paper to use as a reference solution...

[cmprussell]

You're welcome -- I'm glad that this is useful. Here are the verification plots that I mentioned. Note that k refers to the index computed in Eq. A5 and k' refers to the index computed in Eq. A7.

The first two plots compare the two cooling methods -- the old method where k'=k, and the new, corrected method where k' is computed -- to each other for a variety of initial temperatures using the same cooling curve as Townsend 09 (which was from Gnat and Stenberg, 2007). The first plot is linear and the second is log-linear. All computations were done with Phantom as described above. The correct-method plots have more variations than the incorrect-method plots since the corrected method follows all the variations in the cooling curve while the old method was mathematically reduced to constant-power-law cooling.

The next two images over-plot the old and new results on top of Fig. 1 from Townsend 09. The first over-plotted image has small markers for the Phantom calculations while the second over-plotted image has larger markers. Only the top row has over-plotted data. I included the bottom row since the EIS result (solid black line) is duplicated in this row, which makes it easier to see the original Townsend 09 EIS result.