Dust gas-energy Jacobian misses heat capacity

## Summary The decoupled dust gas-energy Newton solve uses `dCooling/dT` directly in a residual whose unknown is gas internal energy, missing the conversion by gas heat capacity. ## Severity `High` ## Affected File `src/radiation/radiation_dust_system.hpp` ## Affected Function / Symbol `SolveGasDustRadiationEnergyExchange(...)`; `SolveGasDustRadiationEnergyExchangeWithPE(...)` ## Audit Metadata - Source log: `issues/likely_real/dust-decoupled-gas-energy-jacobian-missing-cv.md` - Finding tags: correctness/numerics ## Proposed Patch - Convert `dCooling/dT` to `dCooling/dE` by dividing by `ComputeEintTempDerivative(...)` in the gas-energy Jacobian. ## Why This Is a Bug `SolveGasDustRadiationEnergyExchange` and `SolveGasDustRadiationEnergyExchangeWithPE` use `BackwardEulerOneVariable` to solve for `Egas_guess` in the decoupled dust model. Their residual functions are written in terms of gas internal energy, but the Jacobian lambdas add `sum(DefineNetCoolingRateTempDerivative(...))` directly: ```cpp return 1.0 + sum(d_cooling_d_Tgas_); ``` ## Complete Code Patch ```diff diff --git a/src/radiation/radiation_dust_system.hpp b/src/radiation/radiation_dust_system.hpp --- a/src/radiation/radiation_dust_system.hpp +++ b/src/radiation/radiation_dust_system.hpp @@ auto jac = [=](double Egas_) -> double { const double T_gas_ = quokka::EOS<problem_t>::ComputeTgasFromEint(rho, Egas_, massScalars); const auto d_cooling_d_Tgas_ = DefineNetCoolingRateTempDerivative(T_gas_, H_num_den) * dt; - return 1.0 + sum(d_cooling_d_Tgas_); + const double c_v_ = quokka::EOS<problem_t>::ComputeEintTempDerivative(rho, T_gas_, massScalars); + return 1.0 + sum(d_cooling_d_Tgas_) / c_v_; }; @@ auto jac = [=](double Egas_) -> double { const double T_gas_ = quokka::EOS<problem_t>::ComputeTgasFromEint(rho, Egas_, massScalars); const auto d_cooling_d_Tgas_ = DefineNetCoolingRateTempDerivative(T_gas_, H_num_den) * dt; - return 1.0 + sum(d_cooling_d_Tgas_); + const double c_v_ = quokka::EOS<problem_t>::ComputeEintTempDerivative(rho, T_gas_, massScalars); + return 1.0 + sum(d_cooling_d_Tgas_) / c_v_; }; ```