MHD Cooling Issue

[YaoguangPei]

I found a cooling-related issue when testing an MHD run with resampled cooling in QUOKKA.

My setup is:

• is_hydro_enabled = true
• is_mhd_enabled = true
• is_radiation_enabled = false
• periodic box
• uniform initial magnetic field

  [
  \mathbf{B}0 = (5,\mu{\rm G},,0,,0)
  ]

The main issue seems to be in the cooling path:

In cooling/ResampledCooling.hpp, cooling updates use
cpp const Real Eint = RadSystem<problem_t>::ComputeEintFromEgas(rho, x1Mom, x2Mom, x3Mom, Egas);

However, in src/radiation/radiation_system.hpp, the implementation is currently

const double Ekin = p_sq / (2.0 * density); 
const double Eint = Etot - Ekin; 

and the inverse function is

cpp const double Etot = Eint + Ekin; 

This is valid for pure hydro, but not for MHD, because magnetic energy is not included.

For MHD, the gas internal energy should instead be recovered as

[
E{\rm int} = E_{\rm tot} - E_{\rm kin} - E_{\rm mag},
\qquad
E_{\rm mag} = \frac{B^2}{8\pi},
]
and similarly the inverse relation should include magnetic energy:
[
E_{\rm tot} = E_{\rm int} + E_{\rm kin} + E_{\rm mag}.
]

Because ResampledCooling.hpp uses ComputeEintFromEgas() to obtain the thermal energy before cooling integration, this means that in MHD runs the cooling module may be evolving an incorrect internal energy.

This likely causes inconsistent thermodynamics in MHD + cooling runs.

In my outputs, I observed:

• magnetic field initialization looks correct
• but thermal quantities become problematic after evolution
• inferred gas pressure becomes negative
• plasma beta becomes negative
• temperature can collapse to extremely small values

So the cooling-related concern is:

ResampledCooling.hpp currently relies on RadSystem<problem_t>::ComputeEintFromEgas(), but that function appears to assume pure hydro and does not subtract magnetic energy. Therefore, the cooling update may be using the wrong gas internal energy in MHD runs.`

This suggests that the current resampled cooling path is not MHD-safe unless internal energy recovery is modified to include magnetic energy.

[dosubot[bot]]

Your analysis is correct — RadSystem<problem_t>::ComputeEintFromEgas() only subtracts kinetic energy [1], while the hydro system's own ComputeInternalEnergy() properly subtracts both kinetic and magnetic energy [2]. The ComputeEgasFromEint() inverse has the same gap [3].

That said, there's currently a compile-time static_assert preventing MHD + Radiation from being enabled simultaneously [4]:

static_assert(!(Physics_Traits<problem_t>::is_mhd_enabled && 
                Physics_Traits<problem_t>::is_radiation_enabled),
              "MHD + Radiation is not supported yet.");

Since ResampledCooling lives under the radiation infrastructure and uses RadSystem energy functions, this assertion should in principle block the combination you're testing. Could you confirm whether you're bypassing this assertion (e.g., with is_radiation_enabled = false but still linking the cooling module), or if the resampled cooling path is invoked independently of the radiation flag?

If cooling is indeed reachable without is_radiation_enabled = true, then this is a genuine live bug. The fix would mirror what HydroSystem::ComputeMagneticEnergy() already does [5] — subtract E_mag in ComputeEintFromEgas and add it in ComputeEgasFromEint. The challenge is threading the face-centered B-field arrays into those call sites, similar to what was done for particle creation in PR #1481.

Notably, the original MHD implementation (PR #900) explicitly flagged a TODO to "audit the code to make sure magnetic energy is subtracted whenever (Etot - KE) is computed" [6], which suggests this particular path was not yet audited.

To reply, just mention @dosu.

---

Share context across your team and agents. Try Dosu.

[YaoguangPei]

@dosu

Yes — in my test, radiation is not enabled.

The issue is not that I am running MHD + radiation.
Rather, I am running MHD + cooling with is_radiation_enabled = false, but the resampled cooling path still calls RadSystem<problem_t>::ComputeEintFromEgas() internally.

So the problem seems to be that cooling is reachable independently of the radiation flag, while still relying on a helper in RadSystem that is not MHD-safe.

[dosubot[bot]]

Thanks for confirming — that makes the issue clear. You're right: cooling is controlled by a separate runtime flag (cooling.enabled) and is completely independent of is_radiation_enabled [1]. Yet the cooling path calls RadSystem<problem_t>::ComputeEintFromEgas() [2], which only subtracts kinetic energy — never magnetic energy. So this is a live bug for any MHD + cooling run.

The static_assert I mentioned earlier only guards the MHD + radiation solver combination, not MHD + cooling. Cooling slips through because it's a standalone Strang-split source term that just happens to use a RadSystem utility function.

The fix likely involves one of two approaches:

1. Make ComputeEintFromEgas MHD-aware — pass in the magnetic field components (or pre-computed E_mag) and subtract them. This is the most direct fix but requires threading B-field data through the call site in ResampledCooling.hpp.

2. Have ResampledCooling use the hydro system's ComputeInternalEnergy() instead — since HydroSystem already handles the three-component decomposition correctly [3]. This would avoid duplicating logic but may require refactoring the cooling interface.

Either way, the reference implementation exists in HydroSystem::ComputeMagneticEnergy() [4], which averages face-centered B to cell centers and computes 0.5 * (Bx² + By² + Bz²). The same pattern would need to be applied wherever ComputeEintFromEgas is called during cooling.

To reply, just mention @dosu.

---

Share context across your team and agents. Try Dosu.

[BenWibking]

Thank you for reporting this bug. I am working on a fix.

[BenWibking]

@YaoguangPei Can you try this fix in #1886? This is now merged into the development branch.