[cpinte]

Type of PR:
Bug fix / new physics

Description:
Fixes the stochastic heating of non-equilibrium grains (PAHs, nRE) so that it
conserves energy without empirical factors. MCFOST now passes the Camps et al.
(2015) benchmark for all 51 cases: emitted/absorbed power 0.97-1.00 per grain
size, emissivity within 0.88-1.05 of DustEM from 3 to 1000 um (the other codes of
the paper: 0.48-1.42).

User impact: PAH and nRE emission changes in existing models. Mid-IR bands are
lower where energy was not conserved (ref3.0_multi: PAH emission was 2.2x the
absorbed power; SED 0.56-0.89x the previous one at 6-13 um, unchanged elsewhere),
and 25-50 A PAHs are now correctly treated as stochastic.

Physics

• Empirical factors removed (1.3e-2, 7e-8, the 0.5 on the heat capacity, the P(T)
  rescaling). They compensated the errors below and a benchmark field 4 pi too
  weak since aedccef (2017).
• Enthalpy integrated from the heat capacity instead of C(T) T.
• Heating rates integrated over each bin's energy width (Draine & Li 2001), top
  bin collecting all higher transitions: energy conserved on any wavelength grid.
• P(T) normalisation no longer overwrites the lowest bin.
• Cooling and absorption times: C_abs at the photon wavelength, consistent units.
• DustEM heat capacity converted to J/K (was 1e7 too large).

Bug fixes

• Lost updates of the energy absorbed by nRE grains with OpenMP.
• Grains reaching or leaving equilibrium: opacities now accumulated and updated
  in both directions.
• Smallest nLTE grain never selected; wrong J0 index in qRE re-emission.
• calc_Jth: uninitialised constant and grain zone for nRE grains.
• Re-emission CDFs now include the first wavelength bin.
• Misselt files with several atom types misread.

Performance

• Temp_nRE: transition edges located once per grain, enthalpy once per run.
• calc_Jth: once per wavelength instead of per inclination, tabulated Planck
  function, OpenMP loops.
• nLTE/qRE re-emission no longer sums all threads at each event.
• Bit-for-bit identical on 1 thread.
• OzSTAR, 16 threads, W5 (PAH): temperature_nRE 0.72 -> 0.29 s, source_fct
  2.13 -> 1.30 s, total 8.46 -> 7.05 s, despite 2.5x more stochastic grain-cell
  pairs. W5 output 1.3x slower (larger Temperature_nRE.fits.gz). W1-W4 unchanged.

Robustness, options, tests

• Temperature files from another grid, grain distribution or temperature grid
  (-nT, -Tmax_PAH) are refused instead of overflowing the arrays.
• -Mathis_field accepted in any order with -benchmark_SHG.
• Documented -nT, -benchmark_SHG, -Mathis_field; -Tmax_PAH default is 3000 K.
• Camps 2015 driver (test_suite/benchmarks/Camps2015) and slow test; data copied
  to utils/benchmarks/Camps2015 (next utils release).
• ref3.0_multi replaced by ref4.1_multi (same model, 4.1 format); ref3.0 kept for
  backward compatibility. Test data republished.

Testing:

• Camps 2015 benchmark: see above; 1 and 6 threads identical.
• Full local suite including slow tests: 136 passed, 11 skipped.
• New fast CI tests: mismatched temperature files, -Mathis_field option order.

Did you run the botscheck that the code and comments follow the code of conduct? yes

Did you update relevant documentation in the docs directory? yes

[github-actions]

ERROR: # Code Review by Gemini
An unexpected error occurred in Client call: 503 UNAVAILABLE. {'error': {'code': 503, 'message': 'This model is currently experiencing high demand. Spikes in demand are usually temporary. Please try again later.', 'status': 'UNAVAILABLE'}}

[github-actions]

Code Review by Gemini

Here is the code review for the "pah nre" pull request.

---

Physical & Numerical Purpose of Changed Routines

• PAH.f90: Computes heat capacities and integrates enthalpies $U(T) = \int_0^T C(T') dT'$ for PAH, AstroSil, Misselt, and DustEM grain populations.
• thermal_emission.f90: Implements radiative equilibrium and non-equilibrium (stochastic heating) thermal solver. Constructs transition matrices $A_{fi}$ for temperature probability distributions $P(T)$, calculates emission rates, and updates energy absorption/re-emission bounds across grain sizes.
• dust_ray_tracing.f90: Precomputes thermal emissivity $J_{\text{th}}$ across spatial cells and wavelength bins for ray tracing.
• input.f90 / read_opacity.f90: Loads temperature FITS restart files and Misselt optical property files safely.

---

Detailed Review Findings

1. Radiative Transfer Physics & Energy Conservation

• Top Bin Transition Rate & Energy Conservation (Temp_nRE in src/thermal_emission.f90):
  Integrating heating transition rates over discrete frequency intervals $N_{\text{lo}} - N_{\text{hi}}$ (using Draine & Li 2001 formulation) and accumulating all transitions exceeding the highest temperature boundary into the top bin ($A(n_T, T_2) = N_{\text{lo}}$) ensures exact energy conservation regardless of wavelength grid spacing.
• Removal of Empirical Factors & $P(T)$ Normalization (Temp_nRE in src/thermal_emission.f90):
  Removing the legacy reset X(1, id) = X(2, id) and normalizing $X$ directly ensures that low-temperature probability states are mathematically consistent.
• Correct Wavelength Index in Background Field (im_reemission_qRE in src/thermal_emission.f90):
  Replacing J0(icell, lambda) with J0(icell, ilambda) inside the wavelength integration loop fixes a bug where the background field intensity was previously sampled at the incoming photon's wavelength instead of the integration wavelength ilambda.

2. Parallel Computing & Thread Safety

• Thread-Local Energy Accumulation (propagate_packet in src/dust_transfer.f90 & src/thermal_emission.f90):
  Replacing the global scalar E_abs_nRE reduction with thread-indexed array E_abs_nRE_id(id) during photon propagation avoids race conditions and OpenMP atomic/reduction overhead.
• calc_Jth Hoisting (src/dust_ray_tracing.f90 & src/dust_transfer.f90):
  Hoisting calc_Jth(lambda) outside the ray-tracing inclination loop (ibin) avoids redundant thermal emissivity calculations per inclination angle while maintaining thread safety.
• Cache-Friendly Matrix Indexing (Temp_nRE in src/thermal_emission.f90):
  Transposing the temporary transition matrix $B(T_2, T, id)$ aligns the inner recursion loops with Fortran column-major memory layout, improving SIMD/vectorization and CPU cache hit rates.

3. Numerical Precision & Array Bounds Safety

• Enthalpy Low-Temperature Model (enthalpy in src/PAH.f90):
  Using enthalpy(1) = 0.25_dp * C_sub(1) * T(1) correctly assumes a Debye $C(T) \propto T^3$ relation near $T \to 0$, giving exact initial enthalpy $U(T_1) = \int_0^{T_1} C(T') dT' = \frac{1}{4} C(T_1) T_1$.
• Heat Capacity Unit Conversions (Misselt_specific_heat & DustEM_specific_heat in src/PAH.f90):
  Converting specific heat from $\text{erg},\text{g}^{-1}\text{K}^{-1}$ to SI $\text{J/K}$ (1.e-7) corrects the previously overstated values.
• Smallest Grain Selection in nLTE (im_reemission_NLTE in src/thermal_emission.f90):
  Starting the search at kmin = grain_RE_nLTE_start - 1 ensures that the smallest grain size in non-LTE equilibrium (grain_RE_nLTE_start) can be chosen when rand1 <= CDF(grain_RE_nLTE_start).
• FITS Header & Grid Validation (check_temperature_hdu in src/input.f90):
  Adding HDU grid and size checks prevents array overflow when loading precalculated temperature files (Temperature_nRE.fits.gz) generated with mismatched temperature grids (-nT, -Tmax_PAH) or size distributions.

---

Overall Assessment

• Quality: High. The PR systematically eliminates empirical tuning parameters, restores rigorous energy conservation for non-equilibrium grains, resolves multi-threading race conditions, and adds comprehensive benchmark verification (Camps et al. 2015).
• Risks: Low. High-performance paths remain bit-for-bit identical in single-threaded runs while being thread-safe under OpenMP.
• Recommendation: Safe to merge.

[github-actions]

Code Review by Gemini

Code Review Summary

Overall, this is a well-engineered pull request that addresses long-standing energy conservation issues in MCFOST's treatment of non-equilibrium grains (PAHs and nRE dust). The removal of arbitrary empirical normalization factors, the introduction of proper enthalpy integration, frequency-bin edge matching according to Draine & Li (2001), and thread-private tracking of absorbed energy (E_abs_nRE_id) represent significant physics and numerical improvements.

One HIGH severity bug was identified regarding the reset lifecycle of E_abs_nRE_id across multiple thermal Monte Carlo runs.

---

Key Findings

1. [HIGH] Un-reset E_abs_nRE_id across multiple run_thermal_mc() calls

• File: src/dust_transfer.f90 (lines 901–945) & src/thermal_emission.f90

• Description:
  E_abs_nRE_id tracks the energy absorbed by non-equilibrium grains per OpenMP thread during packet propagation (propagate_packet):

  E_abs_nRE_id(id) = E_abs_nRE_id(id) + Stokes(1) * (1.0_dp - proba_abs_RE(icell,lambda))

  E_abs_nRE_id is reset to zero inside emission_nRE():

  E_abs_nRE = sum(E_abs_nRE_id)
  ...
  E_abs_nRE = 0.0 ; E_abs_nRE_id = 0.0_dp

  However, when Temp_nRE determines that the nRE grains have converged (flag_em_nRE = .true.), or when -benchmark_SHG is set, emission_nRE() is not called. As a result, E_abs_nRE_id retains non-zero leftover values from the final Monte Carlo iteration.

  If run_thermal_mc() is called again later in the run (e.g., during hydrostatic equilibrium iterations or multi-step physics loops), E_abs_nRE_id is not zeroed before starting the new thermal photon loop. propagate_packet will accumulate onto these stale values, causing E_abs_nRE to overestimate the absorbed energy in subsequent nRE iterations.

• Impact: Incorrect energy balance and biased PAH/nRE re-emission spectra whenever run_thermal_mc() is executed more than once per run (e.g., in hydrostatic models).

• Suggested Fix:
  Zero out E_abs_nRE and E_abs_nRE_id at the start of run_thermal_mc() (or within reset_temperature()):

  ! In src/dust_transfer.f90 inside run_thermal_mc():
  call reset_radiation_field()
  call reset_temperature()
  call reset_sed()
  if (lnRE) then
     E_abs_nRE = 0.0_dp
     E_abs_nRE_id = 0.0_dp
  endif

---

Other Observations & Highlights

1. Energy Conservation & Enthalpy Calculation (src/PAH.f90):
   The new enthalpy function uses standard Debye low-$T$ behavior ($\int_0^{T_1} C(T) dT = \frac{1}{4} C(T_1) T_1$) and trapezoidal sub-stepping in log-space. Unit conversion (1.e-7) from CGS calorimetry files to SI ($J/K$) is correct.

2. Frequency Transition Edge Pre-calculation (src/thermal_emission.f90):
   The photon_edge helper and k_edge mapping cleanly determine transition bounds $U_{\text{lim}}(T)/h - \nu_{\text{bin}}(T_2)$ without runtime search overhead. Array indices and bounds for N_above (0:n_lambda) and kJnu_o_nu (n_lambda) are safe.

3. Mismatched FITS Temperature Guards (src/input.f90):
   Checking HDU 3 (tab_Temp) and HDU dimensions against n_T and Tmax_PAH before reading prevents buffer overflows when reusing output files with altered parameter flags.

---

Merge Recommendation

Conditional Merge. The PR is physically solid and significantly improves benchmark accuracy. However, the E_abs_nRE_id initialization issue should be fixed before merging to prevent subtle energy non-conservation bugs in multi-iteration / hydrostatic models.

[github-actions]

Code Review by Gemini

Here is a review of the pull request.

---

Physical & Numerical Summary

The changes significantly improve the physical consistency, energy conservation, and numerical stability of non-equilibrium grain (PAHs, nRE) modeling in MCFOST:

1. Energy Conservation: Heating rates and transition matrices ($A_{fi}$) in Temp_nRE are now integrated across exact frequency bin widths (following Draine & Li 2001) rather than evaluated at point samples, ensuring photon/energy conservation across arbitrary wavelength grids without empirical rescaling.
2. Enthalpy Calculation: Grain enthalpy $U(T)$ is properly integrated from the heat capacity $C(T)$ with low-temperature Debye scaling ($C \propto T^3$) rather than approximated by $C(T) T$.
3. Unit & Opacity Consistency: Heat capacities for DustEM and Misselt optical files are corrected to SI units ($\text{J}\cdot\text{K}^{-1}$), and transition probabilities in update_proba_abs_nRE now correctly handle grain transitions in both directions (reaching or leaving equilibrium).

---

Key Finding

1. Accumulator E_abs_nRE_id is not reset prior to photon propagation across run_thermal_mc calls

Severity: HIGH
Location: src/dust_transfer.f90 (lines 783–787) / src/thermal_emission.f90 (emission_nRE)

Problem:
In dust_transfer.f90, the thread-local array E_abs_nRE_id accumulates energy absorbed by non-equilibrium grains during packet propagation:

E_abs_nRE_id(id) = E_abs_nRE_id(id) + Stokes(1) * (1.0_dp - proba_abs_RE(icell,lambda))

E_abs_nRE_id is zeroed inside emission_nRE() when non-equilibrium re-emission iterations occur. However, on the final iteration of run_thermal_mc() (when flag_em_nRE is .true.), emission_nRE() is not called. As a result, E_abs_nRE_id retains the values accumulated during the final iteration's photon loop.

If run_thermal_mc() is subsequently invoked a second time (e.g. during hydrostatic iterations as in test_hydrostatic_nre, or when MCFOST is called via python bindings/C API), E_abs_nRE_id is not zeroed upon entry. mc_photon_loop will then accumulate energy on top of the stale values from the previous thermal calculation, causing overestimation of absorbed PAH energy.

Fix:
Zero E_abs_nRE_id (and E_abs_nRE) at the start of mc_photon_loop when letape_th is true, or inside reset_temperature() / reset_thermal_state():

! In src/dust_transfer.f90 inside mc_photon_loop:
if (letape_th) then
   p_nnfot2 => nnfot2
   if (allocated(E_abs_nRE_id)) E_abs_nRE_id = 0.0_dp
   E_abs_nRE = 0.0_dp
else
...

---

Minor Observations & Code Quality

1. OpenMP Loop Parallelization in calc_Jth (src/dust_ray_tracing.f90):
   The parallelization of calc_Jth over icell (!$omp parallel do schedule(dynamic,100)) pre-computes $B_\lambda(T)$ once per wavelength into B_wl(T) before the grid loop. This significantly speeds up thermal emissivity evaluation for nRE grains without thread contention.

2. Array Bounds Checking in FITS Temperature Reader (src/input.f90):
   The new check_temperature_hdu validation cleanly guards against out-of-bounds reads and memory corruption when loading pre-computed temperature files with mismatched temperature or spatial grids.

---

Summary

• Overall Assessment: Outstanding PR that fixes long-standing energy conservation and unit issues in PAH/nRE stochastic heating.
• Most Important Issue: Ensure E_abs_nRE_id is zeroed before mc_photon_loop begins during thermal structure calculations so that stale absorbed energy does not spill across multiple run_thermal_mc calls.
• Merge Recommendation: Safe to merge once the E_abs_nRE_id reset logic is addressed.

[github-actions]

Code Review by Gemini

Code Review Summary

Overview

This PR is a major physics fix and numerical improvement for the stochastic heating of non-equilibrium grains (PAHs and nRE grains) in MCFOST. It replaces previous empirical correction factors with exact energy-conserving enthalpy integration and heating transition matrices (following Draine & Li 2001 / Guhathakurta & Draine 1989), corrects specific heat unit conversions, fixes OpenMP thread safety during photon re-emission, and introduces bounds validation when reading FITS temperature files.

---

Key Findings & Verification

1. Energy Conservation & Non-Equilibrium Physics

• Enthalpy & Integration (src/PAH.f90, src/thermal_emission.f90):
  • enthalpy(T, igrain) integrates $C(T') dT'$ starting from the $T^3$ Debye low-temperature limit up to $T$.
  • Specific heat capacities in Misselt_specific_heat and DustEM_specific_heat were correctly converted from erg/K to J/K ($10^{-7}$ factor).
  • Frequency transitions in Temp_nRE are computed using exact frequency edge integrations (photon_edge and N_above), ensuring photon energy conservation on any wavelength grid.
• Grain State Transitions (src/thermal_emission.f90):
  • update_proba_abs_nRE now correctly tracks delta_kappa_abs_qRE when grains transition in both directions (from stochastic heating to pseudo-equilibrium qRE, and back to nRE).

2. Parallel Computing & Thread Safety

• OpenMP Reductions & Contention (src/dust_transfer.f90, src/thermal_emission.f90):
  • E_abs_nRE reduction replaced with a per-thread array E_abs_nRE_id(id), preventing OpenMP thread races during photon propagation.
  • im_reemission_NLTE and im_reemission_qRE now estimate mean intensity using thread-local accumulation (xJ_abs(icell, ilambda, id) * nb_proc), eliminating cross-thread cache thrashing during photon propagation loops.
  • calc_Jth and repartition_energie OpenMP parallelizations are indexed cleanly by icell without data races.

3. Input Validation & Robustness

• FITS Header Checks (src/input.f90):
  • check_temperature_hdu and HDU 3 verification in lect_Temperature check temperature grid dimensions (n_T) and bounds (tab_Temp) before loading pre-computed temperature files, preventing array overflows when command-line parameters (-nT, -Tmax_PAH) change.
• Misselt Opacity Header Parsing (src/read_opacity.f90):
  • misselt_load now correctly skips ibuffer + 1 header lines (where ibuffer is the number of atom types), fixing reading errors for files with non-standard atom counts.

---

Recommendation

Status: Safe to merge.

The code is well-structured, mathematically sound, thread-safe, and passes the 51-case Camps et al. (2015) benchmark.