Adding hydrostatic equilibrium

[cpinte]

Type of PR:
This adds new physics and implement an iterative loop to compute the hydrostatic equilibrium

Description:
Iterate between gas hydrostatic equilibrium, dust setlling and thermal calculations

Testing:
No tests yet

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? no

[github-actions]

Code Review by Gemini

Here is the code review for the pull request "Adding hydrostatic equilibrium".

---

Physical & Numerical Understanding

This PR implements an iterative hydrostatic equilibrium loop for the gas disk structure coupled with vertical dust settling redistribution.

1. dust_transfer_sub() in dust_transfer.f90 drives the physics iteration: after computing temperatures via Monte Carlo, it invokes equilibre_hydrostatique() and checks for convergence of $T_{\text{dust}}$ across iterations.
2. equilibre_hydrostatique() in disk_physics.f90 solves vertical hydrostatic balance for the gas in cylindrical columns, smooths vertical temperatures to damp Monte Carlo noise, enforces a floor on density decline rates, under-relaxes updates, and normalizes columns to preserve gas surface density ($\Sigma_{\text{gas}}$).
3. redistribute_dust_hydrostatic() in density.f90 redistributes dust vertically following the newly computed non-Gaussian numerical gas shape $n_{\text{gas}}(z)^{(H_{\text{gas}}/H_{\text{dust}})^2}$, supporting both power-law and Dubrulle settling prescriptions while conserving radial column dust mass.

---

Key Findings

1. CRITICAL: Stale Sigma_gas_target persisted across model runs

• File: src/disk_physics.f90 (lines 249–262)
• Problem: Sigma_gas_target is declared as a SAVE, ALLOCATABLE local variable:

  real(kind=dp), save, allocatable, dimension(:,:) :: Sigma_gas_target
  ...
  if (.not.allocated(Sigma_gas_target)) then
     allocate(Sigma_gas_target(n_rad,n_az), stat=alloc_status)
     ...
  endif

  It is allocated on the first call to equilibre_hydrostatique(), but it is never deallocated or recomputed for subsequent model runs.
• Why it matters: When MCFOST is run in library mode (libmcfost), C/Python wrappers, or parameter study sweeps where dust_transfer_sub() is called multiple times in the same executable process:
  1. If a subsequent run uses different grid dimensions (n_rad or n_az), accessing Sigma_gas_target(i,k) causes an out-of-bounds memory access / segfault.
  2. If a subsequent run uses identical dimensions but a different disk mass or profile, the hydrostatic solver will incorrectly force the new model's gas surface density to match the first model's surface density.
• Suggested Fix: Reset or recompute Sigma_gas_target whenever grid dimensions change or at the start of a physics loop (e.g. when iter == 1), or deallocate Sigma_gas_target when dust_transfer_sub() completes. For example, in equilibre_hydrostatique():

  if (allocated(Sigma_gas_target)) then
     if (size(Sigma_gas_target,1) /= n_rad .or. size(Sigma_gas_target,2) /= n_az) then
        deallocate(Sigma_gas_target)
     endif
  endif

  if (.not.allocated(Sigma_gas_target)) then
     allocate(Sigma_gas_target(n_rad,n_az), stat=alloc_status)
     if (alloc_status > 0) call error("Allocation error Sigma_gas_target")
     do k=1, n_az
        do i=1, n_rad
           total_sum = 0.0_dp
           do j=1, nz
              total_sum = total_sum + gas_density(cell_map(i,j,k)) * cell_height(i,j)
           enddo
           Sigma_gas_target(i,k) = total_sum
        enddo
     enddo
  endif

And expose a public cleanup routine (or deallocate Sigma_gas_target when dust_transfer_sub finishes) so new model parameter sets reset the target surface density.

---

2. MEDIUM: Incorrect step size dz_local used for vertical grid derivative and Euler integration

• File: src/disk_physics.f90 (lines 283–290)
• Problem: In equilibre_hydrostatique(), dz_local is assigned to cell_height(i,j) ($\Delta z_j$):

  do j = 2, nz
     dz_local = cell_height(i,j)
     icell = cell_map(i,j,k)
     dTdz = (Tsmooth(j)-Tsmooth(j-1)) / dz_local
     fac1 = cst * z_grid(icell)/ (r_grid(icell)**3)
     fac2 = -1.0_dp * (dTdz + fac1) / Tsmooth(j)
     fac2 = min(fac2, -min_decline_fraction * fac1 / Tsmooth(j))
     ln_rho(j) = ln_rho(j-1) + fac2 * dz_local
     rho(j) = exp(ln_rho(j))
  enddo

  On non-uniform vertical grids (where cell heights $\Delta z_j$ stretch with altitude), the distance between adjacent cell centers $z_{\text{grid}}(j) - z_{\text{grid}}(j-1)$ is not equal to the cell height cell_height(i,j).
• Why it matters: Using cell height $\Delta z_j$ instead of center-to-center distance $\Delta z_{\text{centers}}$ introduces artificial discretization errors in both $dT/dz$ and the cumulative Euler integration step for $\ln \rho(z)$, distorting the hydrostatic scale height on non-uniform grids.
• Suggested Fix: Compute the distance between cell centers for the integration step:

  real(kind=dp) :: dz_centers

  do j = 2, nz
     icell = cell_map(i,j,k)
     dz_centers = z_grid(icell) - z_grid(cell_map(i,j-1,k))
     dTdz = (Tsmooth(j)-Tsmooth(j-1)) / dz_centers
     fac1 = cst * z_grid(icell) / (r_grid(icell)**3)
     fac2 = -1.0_dp * (dTdz + fac1) / Tsmooth(j)
     fac2 = min(fac2, -min_decline_fraction * fac1 / Tsmooth(j))
     ln_rho(j) = ln_rho(j-1) + fac2 * dz_centers
     rho(j) = exp(ln_rho(j))
  enddo

---

3. LOW: Potential negative base in gas_shape exponentiation

• File: src/density.f90 (lines 1022–1025)
• Problem: In get_gas_shape():

  gas_shape(1) = 1.0_dp
  do j=2, nz
     gas_shape(j) = min(gas_density(cell_map(i,j,k)) / gas_ref, gas_shape(j-1))
  enddo

  If gas_density is negative for any cell (e.g. from numerical precision noise or uninitialized boundary cells), gas_shape(j) will become negative.
• Why it matters: In redistribute_dust_column(), computing new_density(j) = gas_shape(j)**correct_strat_l with a negative base and real exponent raises a floating-point invalid operation (NaN or runtime crash).
• Suggested Fix: Clamp gas_shape(j) to be non-negative:

  gas_shape(1) = 1.0_dp
  do j=2, nz
     gas_shape(j) = max(0.0_dp, min(gas_density(cell_map(i,j,k)) / gas_ref, gas_shape(j-1)))
  enddo

---

Summary

• Overall Assessment: The physical model and implementation structure are sound. The coupling between gas hydrostatic balance, temperature smoothing, and dust settling redistribution is logically clean.
• Most Important Issue: Finding Update density.f90 #1 (Sigma_gas_target persistence) must be fixed before merging to prevent severe bugs when running multiple models or using libmcfost/Python bindings.
• Safe to merge? No, fix Finding Update density.f90 #1 (and ideally correct specie -> species (even singular) #2 and missing documentation: special output data formats #3) first.

[github-actions]

Code Review by Gemini

Overview

This PR implements an iterative hydrostatic equilibrium loop for MCFOST, coupling the gas vertical density calculation to the thermal radiative transfer structure (Tdust) and re-equilibrating the dust vertical distribution consistently.

Overall, the physical model and numerical implementation are sound and well-thought-out. The vertical dust redistribution correctly conserves column mass, dimensional unit conversions are handled accurately, and the restart mechanism via iteration snapshots (Temperature_iter*.fits.gz and gas_density_iter*.fits.gz) provides good usability.

Below are a few issues identified during the review regarding memory persistence and potential floating-point division edge cases.

---

Findings

1. Sigma_gas_target state persistence across multi-run or library invocations

• Severity: MEDIUM
• Location: src/disk_physics.f90, lines ~219–235 (equilibre_hydrostatique)
• Explanation:
  In equilibre_hydrostatique, Sigma_gas_target is declared as:

  real(kind=dp), save, allocatable, dimension(:,:) :: Sigma_gas_target

  and allocated conditionally:

  if (.not.allocated(Sigma_gas_target)) then
     allocate(Sigma_gas_target(n_rad,n_az), stat=alloc_status)
     ...
  endif

  Because of the SAVE attribute, Sigma_gas_target persists across execution calls. If MCFOST is invoked as a library (lmcfost_lib), or if multiple hydrostatic simulations are executed sequentially in the same process with different grid parameters or initial gas density profiles, Sigma_gas_target will retain the target surface density from the first run, corrupting subsequent models or causing array bounds mismatches if n_rad or n_az change.
• Suggested Fix:
  Deallocate Sigma_gas_target when grid/density structures are initialized/reset, or allow equilibre_hydrostatique to re-allocate if grid dimensions change. For instance:

  if (allocated(Sigma_gas_target)) then
     if (size(Sigma_gas_target,1) /= n_rad .or. size(Sigma_gas_target,2) /= n_az) then
        deallocate(Sigma_gas_target)
     endif
  endif
  
  if (.not.allocated(Sigma_gas_target)) then
     allocate(Sigma_gas_target(n_rad,n_az), stat=alloc_status)
     ...
  endif

---

2. Potential division by zero on Tsmooth(j) in equilibre_hydrostatique

• Severity: MEDIUM
• Location: src/disk_physics.f90, lines ~264–266 (equilibre_hydrostatique)
• Explanation:
  The hydrostatic gradient calculation evaluates:

  fac2 = -1.0_dp * (dTdz + fac1) / Tsmooth(j)
  fac2 = min(fac2, -min_decline_fraction * fac1 / Tsmooth(j))

  If Tsmooth(j) evaluates to 0.0_dp (e.g., if unilluminated atmospheric cells or uninitialized temperature values occur), division by Tsmooth(j) will produce NaN or floating-point exceptions.
• Suggested Fix:
  Clamp Tsmooth(j) to a positive minimum floor (e.g., 1.0_dp or tiny_dp):

  real(kind=dp) :: T_cell
  ...
  T_cell = max(Tsmooth(j), 1.0_dp)
  fac2 = -1.0_dp * (dTdz + fac1) / T_cell
  fac2 = min(fac2, -min_decline_fraction * fac1 / T_cell)

---

3. Monte Carlo noise sensitivity in hydrostatic convergence test

• Severity: LOW
• Location: src/dust_transfer.f90, lines ~380–400 (dust_transfer_sub)
• Explanation:
  The convergence check tests relative temperature difference in every cell:

  rel_diff = abs(Tdust(icell)-Tdust_hydro_prev(icell)) / Tdust_hydro_prev(icell)
  if (rel_diff > hydrostatic_precision) lconverged = .false.

  In high-altitude atmospheric cells or optically thick midplane regions receiving very few photon packets, Monte Carlo Poisson noise can cause single-cell temperature fluctuations $> 5%$ (hydrostatic_precision = 0.05), preventing convergence even when the overall disk density and temperature structure have stabilized.
• Suggested Fix:
  Consider restricting the max-change convergence test to cells above a minimum density threshold or minimum packet count, or requiring both mean_dT_o_T and max_dT_o_T (or a 95th percentile) to satisfy tolerance criteria.

---

Summary Assessment

• Overall Quality: High. The PR implements hydrostatic equilibrium and dust vertical re-equilibration in a physically consistent and numerically clean manner.
• Key Issues:
  1. Clear saved Sigma_gas_target allocations when model parameters re-initialize to avoid state leakage in library/multi-run modes.
  2. Add zero-division safeguards around Tsmooth(j) in equilibre_hydrostatique.
• Merge Recommendation: Safe to merge after addressing Findings 1 and 2.

[github-actions]

Code Review by Gemini

Here is a code review of the pull request.

---

Summary of Assessment

This PR adds iterative gas hydrostatic equilibrium coupled with vertical dust settling and Monte Carlo thermal radiative transfer. Overall, the physical formulation and numerical convergence controls (under-relaxation, photon-noise floor adjustments, monotonicity enforcement) are well constructed.

However, there are critical issues regarding silent output corruption for thermal SEDs, state persistence bugs in repeated/library executions, and a potential floating-point exception hazard.

---

Key Findings & Recommendations

---

1. CRITICAL: Thermal SED Output is Wiped in Hydrostatic Mode

• Location: src/dust_transfer.f90 (lines 142–180, 281–295) & src/output.f90

• Problem:
  During hydrostatic iterations, run_thermal_mc(lwrite_output=.false.) is invoked to compute temperatures. Because lwrite_output is .false., ecriture_sed(1) is skipped inside run_thermal_mc. However, at the end of run_thermal_mc, sed and associated flux arrays are unconditionally reset to zero (sed = 0.0; sed_q = 0.0 ...).

  Once the hydrostatic loop finishes in dust_transfer_sub, it calls ecriture_sed(1). But by this point, sed has already been cleared to zero on the final call to run_thermal_mc.

• Impact:
  In hydrostatic mode (lhydrostatic = .true.), .sed_th.fits.gz is written containing all zeroes, silently corrupting the thermal SED output.

• Suggested Fix:
  In run_thermal_mc, do not clear the sed arrays at the end if lwrite_output is .false., or preserve them until ecriture_sed(1) has been called:

  if ((.not. lmcfost_lib) .and. lwrite) then
     call ecriture_temperature(1)
     call ecriture_sed(1)
  endif

  ! Only reset SED arrays if they were written; otherwise let caller read them
  if (lwrite) then
     sed = 0.0; sed_q = 0.0; sed_u = 0.0; sed_v = 0.0
     n_phot_sed = 0.0; n_phot_sed2 = 0.0; n_phot_envoyes = 0.0
     sed_star = 0.0; sed_star_scat = 0.0; sed_disk = 0.0; sed_disk_scat = 0.0
  endif

---

2. HIGH: State Persistence & Out-of-Bounds Hazard in equilibre_hydrostatique

• Location: src/disk_physics.f90 (lines 255–270)

• Problem:
  Sigma_gas_target is declared with the save attribute:

  real(kind=dp), save, allocatable, dimension(:,:) :: Sigma_gas_target

  Because Sigma_gas_target is saved across subroutine calls, it remains allocated across multiple model runs when MCFOST is used via library mode (lmcfost_lib), Python wrappers (pymcfost), or sequential test invocations.

  If a subsequent model run uses different grid dimensions (n_rad or n_az), Sigma_gas_target will not be reallocated and will cause an out-of-bounds memory access or crash. Even if grid dimensions match, the second model will erroneously normalize its gas density to the target surface density of the first model.

• Impact:
  Memory corruption, crashes, or incorrect density normalization when running multiple models in a single process.

• Suggested Fix:
  Deallocate Sigma_gas_target whenever model grid parameters change or when initializing a new run, or re-verify its allocation against current grid dimensions:

  if (allocated(Sigma_gas_target)) then
     if (size(Sigma_gas_target,1) /= n_rad .or. size(Sigma_gas_target,2) /= n_az) then
        deallocate(Sigma_gas_target)
     endif
  endif

  if (.not. allocated(Sigma_gas_target)) then
     allocate(Sigma_gas_target(n_rad, n_az), stat=alloc_status)
     ...

(Also add if (allocated(Sigma_gas_target)) deallocate(Sigma_gas_target) to reset routines in mem.f90 or init_mcfost.f90.)

---

3. MEDIUM: Potential Floating-Point Exception in redistribute_dust_column

• Location: src/density.f90 (lines 1045–1050)
• Problem:
  redistribute_dust_column computes:

  new_density(j) = gas_shape(j)**correct_strat_l

  gas_shape(j) can be 0.0_dp (e.g., above the gas disk). In Fortran, evaluating 0.0_dp ** correct_strat_l when correct_strat_l == 0.0_dp (e.g., if no settling is active) or under strict IEEE compiler checks can trigger an invalid operation domain error or floating-point exception (FPE).
• Impact:
  Potential runtime crash on compilers compiled with strict FPE traps (e.g., gfortran -ffpe-trap=invalid or ifort -fpe0).
• Suggested Fix:
  Guard the exponentiation against zero bases:

  do j=1, nz
     if (gas_shape(j) > 0.0_dp) then
        new_density(j) = gas_shape(j)**correct_strat_l
     else
        new_density(j) = 0.0_dp
     endif
     N_new = N_new + new_density(j) * cell_height(i,j)
  enddo

---

4. LOW / PERFORMANCE: Array Temporary Creation in Hot Hydrostatic Loop

• Location: src/disk_physics.f90 (lines 278–282)
• Problem:
  Inside equilibre_hydrostatique:

  Tsmooth(j) = sum(Tdust(cell_map(i,max(1,j-2):min(nz,j+2),k))) / real(min(nz,j+2)-max(1,j-2)+1,dp)

  cell_map(i, :, k) is non-contiguous in Fortran column-major memory layout. Array section slicing with non-unit strides inside nested loops (n_az * n_rad * nz) forces the compiler to generate heap/stack array temporaries on every cell iteration.
• Impact:
  Unnecessary overhead and allocation thrashing during hydrostatic updates.
• Suggested Fix:
  Use a scalar loop to sum Tdust values:

  do j=1, nz
     tsum = 0.0_dp
     n_pts = 0
     do j_sub = max(1, j-2), min(nz, j+2)
        tsum = tsum + Tdust(cell_map(i, j_sub, k))
        n_pts = n_pts + 1
     enddo
     Tsmooth(j) = tsum / real(n_pts, dp)
  enddo

---

Merge Recommendation

Status: Changes Requested

The hydrostatic equilibrium implementation is physically well designed, but Finding 1 (corrupted/zeroed thermal SED outputs) and Finding 2 (state persistence bug in Sigma_gas_target) should be addressed before merging into the main branch.