Temperature_floor appears ineffective with resampled cooling

[YaoguangPei]

Hi all,
I have met some problems of temperature_floor in my custom problem setup.
I am using a custom TurbulentBox problem with resampled cooling. In my input file, I set:

cooling.enabled = 1
cooling.cooling_table_type = resampled
cooling.hdf5_data_file = /public/home/yaoguangpei/quokka/extern/cooling/resampled_ISM_with_iterate_noco.h5
temperature_floor = 10.0

However, after the run finishes, the temperature in the plotfile still appears to go far below the requested floor. When I check the output with yt, I get:

Temperature min = 1e-10
Temperature max = 837.3544098936275

So my concern is that I set temperature_floor = 10.0, but the reported minimum temperature is still much smaller than 10 K.
To investigate this more carefully, I compared three quantities in post-processing:
the stored plotfile temperature,
the temperature reconstructed from total energy minus kinetic energy,
the temperature reconstructed from gasInternalEnergy.
The results are:

T_stored min = 1e-10
T_stored max = 837.3544098936275
T_stored cells below 10 K = 19747821

T_from_Etot min = -0.032378355164263906
T_from_Etot max = 507.0758234055573
T_from_Etot cells below 10 K = 21934303

T_from_Eint min = -0.03237835516426566
T_from_Eint max = 507.0758234055573
T_from_Eint cells below 10 K = 21934303

total active cells counted = 134217728

Since the temperatures reconstructed from total energy and internal energy agree closely, and both fall below 10 K with slightly negative minimum values, this seems to suggest that the evolved thermal state itself may be dropping below the requested temperature floor, rather than this being only a plotfile or derived-temperature output issue.
The derived temperature I set in my sim is coming from the PopIII.

My questions are:
Is my understanding of temperature_floor correct?
Should temperature_floor = 10.0 guarantee that the evolved gas state itself never corresponds to temperatures below 10 K?
If not, what exactly is the intended behavior of temperature_floor?
Does this behavior suggest that I made an error in my problem configuration or in my custom problem implementation?
If this is not expected behavior, what would be the recommended way to fix it?

If you need more info, please contact me.
Thank you very much for your help.

[dosubot[bot]]

This looks like a known class of issues with temperature_floor enforcement in Quokka, particularly when combined with resampled cooling. There are a few things going on:

The core problem is a timing gap in Strang splitting. Quokka uses Strang operator splitting where cooling is applied in two half-steps around the hydro update. EnforceLimits (which applies the temperature floor) is called after each RK hydro stage and in FixupState, but not immediately after the second half-step of cooling [1]. This means cooling can drive the temperature below the floor at the end of a timestep, and that violation persists until the next timestep's FixupState call. Plotfile output can capture this unphysical state.

The dual energy formalism has an asymmetry. When EnforceLimits detects that the auxiliary internal energy is below the floor, it corrects the auxiliary variable but explicitly does not update total energy [2]. This creates a desynchronization between the two energy representations, which SyncDualEnergy only conditionally reconciles [3].

The resampled cooling ODE tolerance is not density-aware. The absolute tolerance for the ODE integrator is set as reltol_floor * E_floor where E_floor is the raw tempFloor_ value rather than a density-dependent energy computed from tempFloor_ [4]. This can allow the integrator to overshoot below the physical temperature floor.

These issues have been documented in previous bug reports: Issue #1136 reported crashes with negative internal energy when using temperature_floor with resampled cooling [5], and Issue #1456 reported negative internal energies despite temperature_floor = 100.0 in a DiskGalaxy simulation [6].

To answer your specific questions:

• Your understanding is correct—temperature_floor should prevent evolved gas temperatures from dropping below the floor, but in practice it currently doesn't fully guarantee this due to the timing gaps described above.
• The slightly negative temperatures you're seeing from both Etot - Ekin and gasInternalEnergy are consistent with cooling driving the internal energy negative between the second cooling half-step and the next FixupState call.

As a potential workaround, you could try adding an explicit EnforceLimits call immediately after each cooling half-step in your problem setup, or manually clamp the internal energy in your problem's source term functions. But a proper fix would likely need to happen upstream in the Strang splitting sequence.

To reply, just mention @dosu.

---

Docs are dead. Just use Dosu.

[BenWibking]

I think the AI is probably right about this. There have been a few reported issues with the temperature floor and the cooling module, but I haven't been able to reproduce them in my own simulations.

@YaoguangPei Are you able to provide your input files for your simulations? This would help to fix the issue.

[YaoguangPei]

Dear Prof. Wibking,
Thank you for your reply.
In my simulation, I use a newly generated cooling table. However, I suspect that the temperature issue is unlikely to be caused primarily by the cooling table itself. For reference, I have included my input file below.

# *****************************************************************
# Problem size and geometry
# *****************************************************************
geometry.prob_lo     =  0.0  0.0  0.0
geometry.prob_hi     =  2.4e19  2.4e19  2.4e19
geometry.is_periodic =  1    1    1

# *****************************************************************
# VERBOSITY
# *****************************************************************
amr.v              = 2       # verbosity in Amr

# *****************************************************************
# Resolution and refinement
# *****************************************************************
amr.n_cell          = 256 256 256
amr.max_level       = 3     # number of levels = max_level + 1
amr.blocking_factor = 32 #grid size must be divisible by this
amr.max_grid_size   = 1024 # at least 128 for GPUs
density_floor       = 1.672621777e-28
derived_vars        = temperature

stop_time = 2*3.1e14
plottime_interval = 3.1e13
max_timesteps = 2000000
checkpointtime_interval = 3.1e14



cooling.enabled = 1
cooling.cooling_table_type = resampled
cooling.hdf5_data_file = "/public/home/yaoguangpei/quokka/extern/cooling/resampled_ISM_with_iterate_noco.h5"
temperature_floor = 10 #K
do_reflux = 1
do_subcycle = 0
hydro.artificial_viscosity_coefficient = 0.1


turbulence.enabled           = 1
turbulence.length            = 1.2e19
turbulence.target_vdisp      = 1e5

turbulence.ampl_factor       = 1.0
turbulence.ampl_auto_adjust  = 1
turbulence.k_driv            = 1.0
turbulence.k_min             = 1.0
turbulence.k_max             = 2.0
turbulence.sol_weight        = 0.5
turbulence.spect_form        = 1

turbulence.random_seed       = -1

turbulence.nsteps_per_t_turb = 10

If needed, I would also be happy to provide more information about the cooling table or other simulation setup details.

Best regards,
Yaoguang Pei

[BenWibking]

I think I have a fix in PR #1829. Can you test the code in that branch and see if it fixes your issue?

[YaoguangPei]

Dear Prof. Wibking,
Thank you for your letter.
I have tested the new cooling code, and it appears to work well. I no longer see any breaks related to temperature_floor.
For example, in one test the temperature statistics are:

    T_min                           = 2.49229854e+01 K
    T_mean_volume                   = 7.41756812e+01 K
    T_mean_mass / density-weighted  = 4.94595904e+01 K
    T_max                           = 4.50778170e+02 K

Thank you very much for taking the time to fix this problem.
Best regards,
Yaoguang Pei

[BenWibking]

That is great to hear. You're welcome. Best wishes in using Quokka!