Temperature increasing from Radiative Cooling

[mabruzzo]

I've found a weird bug in which the Temperature of a ~7.5e5 K gas at a mass density of 1.3557e-27 g/cm^3 increases with time (the UV background has been turned off). I first encountered this problem while using Enzo-E and have been able to reproduce the error with the python bindings of the latest version. Specifically, the bug can be reproduced with this script.

The script evolves the gas over a single time-step (which is ~10 orders of magnitude shorter than the cooling time). Interestingly, the metal density appears to increase over the timestep by 3 orders of magnitudes.

As an aside: if I switch my densities to be in units of hydrogen masses, the Temperature does decrease over the timestep, but the pressure is 14 orders of magnitudes higher (However, this might be because I am misunderstanding how the units work).

[brittonsmith]

@mabruzzo, the cause of the change in metal density comes from the following line in the ceiling_species function, called by solve_rate_cool:

                  metal(i,j,k) = min(max(metal(i,j,k), tiny),
     &                 0.9_RKIND*d(i,j,k))

This routine is meant to keep species densities from getting too low, but this seems like a bit of a strange choice for what to do with the metal density. Maybe there is something more sensible that can be done there. In general, though, one only hits these limits when the chosen unit system is not keeping the densities in code units near 1.0. This is why you don't see this behavior when setting the density units to the hydrogen mass.

Independently of this, the reason you're seeing the pressure values change so much is that conversion to CGS units is not being done correctly. The values coming out are still in code units and should be multiplied by density_units * velocity_units^2 in the Python functions to be in CGS units. That appears in a few places, but should be quick to fix up.

For the first issue, perhaps we should think about altering the behavior when numerical values approach set limits. Maybe it's better to issue warnings or even errors than to make silent alterations and keep running.

[brittonsmith]

PR #49 fixes the issues relating to the pressure values changing, but the other issue of how better to handle small field values remains. I think I'm leaning toward issuing warnings, but I'm interested what other people think.

[gregbryan]

Thanks all for raising this! I agree that issuing a warning seems better than silently changing the metallicity (particularly in the current formulation). Maybe the warning should be based on the ratio of metal/d, which should be independent of units?

[mabruzzo]

Sounds good to me. Thanks for getting back to me so quickly.

[brittonsmith]

It was agreed at today's meeting that the main thing to do here was to improve the documentation around choosing appropriate unit systems. This can be done for the next release.

However, in thinking about this further, the issue of placing a ceiling on the metal density of 90% of the total density is a bit separate. Is this something we really need to be doing? It is an extreme situation, but maybe not one that requires intervention. @gregbryan, @mabruzzo do you have any thoughts on this?

[gregbryan]

Yeah, placing that ceiling on the metal density seems a bit arbitrary. I suspect it came out of some weird corner case that someone ran into once and then was put into the code. It seems like something we should not be doing, at least not silently -- modifying the metallicity field (violates principle of least surprise). Maybe remove the ceiling and issue a warning or abort if the ratio exceeds 1? That seems likely to be indicative of user error (e.g., incorrect initialization of the metal density field).

[brittonsmith]

Ok, I agree with this. Let's go with this plan. Thanks, Greg!

[brittonsmith]

Resolved in PR #121.