mod_eos_LTE.t

Go to the documentation of this file.

!=============================================================================
!> LTE (Saha-table) EoS kernels and finalise for the eos% family.
!>
!> Carved out of mod_eos.t (the per-type split). Owns everything LTE: the
!> per-substep update_eos_LTE (caches Te_/Ne_), the cached-Te/Rfactor getters,
!> the scalar table kernels (T, y, p2eint, gamma1, p/nH, eint, bisections), the
!> per-method table finalisers, and eos_finalise_LTE (the LTE arm of the
!> eos_finalise dispatcher). Depends only on the EoS leaf modules + shared
!> accessors, never on mod_eos -> no circular use. mod_eos re-exports the public
!> kernels so existing `use mod_eos` callers are unaffected.
!=============================================================================
module mod_eos_lte
    use mod_global_parameters
    use mod_eos_container
    use mod_eos_shared_functions
    use mod_eos_interp
    use mod_eos_lte_saha
    use mod_eos_lte_entropy
    use mod_eos_lte_tables
    use mod_eos_lte_state
    use mod_timing
    use mod_comm_lib, only: mpistop
 
    implicit none
    private
 
    !> update_eos_LTE, get_Te_LTE, get_temperature_from_eint_LTE and
    !> Rfactor_from_LTE are PRIVATE: all bound to eos% pointers inside
    !> eos_finalise_LTE, so no external module names them.
    !> get_temperature_from_eint_fast_LTE stays public -- the seam binds it to
    !> tc_fl%get_temperature_from_eint (fast TC path).
    public :: get_temperature_from_eint_fast_lte
    !> Scalar table kernels (re-exported by mod_eos for the hd/mhd/ffhd seams).
    !> T_and_y_from_nH_eint is internal-only (entropy update path).
    public :: y_from_nh_eint, t_from_nh_eint
    public :: p2eint_from_nh_p, gamma1_from_nh_p, p_nh_from_eint
    public :: eint_nh_from_t, eint_from_p_bisect, eos_get_eintt_grid
    !> LTE arms of the eos_init / eos_finalise dispatchers
    public :: eos_init_lte, eos_finalise_lte
 
contains
 
    !> LTE arm of eos_init (before units are known): dispatch on eos_method to
    !> the chosen method's loader, each owned by its method module (analytic loads
    !> nothing). The LTE temperature-from-pressure getter is a later pass, so
    !> unlike FI/PI it is left unset here.
    subroutine eos_init_lte()
        select case (eos%method_id)
        case (eos_analytic); call load_analytic_lte()
        case (eos_entropy);  call load_entropy_lte()
        case default;        call load_state_lte()   ! EOS_STATE
        end select
    end subroutine eos_init_lte
 
    !> LTE arm of eos_finalise: wire the LTE runtime pointer targets, then
    !> dispatch on eos_method to the per-method finaliser, each owned by its
    !> method module (analytic->mod_eos_LTE_saha, entropy->mod_eos_LTE_entropy,
    !> tables->mod_eos_LTE_tables).
    subroutine eos_finalise_lte()
        eos%update_eos => update_eos_lte
        ! eos%get_thermal_pressure is set by (m)hd_link_eos
        eos%get_temperature_from_eint => get_temperature_from_eint_lte
        eos%get_Te => get_te_lte
        !> R-factor: physics-independent (cached-Ne lookup), bound here rather
        !> than in the hd/mhd/ffhd link arms so the target stays private.
        eos%get_Rfactor => rfactor_from_lte
        select case (eos%method_id)
        case (eos_analytic); call finalise_analytic_lte()
        case (eos_entropy);  call finalise_entropy_lte()
        case default;        call finalise_state_lte()   ! EOS_STATE
        end select
    end subroutine eos_finalise_lte
 
 
    subroutine update_eos_lte(ixI^L, ixO^L, w, x)
        use mod_physics
        !> This routine is called before each RK substep
        integer, intent(in)             :: ixi^l,ixo^l
        double precision, intent(in)    :: x(ixi^s,1:ndim)
        double precision, intent(inout) :: w(ixi^s,1:nw)
 
        double precision :: wlocal(ixi^s,1:nw)
        double precision :: deltaefactor(ixi^s)
        double precision :: prev_y(ixi^s), new_y(ixi^s)
        double precision :: pth(ixi^s)
        double precision :: nh_in(ixi^s), nh(ixi^s), eint_in(ixi^s)
        double precision :: rfactor_fi
        double precision :: yy, eint_nh_floor
        double precision :: time0
        integer :: ix^d
 
        timeeos0 = mpi_wtime() !> For monitoring cost of eos module
 
        wlocal(ixi^s,1:nw)=w(ixi^s,1:nw)
        call phys_e_to_ei(ixi^l,ixo^l,wlocal,x) !> wlocal now contains internal energy, NOT total energy
 
        pth(ixo^s) = eos%gamma_minus_1 * wlocal(ixo^s,iw_e) !> pressure from internal energy only - SHOULD ONLY BE USED WHEN IonE IS UNIMPORTANT
 
        call eos%get_nH(w, x, ixi^l, ixo^l, nh)
        ! Enforce internal energy floor for EoS lookup only.
        ! Prevents NaN from dlog10(eint<0) after strong rarefactions where
        ! kinetic energy can numerically exceed total energy.
        ! The conserved energy w(iw_e) is NOT modified: the physical fluxes
        ! based on the small positive pressure will naturally restore the cell.
        if (eos%method_id == eos_analytic) then
            ! Analytic: floor to ~100 K equivalent (no tables loaded)
            ! In code units: eint/nH = T_code / (gamma-1) for neutral gas
            {do ix^db=ixomin^db,ixomax^db\}
                eint_nh_floor = nh(ix^d) * eos%inv_gamma_minus_1 * 100.0d0 / unit_temperature
                wlocal(ix^d,iw_e) = max(wlocal(ix^d,iw_e), eint_nh_floor)
            {end do\}
        else
            nh_in(ixo^s) = dlog10(nh(ixo^s))
            {do ix^db=ixomin^db,ixomax^db\}
                !> eos%T%var2_min holds the forward log(eint/nH) lower bound for
                !> all table methods (entropy sources it from Tfwd in eos_finalise).
                eint_nh_floor = nh(ix^d) * 10.0d0**eos%T%var2_min
                if (wlocal(ix^d,iw_e) < eint_nh_floor) then
                    wlocal(ix^d,iw_e) = eint_nh_floor
                end if
            {end do\}
            eint_in(ixo^s) = dlog10(wlocal(ixo^s,iw_e)) - nh_in(ixo^s)
        end if
 
        !> Constant FI Rfactor for the fully-ionised fast path.
        !> Must NOT use eos%get_Rfactor which reads iw_ne (uninitialised at IC).
        rfactor_fi = eos%n_per_nH_FI / eos%nH2rhoFactor
 
        {do ix^db=ixomin^db,ixomax^db\}
            if (wlocal(ix^d,iw_e) / w(ix^d,iw_rho) > eos%eint_rho_FI_threshold) then
                !> Fully ionised: analytical formulae, no table lookups
                new_y(ix^d) = eos%neOnH_FI
                if (eos%ionE) then
                    !> T = (eint - eion*nH) * (gamma-1) / (Rfactor_FI * rho)
                    w(ix^d,iw_te) = eos%gamma_minus_1 &
                        * (wlocal(ix^d,iw_e) - eos%eion_per_nH * nh(ix^d)) &
                        / (rfactor_fi * w(ix^d,iw_rho))
                else
                    w(ix^d,iw_te) = pth(ix^d) / (nh(ix^d) * (1.0d0 + eos%He_abundance + new_y(ix^d)))
                end if
            else
                !> Ionisation zone
                if (eos%method_id == eos_analytic) then
                    !> Analytical Saha: solve quadratic for T and y
                    call saha_t_from_nh_eint(nh(ix^d), &
                        wlocal(ix^d,iw_e) / nh(ix^d), &
                        w(ix^d,iw_te), yy)
                    new_y(ix^d) = yy
                else
                    !> Table lookups: fused T+y for shared index computation
                    if (.not. eos%ionE) then
                        new_y(ix^d) = y_from_nh_eint(nh_in(ix^d),eint_in(ix^d))
                        w(ix^d,iw_te) = pth(ix^d) / (nh(ix^d) * (1.0d0 + eos%He_abundance + new_y(ix^d)))
                    else if (eos%method_id == eos_entropy) then
                        !> Entropy method dispatches T_and_y_from_nH_eint to the
                        !> biquintic-Hermite forward at any eint; no FI fallback
                        !> needed (and eos%T isn't loaded so var2_max is 0 here).
                        call t_and_y_from_nh_eint(nh_in(ix^d), eint_in(ix^d), &
                            w(ix^d,iw_te), yy)
                        new_y(ix^d) = yy
                    else
                        if (eint_in(ix^d) < eos%T%var2_max) then
                            call t_and_y_from_nh_eint(nh_in(ix^d), eint_in(ix^d), &
                                w(ix^d,iw_te), yy)
                            new_y(ix^d) = yy
                        else
                            !> Above-table fallback with ionisation energy correction
                            new_y(ix^d) = eos%neOnH_FI
                            w(ix^d,iw_te) = eos%gamma_minus_1 &
                                * (wlocal(ix^d,iw_e) - eos%eion_per_nH * nh(ix^d)) &
                                / (rfactor_fi * w(ix^d,iw_rho))
                        end if
                    end if
                end if
            end if
        {end do\}
 
        w(ixo^s,iw_ne) = new_y(ixo^s) * nh(ixo^s)
 
        timeeos_update=timeeos_update+(mpi_wtime()-timeeos0)
 
    end subroutine update_eos_lte
 
    !> Thermodynamic getters shared verbatim by the hd/mhd seams (eos%/phys
    !> pointer targets). Only routines free of physics-module symbols live here:
    !> Rfactor_from_LTE uses the cached Ne_ field; get_gamma1_FI returns
    !> eos%gamma. The FI/LTE variants that read primitive indices (rho_,p_) or
    !> the physics RR stay in each seam.
    !> Rfactor from LTE EOS: R = (1 + He + ne/nH) / nH2rhoFactor
    !> Uses stored Ne_ wextra field from update_eos_LTE.
    subroutine rfactor_from_lte(w,x,ixI^L,ixO^L,Rfactor)
        use mod_global_parameters
        integer, intent(in) :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s,1:nw)
        double precision, intent(in) :: x(ixi^s,1:ndim)
        double precision, intent(out):: rfactor(ixi^s)
 
        double precision :: y_nh
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            y_nh = w(ix^d, iw_ne) * eos%nH2rhoFactor / w(ix^d, iw_rho)
            rfactor(ix^d) = (1.0d0 + eos%He_abundance + y_nh) / eos%nH2rhoFactor
        {end do\}
 
    end subroutine rfactor_from_lte
 
    !>#######################################################################
    !> Temperature getters: cached Te accessors and temperature from the
    !> energy/pressure variable. The from-eint variants are split by eos_type
    !> (needed by the STS path: dt sees conserved, set_source sees eint -- see
    !> mod_thermal_conduction).
    !>#######################################################################
 
    subroutine get_te_lte(w,x,ixI^L,ixO^L,T)
        use mod_global_parameters
        integer, intent(in)           :: ixi^l, ixo^l
        double precision, intent(in)  :: w(ixi^s,1:nw)
        double precision, intent(in)  :: x(ixi^s,1:ndim)
        double precision, intent(out) :: t(ixi^s)
 
        t(ixo^s) = w(ixo^s,iw_te)
 
    end subroutine get_te_lte
 
    subroutine get_temperature_from_eint_lte(w, x, ixI^L, ixO^L, res)
        !> Assumes input energy is internal energy.
        !> Includes FI bypass: cells with eint/rho above threshold skip table lookups.
        use mod_physics
        integer, intent(in)             :: ixi^l,ixo^l
        double precision, intent(in)    :: x(ixi^s,1:ndim)
        double precision, intent(in) :: w(ixi^s,1:nw)
        double precision, intent(out)   :: res(ixi^s)
 
        double precision :: nh(ixi^s),nh_in(ixi^s), eint_in(ixi^s)
        double precision :: rfactor(ixi^s), rfactor_fi, t_loc, y_loc
        integer :: ix^d
 
        timeeos0 = mpi_wtime()
 
        rfactor_fi = eos%n_per_nH_FI / (1.0d0 + 4.0d0*eos%He_abundance)
 
        call eos%get_nH(w, x, ixi^l, ixo^l, nh)
        nh_in(ixo^s) = dlog10(nh(ixo^s))
        eint_in(ixo^s) = dlog10(w(ixo^s,iw_e)) - nh_in(ixo^s)
 
        {do ix^db=ixomin^db,ixomax^db\}
            if (w(ix^d,iw_e) / w(ix^d,iw_rho) > eos%eint_rho_FI_threshold) then
                !> Fully ionised: T = (eint - eion*nH) * (gamma-1) / (Rfactor_FI * rho)
                res(ix^d) = eos%gamma_minus_1 &
                    * (w(ix^d,iw_e) - eos%eion_per_nH * nh(ix^d)) &
                    / (rfactor_fi * w(ix^d,iw_rho))
            else if (eos%method_id == eos_analytic) then
                call saha_t_from_nh_eint(nh(ix^d), &
                    w(ix^d,iw_e) / nh(ix^d), t_loc, y_loc)
                res(ix^d) = t_loc
            else
                res(ix^d) = t_from_nh_eint(nh_in(ix^d),eint_in(ix^d))
            endif
        {end do\}
 
        timeeos_tfromei=timeeos_tfromei+(mpi_wtime()-timeeos0)
    end subroutine get_temperature_from_eint_lte
 
    subroutine get_temperature_from_eint_fast_lte(w, x, ixI^L, ixO^L, res)
        !> Fast TC variant: two-pass regime-aware bypass.
        !>
        !> Pass 1 (vectorised): compute FI formula for ALL cells as array ops.
        !>   The compiler vectorises this with SVML (no branches, pure arithmetic).
        !>
        !> Pass 2 (scalar): overwrite ionisation-zone cells with
        !>   bilinear table lookup + dexp.  Threshold check uses multiply
        !>   (eint <= threshold * rho) to avoid a 14-cycle scalar division.
        use mod_physics
        use mod_eos_lte_entropy, only: entropy_t_from_nh_eint
        integer, intent(in)             :: ixi^l,ixo^l
        double precision, intent(in)    :: x(ixi^s,1:ndim)
        double precision, intent(in)    :: w(ixi^s,1:nw)
        double precision, intent(out)   :: res(ixi^s)
 
        double precision :: inv_rfactor_fi, eion_rho_inv
        double precision :: log_nh_val, log_eint_nh_val
        double precision :: fx, fy, rx, ry, nh_loc, t_loc, y_loc
        integer :: jx, jy, jx1, jy1, ix^d
        double precision, parameter :: ln10 = 2.302585092994046d0
 
        if (eos%type_id /= eos_type_lte) call mpistop("get_temperature_from_eint_fast_LTE called outside its eos_type (LTE)")
        timeeos0 = mpi_wtime()
 
        !> Precompute scalar constants (avoid per-cell divisions)
        inv_rfactor_fi = (1.0d0 + 4.0d0*eos%He_abundance) / eos%n_per_nH_FI
        eion_rho_inv = eos%eion_per_nH / eos%nH2rhoFactor
 
        !> Pass 1: FI formula for ALL cells (vectorisable array operations).
        !> T = (gamma-1) * (eint - eion_per_rho * rho) / (Rfactor_FI * rho)
        res(ixo^s) = eos%gamma_minus_1 * inv_rfactor_fi &
            * (w(ixo^s,iw_e) - eion_rho_inv * w(ixo^s,iw_rho)) &
            / w(ixo^s,iw_rho)
 
        !> Pass 2: overwrite ionisation-zone cells.
        {do ix^db=ixomin^db,ixomax^db\}
            if (w(ix^d,iw_e) <= eos%eint_rho_FI_threshold * w(ix^d,iw_rho)) then
                if (eos%method_id == eos_analytic) then
                    nh_loc = w(ix^d, iw_rho) / eos%nH2rhoFactor
                    call saha_t_from_nh_eint(nh_loc, &
                        w(ix^d,iw_e) / nh_loc, t_loc, y_loc)
                    res(ix^d) = t_loc
                else if (eos%method_id == eos_entropy) then
                    if (iw_log_nh > 0) then
                        log_nh_val = block%wextra(ix^d, iw_log_nh)
                    else
                        log_nh_val = dlog10(w(ix^d, iw_rho) / eos%nH2rhoFactor)
                    end if
                    log_eint_nh_val = dlog10(w(ix^d, iw_e)) - log_nh_val
                    res(ix^d) = entropy_t_from_nh_eint(eos%Tfwd, eos%Tfwd_x, &
                        eos%Tfwd_y, eos%Tfwd_xy, log_nh_val, log_eint_nh_val)
                else
                    if (iw_log_nh > 0) then
                        log_nh_val = block%wextra(ix^d, iw_log_nh)
                    else
                        log_nh_val = dlog10(w(ix^d, iw_rho) / eos%nH2rhoFactor)
                    end if
                    log_eint_nh_val = dlog10(w(ix^d, iw_e)) - log_nh_val
 
                    if (eos%T%is_uniform) then
                        !> Inlined uniform bilinear (hot path: no function call).
                        ry = max(0.0d0, min((log_nh_val - eos%T%var1_min) * eos%T%step_inv_1, &
                                             dble(eos%T%dim1-1)))
                        rx = max(0.0d0, min((log_eint_nh_val - eos%T%var2_min) * eos%T%step_inv_2, &
                                             dble(eos%T%dim2-1)))
                        jy = int(ry); jx = int(rx)
                        jy1 = min(jy+1, eos%T%dim1-1)
                        jx1 = min(jx+1, eos%T%dim2-1)
                        fy = ry - dble(jy); fx = rx - dble(jx)
                        res(ix^d) = dexp(ln10 * ( &
                            (1.0d0-fy)*((1.0d0-fx)*eos%T%table(jy+1,jx+1)  &
                                              + fx *eos%T%table(jy+1,jx1+1)) &
                           +      fy *((1.0d0-fx)*eos%T%table(jy1+1,jx+1)   &
                                              + fx *eos%T%table(jy1+1,jx1+1))))
                    else
                        !> Non-uniform: dispatch to the binary-search bilinear.
                        !> Slightly slower per cell (binary search adds ~8 cmps per axis)
                        !> but unavoidable when the grid is non-uniform.
                        res(ix^d) = dexp(ln10 * bilinear_lookup( &
                            log_nh_val, log_eint_nh_val, eos%T))
                    end if
                end if
            end if
        {end do\}
 
        timeeos_tfromei = timeeos_tfromei + (mpi_wtime()-timeeos0)
    end subroutine get_temperature_from_eint_fast_lte
 
    !>#######################################################################
    !> Scalar EoS kernels: pointwise lookups on (log10 nH, log10 eint/nH) or
    !> (log10 nH, log10 p/nH) in code units, called by the hd/mhd EoS layer.
    !> Each dispatches analytic / entropy / table. Ionisation fraction and
    !> temperature first, then pressure / Gamma_1 / eint.
    !>#######################################################################
 
    !> Ionization fraction from (log10 nH, log10 eint/nH) in code units.
    !> Dispatches: analytic -> Saha quadratic, tables -> PCHIP interpolation.
    double precision function y_from_nh_eint(nH, eint_nh) result(result_val)
        use mod_eos_lte_entropy, only: entropy_y_from_nh_eint
        double precision, intent(in) :: nh, eint_nh
        double precision, parameter :: ln10 = 2.302585092994046d0
        double precision :: t_loc, y_loc, eint_rho
 
        if (eos%type_id /= eos_type_lte) call mpistop("y_from_nH_eint called outside its eos_type (LTE)")
        if (eos%method_id == eos_analytic) then
            ! FI bypass: skip Saha solve for fully ionized cells
            eint_rho = 10.0d0**eint_nh / eos%nH2rhoFactor
            if (eint_rho > eos%eint_rho_FI_threshold) then
                result_val = eos%neOnH_FI
                return
            end if
            call saha_t_from_nh_eint(10.0d0**nh, 10.0d0**eint_nh, t_loc, y_loc)
            result_val = y_loc
        else if (eos%method_id == eos_entropy) then
            ! neOnH stored linearly (not log10); single bicubic Hermite lookup.
            result_val = entropy_y_from_nh_eint(eos%neOnH, eos%neOnH_x, &
                eos%neOnH_y, eos%neOnH_xy, nh, eint_nh)
        else
            result_val = dexp(ln10 * bicubic_lookup(nh, eint_nh, eos%neOnH))
        end if
    end function y_from_nh_eint
 
    !> Temperature from (log10 nH, log10 eint/nH) in code units.
    !> Dispatches: analytic -> Saha bisection/Newton, tables -> PCHIP interpolation.
    double precision function t_from_nh_eint(nH, eint_nh) result(result_val)
        use mod_eos_lte_entropy, only: entropy_t_from_nh_eint
        double precision, intent(in) :: nh, eint_nh
        double precision, parameter :: ln10 = 2.302585092994046d0
        double precision :: t_loc, y_loc, eint_rho, rfactor_fi
 
        if (eos%type_id /= eos_type_lte) call mpistop("T_from_nH_eint called outside its eos_type (LTE)")
        if (eos%method_id == eos_analytic) then
            ! FI bypass: skip Saha solve for fully ionized cells
            eint_rho = 10.0d0**eint_nh / eos%nH2rhoFactor
            if (eint_rho > eos%eint_rho_FI_threshold) then
                rfactor_fi = eos%n_per_nH_FI / (1.0d0 + 4.0d0*eos%He_abundance)
                result_val = eos%gamma_minus_1 &
                    * (10.0d0**eint_nh - eos%eion_per_nH) / rfactor_fi
                return
            end if
            call saha_t_from_nh_eint(10.0d0**nh, 10.0d0**eint_nh, t_loc, y_loc)
            result_val = t_loc
        else if (eos%method_id == eos_entropy) then
            result_val = entropy_t_from_nh_eint(eos%Tfwd, eos%Tfwd_x, &
                eos%Tfwd_y, eos%Tfwd_xy, nh, eint_nh)
        else
            result_val = dexp(ln10 * bicubic_lookup(nh, eint_nh, eos%T))
        end if
    end function t_from_nh_eint
 
    !> Fused T+y lookup from (log10 nH, log10 eint/nH) in code units.
    !> Computes grid indices once, evaluates both T and y tables.
    !> Saves one index computation + better cache utilisation vs separate calls.
    subroutine t_and_y_from_nh_eint(log_nH, log_eint_nH, T_out, y_out)
        use mod_eos_lte_entropy, only: entropy_t_and_y_from_nh_eint
        double precision, intent(in) :: log_nh, log_eint_nh
        double precision, intent(out) :: t_out, y_out
        double precision, parameter :: ln10 = 2.302585092994046d0
        double precision :: results(3)
 
        if (eos%method_id == eos_analytic) then
            call saha_t_from_nh_eint(10.0d0**log_nh, 10.0d0**log_eint_nh, t_out, y_out)
        else if (eos%method_id == eos_entropy) then
            call entropy_t_and_y_from_nh_eint(eos%Tfwd, eos%Tfwd_x, &
                eos%Tfwd_y, eos%Tfwd_xy, &
                eos%neOnH, eos%neOnH_x, eos%neOnH_y, eos%neOnH_xy, &
                log_nh, log_eint_nh, t_out, y_out)
        else if (allocated(eos%table_eint_il)) then
            ! Interleaved PCHIP fast-path.
            ! (T, n_e/n_H, p/n_H), so cache behaviour is preserved.
            if (eos%T%is_uniform) then
                call interp_pchip_interleaved(log_nh, log_eint_nh, &
                    eos%table_eint_il, 3, eos%T%dim2, eos%T%dim1, &
                    eos%T%var1_min, eos%T%var1_max, &
                    eos%T%var2_min, eos%T%var2_max, results)
            else
                call interp_pchip_interleaved_nu(log_nh, log_eint_nh, &
                    eos%table_eint_il, 3, eos%T%dim2, eos%T%dim1, &
                    eos%T%var1_nodes, eos%T%var2_nodes, &
                    eos%T%guard_1, eos%T%guard_M_1, eos%T%guard_scale_1, &
                    eos%T%guard_2, eos%T%guard_M_2, eos%T%guard_scale_2, &
                    results)
            end if
            t_out = dexp(ln10 * results(1))
            y_out = dexp(ln10 * results(2))
        else
            ! No interleaved table built -- separate dispatcher calls
            t_out = dexp(ln10 * bicubic_lookup(log_nh, log_eint_nh, eos%T))
            y_out = dexp(ln10 * bicubic_lookup(log_nh, log_eint_nh, eos%neOnH))
        end if
    end subroutine t_and_y_from_nh_eint
 
    !> Pressure-to-eint ratio from (log10 nH, log10 p/nH) in code units.
    !> Dispatches: analytic -> Saha solve for eint/p, tables -> PCHIP interpolation.
    double precision function p2eint_from_nh_p(nH, ponH) result(result_val)
        use mod_eos_lte_entropy, only: entropy_eint_from_nh_p
        double precision, intent(in) :: nh, ponh
        double precision :: nh_code, p_code, t_loc, y_loc, eint_nh_loc, p_rho
 
        if (eos%type_id /= eos_type_lte) call mpistop("p2eint_from_nH_p called outside its eos_type (LTE)")
        if (eos%method_id == eos_analytic) then
            nh_code = 10.0d0**nh
            p_code = nh_code * 10.0d0**ponh
            ! FI bypass: skip Saha solve for fully ionized cells
            p_rho = p_code / (nh_code * eos%nH2rhoFactor)
            if (p_rho > eos%p_rho_FI_threshold) then
                result_val = eos%inv_gamma_minus_1 &
                    + eos%eion_per_nH * nh_code / p_code
                return
            end if
            call saha_state_from_nh_p(nh_code, p_code, t_loc, y_loc, eint_nh_loc)
            result_val = eint_nh_loc * nh_code / p_code
        else if (eos%method_id == eos_entropy) then
            ! One bicubic-Hermite lookup on the eintP table.
            result_val = entropy_eint_from_nh_p(eos%eintP, eos%eintP_x, &
                eos%eintP_y, eos%eintP_xy, nh, ponh)
        else
            result_val = bicubic_lookup(nh, ponh, eos%p2eint)
        end if
    end function p2eint_from_nh_p
 
    !> Gamma_1 from pressure-indexed table: (log10 nH, log10 p/nH) -> Gamma_1.
    !> For 'entropy' the conversion p -> eint -> Gamma_1 via formula is intended
    !> to keep Maxwell consistency with the forward at runtime. The p->eint inverse 
    !> is one bisection per cell -- non-trivial cost; this function is in the hot path
    !> via hd_get_csound2_LTE.
    double precision function gamma1_from_nh_p(log_nH, log_p_nH) result(g1)
        use mod_eos_lte_entropy, only: entropy_gamma1_from_nh_p
        double precision, intent(in) :: log_nh, log_p_nh
        if (eos%type_id /= eos_type_lte) call mpistop("gamma1_from_nH_p called outside its eos_type (LTE)")
        if (eos%method_id == eos_entropy) then
            ! Single bicubic Hermite lookup on the pre-built Gamma_1(rho, p) table.
            ! ZERO runtime iterations; no p->eint intermediate inversion.
            g1 = entropy_gamma1_from_nh_p(eos%g1p, eos%g1p_x, eos%g1p_y, &
                eos%g1p_xy, log_nh, log_p_nh)
        else
            g1 = bicubic_lookup(log_nh, log_p_nh, eos%gamma1_p)
        end if
    end function gamma1_from_nh_p
 
    !> Merged log10(p/nH) lookup: (log10 nH, log10 eint/nH) -> log10(p/nH)
    !> Single PCHIP evaluation replacing separate T + neOnH lookups.
    double precision function log_p_from_nh_eint(log_nH, log_eint_nH) result(lp)
        double precision, intent(in) :: log_nh, log_eint_nh
        lp = bicubic_lookup(log_nh, log_eint_nh, eos%log_p)
    end function log_p_from_nh_eint
 
    !> p/nH from (log10 nH, log10 eint/nH) in code units.
    !> Returns (1+He+y)*T directly -- single lookup replaces T + y lookups.
    double precision function p_nh_from_eint(log_nH, log_eint_nH) result(p_nH)
        use mod_eos_lte_entropy, only: entropy_p_nh_from_eint
        double precision, intent(in) :: log_nh, log_eint_nh
        double precision, parameter :: ln10 = 2.302585092994046d0
        double precision :: t_loc, y_loc
 
        if (eos%type_id /= eos_type_lte) call mpistop("p_nH_from_eint called outside its eos_type (LTE)")
        if (eos%method_id == eos_analytic) then
            call saha_t_from_nh_eint(10.0d0**log_nh, 10.0d0**log_eint_nh, t_loc, y_loc)
            p_nh = (1.0d0 + eos%He_abundance + y_loc) * t_loc
        else if (eos%method_id == eos_entropy) then
            p_nh = entropy_p_nh_from_eint(eos%pfwd, eos%pfwd_x, eos%pfwd_y, &
                                           eos%pfwd_xy, log_nh, log_eint_nh)
        else
            p_nh = dexp(ln10 * bicubic_lookup(log_nh, log_eint_nh, eos%p_over_nH))
        end if
    end function p_nh_from_eint
 
    !> Internal energy per nH from (log10 nH, log10 T) in code units.
    !> Uses the bisection-built inverse table (H+He, machine precision).
    !> Fallback: H-only Saha if table not built.
    double precision function eint_nh_from_t(log_nH, log_T) result(eint_nH)
        use mod_eos_lte_entropy, only: entropy_eint_from_nh_t
        double precision, intent(in) :: log_nh, log_t
        double precision, parameter :: ln10 = 2.302585092994046d0
        double precision :: log_e_nh
 
        if (eos%type_id /= eos_type_lte) call mpistop("eint_nH_from_T called outside its eos_type (LTE)")
        if (eos%method_id == eos_entropy) then
            ! Single bicubic Hermite lookup on the pre-built eint(rho, T)
            ! inverse table. ZERO runtime iterations.
            log_e_nh = entropy_eint_from_nh_t(eos%eintT, eos%eintT_x, &
                eos%eintT_y, eos%eintT_xy, log_nh, log_t)
            eint_nh = 10.0d0**log_e_nh
        else if (allocated(eos%eint_from_T%table)) then
            eint_nh = dexp(ln10 * bicubic_lookup(log_nh, log_t, eos%eint_from_T))
        else
            eint_nh = saha_eint_from_nh_t(10.0d0**log_nh, 10.0d0**log_t)
        end if
    end function eint_nh_from_t
 
    !>#######################################################################
    !> Iterative solvers: table-guessed bisection for the WB pressure->eint
    !> inversion since the method needs high accuracy to not self-seed perturbations.
    !>#######################################################################
 
    !> Cached bisection on the log_p table for WB p->eint inversion.
    !> Precomputes nH-direction indices and table values once, then
    !> bisects using only cheap PCHIP evaluations with varying tx.
    !> Expects a narrow initial bracket [lo, hi] (e.g. from p2eint guess).
    !>
    !> Adaptive-grid support: when eos%log_p%is_uniform is .false., the
    !> y- and x-direction index calculations switch to binary search on
    !> the explicit node arrays. The cached 4x4 stencil mechanism is
    !> unchanged; only the index-to-cell mapping changes.
    subroutine log_p_bisect_cached(log_nH, log_p_target, &
        log_eint_lo, log_eint_hi, max_iter, log_eint_result)
        use mod_lookup_table, only: find_index_bsearch
        double precision, intent(in)    :: log_nh, log_p_target
        double precision, intent(inout) :: log_eint_lo, log_eint_hi
        integer, intent(in)             :: max_iter
        double precision, intent(out)   :: log_eint_result
 
        integer :: nx, ny, ix, iy, iter, ii
        integer :: i0, i1, i2, i3, j0, j1, j2, j3
        double precision :: tx, ty, rx, ry
        double precision :: xstep_inv, ystep_inv
        double precision :: vmin_x, vmax_x, vmin_y, vmax_y
        double precision :: log_eint_mid, log_p_eval
        logical :: is_unif
 
        !> Table values: tv(y_row, x_col) for 4 y-rows x 4 x-cols
        double precision :: tv(4,4)
        !> Cached ix for detecting grid cell change
        integer :: ix_cached
 
        nx = eos%log_p%dim2   ! eint axis (varx)
        ny = eos%log_p%dim1   ! nH axis (vary)
        is_unif = eos%log_p%is_uniform
 
        if (is_unif) then
            vmin_x = eos%log_p%var2_min
            vmax_x = eos%log_p%var2_max
            vmin_y = eos%log_p%var1_min
            vmax_y = eos%log_p%var1_max
            xstep_inv = dble(nx-1) / (vmax_x - vmin_x)
            ystep_inv = dble(ny-1) / (vmax_y - vmin_y)
        end if
 
        !> Precompute nH indices (constant across all iterations)
        if (is_unif) then
            ry = (log_nh - vmin_y) * ystep_inv
            ry = max(0.0d0, min(ry, dble(ny-1)))
            iy = int(ry)
            ty = ry - dble(iy)
        else
            !> Non-uniform: binary search on var1_nodes (length ny = dim1)
            if (log_nh <= eos%log_p%var1_nodes(1)) then
                iy = 0; ty = 0.0d0
            else if (log_nh >= eos%log_p%var1_nodes(ny)) then
                iy = ny - 2; ty = 1.0d0
            else
                ii = find_index_guard(eos%log_p%var1_nodes, ny, log_nh, &
                    eos%log_p%guard_1, eos%log_p%guard_M_1, eos%log_p%guard_scale_1)
                iy = max(0, min(ii - 2, ny - 2))
                ty = (log_nh - eos%log_p%var1_nodes(iy+1)) &
                   / (eos%log_p%var1_nodes(iy+2) - eos%log_p%var1_nodes(iy+1))
                ty = max(0.0d0, min(ty, 1.0d0))
            end if
        end if
 
        !> Clamped y-row indices
        j0 = max(0, min(ny-1, iy-1))
        j1 = max(0, min(ny-1, iy  ))
        j2 = max(0, min(ny-1, iy+1))
        j3 = max(0, min(ny-1, iy+2))
 
        ix_cached = -1
 
        do iter = 1, max_iter
            log_eint_mid = 0.5d0 * (log_eint_lo + log_eint_hi)
 
            !> Compute eint-direction index
            if (is_unif) then
                rx = (log_eint_mid - vmin_x) * xstep_inv
                rx = max(0.0d0, min(rx, dble(nx-1)))
                ix = int(rx)
                tx = rx - dble(ix)
            else
                !> Non-uniform: binary search on var2_nodes (length nx = dim2)
                if (log_eint_mid <= eos%log_p%var2_nodes(1)) then
                    ix = 0; tx = 0.0d0
                else if (log_eint_mid >= eos%log_p%var2_nodes(nx)) then
                    ix = nx - 2; tx = 1.0d0
                else
                    ii = find_index_guard(eos%log_p%var2_nodes, nx, log_eint_mid, &
                        eos%log_p%guard_2, eos%log_p%guard_M_2, eos%log_p%guard_scale_2)
                    ix = max(0, min(ii - 2, nx - 2))
                    tx = (log_eint_mid - eos%log_p%var2_nodes(ix+1)) &
                       / (eos%log_p%var2_nodes(ix+2) - eos%log_p%var2_nodes(ix+1))
                    tx = max(0.0d0, min(tx, 1.0d0))
                end if
            end if
 
            !> Load 16 table values only when ix changes
            if (ix /= ix_cached) then
                i0 = max(0, min(nx-1, ix-1))
                i1 = max(0, min(nx-1, ix  ))
                i2 = max(0, min(nx-1, ix+1))
                i3 = max(0, min(nx-1, ix+2))
 
                tv(1,1) = eos%log_p%table(j0+1, i0+1)
                tv(1,2) = eos%log_p%table(j0+1, i1+1)
                tv(1,3) = eos%log_p%table(j0+1, i2+1)
                tv(1,4) = eos%log_p%table(j0+1, i3+1)
                tv(2,1) = eos%log_p%table(j1+1, i0+1)
                tv(2,2) = eos%log_p%table(j1+1, i1+1)
                tv(2,3) = eos%log_p%table(j1+1, i2+1)
                tv(2,4) = eos%log_p%table(j1+1, i3+1)
                tv(3,1) = eos%log_p%table(j2+1, i0+1)
                tv(3,2) = eos%log_p%table(j2+1, i1+1)
                tv(3,3) = eos%log_p%table(j2+1, i2+1)
                tv(3,4) = eos%log_p%table(j2+1, i3+1)
                tv(4,1) = eos%log_p%table(j3+1, i0+1)
                tv(4,2) = eos%log_p%table(j3+1, i1+1)
                tv(4,3) = eos%log_p%table(j3+1, i2+1)
                tv(4,4) = eos%log_p%table(j3+1, i3+1)
 
                ix_cached = ix
            end if
 
            !> Evaluate PCHIP: 4 x-rows then 1 y-interp
            log_p_eval = pchip_2d_from_cache(tv, tx, ty)
 
            if (log_p_eval < log_p_target) then
                log_eint_lo = log_eint_mid
            else
                log_eint_hi = log_eint_mid
            end if
 
            if (dabs(log_eint_hi - log_eint_lo) < 1.0d-14) exit
        end do
 
        log_eint_result = 0.5d0 * (log_eint_lo + log_eint_hi)
 
    contains
 
        pure double precision function pchip_1d(p0, p1, p2, p3, t) result(v)
            double precision, intent(in) :: p0, p1, p2, p3, t
            double precision :: d0, d1, d2, m1, m2, s, a1, a2, lim
            double precision :: tt, ttt, h00, h10, h01, h11
 
            d0 = p1 - p0;  d1 = p2 - p1;  d2 = p3 - p2
 
            if (d1 == 0.0d0) then
                m1 = 0.0d0;  m2 = 0.0d0
            else
                if (d0*d1 <= 0.0d0) then
                    m1 = 0.0d0
                else
                    m1 = 2.0d0*d0*d1 / (d0 + d1)
                end if
                if (d1*d2 <= 0.0d0) then
                    m2 = 0.0d0
                else
                    m2 = 2.0d0*d1*d2 / (d1 + d2)
                end if
                s = sign(1.0d0, d1)
                a1 = s*m1;  a2 = s*m2
                if (a1 < 0.0d0) a1 = 0.0d0
                if (a2 < 0.0d0) a2 = 0.0d0
                lim = 3.0d0*abs(d1)
                if (a1 > lim) a1 = lim
                if (a2 > lim) a2 = lim
                m1 = s*a1;  m2 = s*a2
            end if
 
            tt = t*t;  ttt = tt*t
            h00 = 2.0d0*ttt - 3.0d0*tt + 1.0d0
            h10 = ttt - 2.0d0*tt + t
            h01 = -2.0d0*ttt + 3.0d0*tt
            h11 = ttt - tt
            v = h00*p1 + h10*m1 + h01*p2 + h11*m2
        end function pchip_1d
 
        pure double precision function pchip_2d_from_cache(c, tx, ty) result(z)
            double precision, intent(in) :: c(4,4), tx, ty
            double precision :: g0, g1, g2, g3
            !> 4 x-direction interpolations (one per y-row)
            g0 = pchip_1d(c(1,1), c(1,2), c(1,3), c(1,4), tx)
            g1 = pchip_1d(c(2,1), c(2,2), c(2,3), c(2,4), tx)
            g2 = pchip_1d(c(3,1), c(3,2), c(3,3), c(3,4), tx)
            g3 = pchip_1d(c(4,1), c(4,2), c(4,3), c(4,4), tx)
            !> 1 y-direction interpolation
            z = pchip_1d(g0, g1, g2, g3, ty)
        end function pchip_2d_from_cache
 
    end subroutine log_p_bisect_cached
 
    !> Given log10(nH) and log10(p), find log10(eint/nH) by table-guessed
    !> bisection on the forward pressure table.  Wraps the bracketing logic
    !> and log_p_bisect_cached into a single call for use by both HD and MHD
    !> from_conserved routines.
    subroutine eint_from_p_bisect(log_nH_val, log_p_val, log_eint_nH_out)
        double precision, intent(in)  :: log_nh_val, log_p_val
        double precision, intent(out) :: log_eint_nh_out
 
        double precision :: log_p_target, p2eint_ratio, log_eint_guess
        double precision :: log_p_at_guess, log_eint_lo, log_eint_hi
        double precision :: f_bracket, margin
        integer :: max_iter
 
        if (eos%type_id /= eos_type_lte) call mpistop("eint_from_p_bisect called outside its eos_type (LTE)")
        log_p_target = log_p_val - log_nh_val
 
        ! Initial guess from p2eint table
        p2eint_ratio = p2eint_from_nh_p(log_nh_val, log_p_target)
        log_eint_guess = dlog10(p2eint_ratio) + log_p_target
 
        log_eint_guess = max(log_eint_guess, eos%log_p%var2_min)
        log_eint_guess = min(log_eint_guess, eos%log_p%var2_max)
 
        log_p_at_guess = log_p_from_nh_eint(log_nh_val, log_eint_guess)
 
        ! Establish bracket around the guess
        margin = 5.0d-4
        if (log_p_at_guess < log_p_target) then
            log_eint_lo = log_eint_guess
            log_eint_hi = min(log_eint_guess + margin, eos%log_p%var2_max)
            f_bracket = log_p_from_nh_eint(log_nh_val, log_eint_hi) - log_p_target
        else
            log_eint_lo = max(log_eint_guess - margin, eos%log_p%var2_min)
            log_eint_hi = log_eint_guess
            f_bracket = -(log_p_from_nh_eint(log_nh_val, log_eint_lo) - log_p_target)
        end if
 
        if (f_bracket >= 0.0d0) then
            max_iter = 8
        else
            log_eint_lo = eos%log_p%var2_min
            log_eint_hi = eos%log_p%var2_max
            max_iter = 20
        end if
 
        call log_p_bisect_cached(log_nh_val, log_p_target, &
            log_eint_lo, log_eint_hi, max_iter, log_eint_nh_out)
 
    end subroutine eint_from_p_bisect
 
    !>#######################################################################
    !> Table-metadata helpers: expose inverse-table grid bounds without
    !> reaching into EoS table internals.
    !>#######################################################################
 
    !> log_nH grid metadata of the (log_nH, log_T) inverse table (eint from T),
    !> choosing the container by method ('tables'->eint_from_T, 'entropy'->eintT).
    !> n_nH=0 if no such table (analytic/FI). Lets cooling's build_Y_mod_table
    !> get its grid without reaching into EoS table internals.
    subroutine eos_get_eintt_grid(n_nH, lg_nH_min, lg_nH_max)
        integer, intent(out)          :: n_nh
        double precision, intent(out) :: lg_nh_min, lg_nh_max
 
        if (eos%type_id /= eos_type_lte) call mpistop("eos_get_eintT_grid called outside its eos_type (LTE)")
        if (allocated(eos%eint_from_T%table)) then
            n_nh = eos%eint_from_T%dim1
            lg_nh_min = eos%eint_from_T%var1_min
            lg_nh_max = eos%eint_from_T%var1_max
        else if (allocated(eos%eintT%table)) then
            n_nh = eos%eintT%dim1
            lg_nh_min = eos%eintT%var1_min
            lg_nh_max = eos%eintT%var1_max
        else
            n_nh = 0
            lg_nh_min = 0.0d0
            lg_nh_max = 0.0d0
        end if
    end subroutine eos_get_eintt_grid
 
end module mod_eos_lte