mod_eos_lte_saha.t

Go to the documentation of this file.

!=============================================================================
!> Analytic hydrogen Saha EoS with neutral helium (eos_method == 'analytic')
!>
!> Ground-state Saha for hydrogen, solved per call without tables:
!>   X(T)        = saha_pf * T^{3/2} * exp(-chi/kB/T)
!>   y(=ne/nH)   = 2X / (X + sqrt(X^2 + 4 nH X))
!> Public entry points take and return code units; the CGS conversion is internal.
!> T(nH, eint/nH) is found by bisection and p(nH, T) closes the state.
!> build_gamma1_analytic_table tabulates Gamma_1(nH, p/nH) once for the 'exact' gamma1 path,
!> on the axes that gamma1_from_nH_p reads for every method
!=============================================================================
module mod_eos_lte_saha
    use mod_global_parameters
    use mod_eos_container, only: eos, eos_analytic
    use mod_eos_interp,    only: precompute_step_inv
    use mod_comm_lib,      only: mpistop
    use mod_eos_lte_tables, only: precompute_fi_bypass_constants
    implicit none
    private
 
    public :: saha_eint_from_nh_t, saha_state_from_nh_p
    public :: saha_t_from_nh_eint
    public :: load_analytic_lte, finalise_analytic_lte
 
    ! Prefactor (2 pi m_e k_B/h^2)^{3/2} in CGS (cm^-3)
    ! SI value is 2.4146830395719654e21 m^-3; divide by 1e6 for CGS
    double precision, parameter :: saha_pf = 2.4146830395719654d15
    double precision, parameter :: saha_chi_kB = 157763.42386247337d0
    double precision, parameter :: saha_kB_cgs = 1.380649d-16
    double precision, parameter :: saha_chi_H_cgs = 2.178710282685096d-11
 
    !> One-shot guard: a bisection that exhausts max_iter without meeting the
    !> tolerance warns once (rank 0) rather than silently returning an
    !> under-converged T every offending cell.
    logical :: saha_nonconv_warned = .false.
 
contains
 
    !> Analytic-method load: no binary tables (on-the-fly Saha solve); just
    !> announce the method.
    subroutine load_analytic_lte()
        if (mype == 0) then
            write(*,*) 'EoS method: analytic (hydrogen Saha, helium neutral)'
            write(*,*) '  gamma1_method: ', trim(eos%gamma1_method)
        end if
    end subroutine load_analytic_lte
 
    !> Analytic-method finalise: no table unit conversion needed;
    !> build the Gamma1 table only for 'exact' mode, then the shared FI-bypass
    !> constants.
    subroutine finalise_analytic_lte()
        if (eos%gamma1_method == 'exact') then
            call build_gamma1_analytic_table()
        end if
        call precompute_fi_bypass_constants(saha_chi_h_cgs)
    end subroutine finalise_analytic_lte
 
    !> Ionization fraction from Saha in CGS (no unit conversions).
    !> Use this inside bisection loops where nH_cgs is constant.
    double precision function saha_y_cgs(nH_cgs, T_cgs) result(y)
        double precision, intent(in) :: nh_cgs, t_cgs
        double precision :: x_saha, disc
 
        x_saha = saha_pf * t_cgs * dsqrt(t_cgs) * dexp(-saha_chi_kb / t_cgs)
        disc = dsqrt(x_saha * x_saha + 4.0d0 * nh_cgs * x_saha)
        ! Use numerically stable form to avoid catastrophic cancellation:
        ! When X >> nH (high ionization): standard form loses precision in -X + sqrt(...)
        ! Stable form: y = 2*X / (X + sqrt(X^2 + 4*nH*X)) -> 1 as X -> inf
        y = 2.0d0 * x_saha / (x_saha + disc)
 
    end function saha_y_cgs
 
    !> Ionisation fraction y = ne/nH from Saha at (nH, T) in code units
    double precision function saha_y_from_nh_t(nH_code, T_code) result(y)
        double precision, intent(in) :: nh_code, t_code
 
        y = saha_y_cgs(nh_code * unit_numberdensity, t_code * unit_temperature)
 
    end function saha_y_from_nh_t
 
    !> eint/nH in CGS with helium neutral, 1.5 (1+He+y) kB T, plus y chi_H when ionE
    double precision function saha_eint_cgs(y, T_cgs) result(e)
        double precision, intent(in) :: y, t_cgs
        e = 1.5d0 * (1.0d0 + eos%He_abundance + y) * saha_kb_cgs * t_cgs
        if (eos%ionE) e = e + y * saha_chi_h_cgs
    end function saha_eint_cgs
 
    !> eint/nH from Saha at (nH, T), in code units
    double precision function saha_eint_from_nh_t(nH_code, T_code) result(eint_nH_code)
        double precision, intent(in) :: nh_code, t_code
 
        if (eos%method_id /= eos_analytic) call mpistop("saha_eint_from_nH_T requires eos_method=analytic")
        eint_nh_code = saha_eint_cgs(saha_y_from_nh_t(nh_code, t_code), t_code * unit_temperature) &
            * unit_numberdensity / unit_pressure
 
    end function saha_eint_from_nh_t
 
    !> Given (nH, p) in code units, T, y and optionally eint/nH: bisection on
    !> p = nH*(1+He+y(T))*T, with kB absorbed in the code units
    subroutine saha_state_from_nh_p(nH_code, p_code, T_out, y_out, eint_nH_out)
        double precision, intent(in)  :: nh_code, p_code
        double precision, intent(out) :: t_out, y_out
        double precision, intent(out), optional :: eint_nh_out
 
        double precision :: nh_cgs, ut
        double precision :: t_lo, t_hi, t_mid, y_mid, p_eval
        integer :: iter
        integer, parameter :: max_iter = 30
        double precision, parameter :: tol = 1.0d-8
 
        if (eos%method_id /= eos_analytic) call mpistop("saha_state_from_nH_p requires eos_method=analytic")
        nh_cgs = nh_code * unit_numberdensity
        ut = unit_temperature
 
        ! brackets: y = 0 gives T = p/(nH (1+He)), y = 1 gives T = p/(nH (2+He))
        t_hi = p_code / (nh_code * (1.0d0 + eos%He_abundance))
        t_lo = p_code / (nh_code * (2.0d0 + eos%He_abundance))
        t_lo = max(t_lo, 100.0d0 / ut)
 
        do iter = 1, max_iter
            t_mid = 0.5d0 * (t_lo + t_hi)
            y_mid = saha_y_cgs(nh_cgs, t_mid * ut)
            p_eval = nh_code * (1.0d0 + eos%He_abundance + y_mid) * t_mid
 
            if (p_eval < p_code) then
                t_lo = t_mid
            else
                t_hi = t_mid
            end if
 
            if (dabs(t_hi - t_lo) < tol * t_mid) exit
        end do
 
        if (iter > max_iter .and. .not. saha_nonconv_warned .and. mype == 0) then
            saha_nonconv_warned = .true.
            write(*,*) "WARNING: saha_state_from_nH_p bisection reached max_iter ", &
                 "without converging; using last iterate. Not reported again."
        end if
 
        t_out = t_mid
        y_out = y_mid
 
        if (present(eint_nh_out)) &
            eint_nh_out = saha_eint_cgs(y_mid, t_mid * ut) * unit_numberdensity / unit_pressure
 
    end subroutine saha_state_from_nh_p
 
    !> Temperature and ionisation fraction from (nH, eint/nH) in code units, by bisection
    subroutine saha_t_from_nh_eint(nH_code, eint_nH_code, T_out, y_out)
        double precision, intent(in)  :: nh_code, eint_nh_code
        double precision, intent(out) :: t_out, y_out
 
        if (eos%method_id /= eos_analytic) call mpistop("saha_T_from_nH_eint requires eos_method=analytic")
        call saha_t_bisection(nh_code, eint_nh_code, t_out, y_out)
 
    end subroutine saha_t_from_nh_eint
 
    !> Bisection for T(nH, eint/nH) between the fully ionised (y = 1) and neutral (y = 0) limits
    subroutine saha_t_bisection(nH_code, eint_nH_code, T_out, y_out)
        double precision, intent(in)  :: nh_code, eint_nh_code
        double precision, intent(out) :: t_out, y_out
 
        double precision :: nh_cgs, eint_nh_cgs, t_lo, t_hi, t_mid
        double precision :: y_mid, cv_ion, cv_neu
        integer :: iter
        integer, parameter :: max_iter = 30
        double precision, parameter :: tol = 1.0d-8
 
        nh_cgs = nh_code * unit_numberdensity
        eint_nh_cgs = eint_nh_code * unit_pressure / unit_numberdensity
 
        ! temperature bounds from the fully ionised (y = 1) and neutral (y = 0) limits
        cv_ion = 1.5d0 * (2.0d0 + eos%He_abundance) * saha_kb_cgs
        cv_neu = 1.5d0 * (1.0d0 + eos%He_abundance) * saha_kb_cgs
        if (eos%ionE) then
            t_lo = (eint_nh_cgs - saha_chi_h_cgs) / cv_ion
        else
            t_lo = eint_nh_cgs / cv_ion
        end if
        t_lo = max(t_lo, 100.0d0)
        t_hi = max(eint_nh_cgs / cv_neu, t_lo + 1.0d0)
 
        ! eint(T) is monotonic, so the bracket always holds; where hydrogen is already fully
        ! ionised at T_lo (to 1e-10), T_lo = (eint - chi_H)/cv_ion is the solution
        if (eos%ionE) then
            if (saha_y_cgs(nh_cgs, t_lo) > 1.0d0 - 1.0d-10) then
                t_out = t_lo / unit_temperature
                y_out = 1.0d0
                return
            end if
        end if
 
        ! Bisect in CGS
        do iter = 1, max_iter
            t_mid = 0.5d0 * (t_lo + t_hi)
            y_mid = saha_y_cgs(nh_cgs, t_mid)
            if (saha_eint_cgs(y_mid, t_mid) < eint_nh_cgs) then
                t_lo = t_mid
            else
                t_hi = t_mid
            end if
 
            if (dabs(t_hi - t_lo) < tol * t_mid) exit
        end do
 
        if (iter > max_iter .and. .not. saha_nonconv_warned .and. mype == 0) then
            saha_nonconv_warned = .true.
            write(*,*) "WARNING: saha_T_bisection reached max_iter without ", &
                 "converging; using last iterate. Not reported again."
        end if
 
        t_out = t_mid / unit_temperature
        y_out = y_mid
 
    end subroutine saha_t_bisection
 
    !> Build the analytic Gamma_1 table in eos%gamma1_p on (log10 nH, log10 p/nH), the axes of the
    !> LTE tables, so that gamma1_from_nH_p serves every method. Each node takes T and y from the
    !> Saha solve at (nH, p), and Gamma_1 = chi_rho + chi_T^2 (1+He+y)/c_v, with c_v per nH in kB
    subroutine build_gamma1_analytic_table()
        integer, parameter :: ng = 256
        integer :: i, j
        double precision :: log_nh_min, log_nh_max, log_pn_min, log_pn_max, t_min, t_max
        double precision :: dx1, dx2, nh_code, p_code, t_code, y, chi, dlny, npart, cv
        double precision :: chi_t, chi_rho, g1_val, g1_min, g1_max
 
        ! log10(nH [CGS]) from 5 to 19, and p/nH = (1+He+y) T over T from 10^2.8 K neutral to
        ! 10^6.3 K ionised, all in code units
        log_nh_min = 5.0d0  - dlog10(unit_numberdensity)
        log_nh_max = 19.0d0 - dlog10(unit_numberdensity)
        t_min = 10.0d0**2.8d0 / unit_temperature
        t_max = 10.0d0**6.3d0 / unit_temperature
        log_pn_min = dlog10((1.0d0 + eos%He_abundance) * t_min)
        log_pn_max = dlog10((2.0d0 + eos%He_abundance) * t_max)
 
        eos%gamma1_p%dim1 = ng
        eos%gamma1_p%dim2 = ng
        eos%gamma1_p%var1_min = log_nh_min
        eos%gamma1_p%var1_max = log_nh_max
        eos%gamma1_p%var2_min = log_pn_min
        eos%gamma1_p%var2_max = log_pn_max
        eos%gamma1_p%filename = 'computed_gamma1_analytic'
        allocate(eos%gamma1_p%table(ng, ng))
 
        dx1 = (log_nh_max - log_nh_min) / dble(ng - 1)
        dx2 = (log_pn_max - log_pn_min) / dble(ng - 1)
        g1_min = 1.0d30
        g1_max = -1.0d30
        do j = 1, ng
            do i = 1, ng
                nh_code = 10.0d0**(log_nh_min + (i - 1) * dx1)
                p_code  = nh_code * 10.0d0**(log_pn_min + (j - 1) * dx2)
                call saha_state_from_nh_p(nh_code, p_code, t_code, y)
                chi = saha_chi_kb / (t_code * unit_temperature)
                ! dln y/dln T at fixed nH from y^2/(1-y) ~ T^1.5 exp(-chi_H/kT)/nH
                dlny = (1.0d0 - y) / (2.0d0 - y) * (1.5d0 + chi)
                npart = 1.0d0 + eos%He_abundance + y
                cv = 1.5d0 * npart + 1.5d0 * y * dlny
                if (eos%ionE) cv = cv + y * dlny * chi
                chi_t = 1.0d0 + y * dlny / npart
                chi_rho = 1.0d0 - y * (1.0d0 - y) / ((2.0d0 - y) * npart)
                g1_val = chi_rho + chi_t**2 * npart / cv
                eos%gamma1_p%table(i, j) = g1_val
                g1_min = min(g1_min, g1_val)
                g1_max = max(g1_max, g1_val)
            end do
        end do
        call precompute_step_inv(eos%gamma1_p)
 
        if (mype == 0) then
            write(*, '(A,F8.4,A,F8.4)') &
                ' Gamma1 analytic table built (nH x p/nH): min = ', g1_min, ', max = ', g1_max
        end if
 
    end subroutine build_gamma1_analytic_table
 
end module mod_eos_lte_saha
!> Needs a line after to pass the preprocessor