mod_hd_eos.t

Go to the documentation of this file.

!=============================================================================
!> HD <-> EoS seam: binds the eos% authority into hydrodynamics.
!>
!> hd_link_eos wires the eos%/phys_ procedure pointers (conversions, pthermal,
!> csound2/gamma1, get_Rfactor, temperature, update_temperature) per eos_type
!> (FI / LTE / PI, with PI energy- and prominence-mode overrides).
!> bind_eos_to_source wires the thermal-conduction / radiative-cooling / FLD
!> fluid-port callbacks from eos%.
!>
!> Holds the HD block routines that wrap the shared thermodynamics with HD's
!> mechanical-energy bookkeeping -- total-energy ("origin") formulation,
!> dust-aware: FI/LTE/PI conversions, p_to_e, pthermal, csound2/gamma1, and the
!> PI energy + prominence R-factor / temperature kernels.
!=============================================================================
module mod_hd_eos
    use mod_global_parameters
    use mod_physics
    use mod_eos
    !> Mode-specific kernels come from their sub-modules (the facade no longer
    !> re-exports them); each mpistops if called under the wrong eos_type/method.
    use mod_eos_lte
    use mod_eos_lte_entropy, only: entropy_eint_from_p_bisect
    use mod_eos_lte_saha
    use mod_eos_pi
    use mod_eos_container
    use mod_hd_phys
    use mod_timing
    use mod_radiative_cooling, only: build_y_mod_table
 
    use mod_comm_lib, only: mpistop
 
    implicit none
    private
 
    procedure(sub_convert), pointer, public :: hd_to_primitive  => null()
    procedure(sub_convert), pointer, public :: hd_to_conserved  => null()
 
    public :: hd_link_eos, hd_get_pthermal
 
contains
 
    !> Link the appropriate EOS conversion routines based on the selected EoS type
    subroutine hd_link_eos()
        use mod_usr_methods
 
        ! PI (partial approximations) takes the FI ideal-gas conversions as
        ! its BASE (e_int=p/(gamma-1); variable R via get_Rfactor). Energy
        ! mode (ionE) overrides them with the PI-energy variants at the end
        ! of this routine; FI conversions stay pristine.
        if (eos%eos_type == 'FI' .or. eos%eos_type == 'PI') then
            eos%to_conserved => hd_to_conserved_origin
            eos%to_primitive => hd_to_primitive_origin
        else if (eos%eos_type == 'LTE') then
            eos%to_conserved => hd_to_conserved_lte
            eos%to_primitive => hd_to_primitive_lte
        else
            call mpistop('Error: Unknown HD EOS type: ' // trim(eos%eos_type))
        end if
 
        phys_to_primitive => eos%to_primitive
        phys_to_conserved => eos%to_conserved
        phys_get_rho      => eos%get_rho
        phys_get_pthermal => hd_get_pthermal
        eos%get_thermal_pressure => hd_get_pthermal
        phys_bind_eos_to_source  => bind_eos_to_source
 
        eos%p_to_e => p_to_e !> suitable for both FI and LTE
 
        ! Link sound speed and gamma1 computation
        if (eos%eos_type == 'LTE' .and. eos%ionE) then
            eos%get_csound2 => hd_get_csound2_lte
            phys_get_gamma1 => hd_get_gamma1_lte
            ! EoS-aware prolongation: interpolate in (rho, v, T) space
            phys_to_prolong   => hd_to_prolong_lte
            phys_from_prolong => hd_from_prolong_lte
        else
            eos%get_csound2 => hd_get_csound2_fi
            phys_get_gamma1 => get_gamma1_fi
        end if
 
        ! Rfactor: only the FI case is physics-dependent (the usr_Rfactor user
        ! hook). LTE and PI bind pure eos% routines in eos_finalise_{LTE,PI},
        ! which run after this and before bind_eos_to_source (so they win and
        ! their targets stay private). Note: this makes usr_Rfactor apply only to
        ! FI for HD -- consistent with MHD, where the EoS owns R for LTE/PI.
        if(associated(usr_rfactor)) then
            eos%get_Rfactor=>usr_rfactor
        else
            eos%get_Rfactor=>rfactor_from_constant_ionization
        end if
 
        !> PI energy-mode overrides (eos_type='PI', ionE=.true.)
        ! Energy mode differs from no-energy PI only in the eint<->p relation,
        ! so override exactly the routines that touch it (conversions, p_to_e,
        ! csound2/gamma1, Te update), all via the portable backend. pthermal
        ! is unchanged: hd_get_pthermal already carries a PI-energy branch.
        if (eos%eos_type == 'PI' .and. eos%ionE) then
            if (.not. hd_energy) &
                call mpistop('PI energy EoS requires hd_energy=.true.')
            eos%to_conserved  => hd_to_conserved_pi
            eos%to_primitive  => hd_to_primitive_pi
            phys_to_primitive => eos%to_primitive
            phys_to_conserved => eos%to_conserved
            eos%p_to_e        => p_to_e_pi
            !> eos%get_csound2 => get_csound2_PI set in eos_finalise_PI (private target)
            phys_get_gamma1   => get_gamma1_pi
            ! Te_ from eint each substep (uses HD's Te_ index: PI registers
            ! Te via var_set_auxvar, so the generic iw_te is unset).
        end if
        hd_to_primitive          => eos%to_primitive
        hd_to_conserved          => eos%to_conserved
 
    end subroutine hd_link_eos
 
    subroutine bind_eos_to_source() !> this is called in eos_finalise through mod_physics procedure linking
        if (allocated(tc_fl)) then
            tc_fl%get_temperature_from_conserved => eos%get_temperature_from_etot
            !> Use fast bilinear T lookup for TC STS substeps (density fixed,
            !> TC flux dominated by corona where T(eint) is smooth)
            if (eos%eos_type == 'LTE' .and. eos%ionE) then
                tc_fl%get_temperature_from_eint => get_temperature_from_eint_fast_lte
            else
                tc_fl%get_temperature_from_eint => eos%get_temperature_from_eint
            end if
            tc_fl%get_rho => eos%get_rho
            tc_fl%get_ne_nH => eos%get_ne_nH
            tc_fl%get_var_Rfactor => eos%get_Rfactor
            tc_fl%inv_gamma_minus_1 =  eos%inv_gamma_minus_1
            tc_fl%nH2rhoFactor      =  eos%nH2rhoFactor
            tc_fl%log_T_floor       =  eos_get_log_t_floor()
            tc_fl%eint_from_T       => eint_nh_from_t
        end if
 
        if (allocated(rc_fl)) then
            rc_fl%get_rho => eos%get_rho
            rc_fl%get_pthermal => eos%get_thermal_pressure
            rc_fl%get_var_Rfactor => eos%get_Rfactor
            rc_fl%get_Te => eos%get_Te
            rc_fl%get_ne_nH => eos%get_ne_nH
            rc_fl%ionE              =  eos%ionE
            rc_fl%method            =  eos%method
            rc_fl%inv_gamma_minus_1 =  eos%inv_gamma_minus_1
            rc_fl%nH2rhoFactor      =  eos%nH2rhoFactor
            rc_fl%eion_per_nH       =  eos%eion_per_nH
            rc_fl%eint_from_T       => eint_nh_from_t
            rc_fl%p2eint            => p2eint_from_nh_p
            rc_fl%T_from_eint       => t_from_nh_eint
            rc_fl%y_from_eint       => y_from_nh_eint
            !> Build the variable-c_V Townsend Y_mod table now that all
            !> EoS tables (eint_from_T, T, neOnH) are in code units.
            !> build_Y_mod_table checks coolmethod=='exact' and .not.isPPL
            !> internally and early-returns otherwise. LTE only:
            !> eos_get_eintT_grid reads the LTE eintT grid; PI runs classical
            !> Townsend (Y_mod-for-PI is a later refinement).
            if (eos%ionE .and. eos%eos_type == 'LTE') then
                call eos_get_eintt_grid(rc_fl%Y_mod_n_nH, &
                     rc_fl%Y_mod_lg_nH_min, rc_fl%Y_mod_lg_nH_max)
                call build_y_mod_table(rc_fl)
            end if
        end if
 
        !> PI energy mode: repoint cooling/conduction eint<->T callbacks at
        !> the ionisation backend (the LTE-table versions above assume LTE
        !> tables PI does not load). Classical Townsend; Y_mod not built (above).
        if (eos%eos_type == 'PI' .and. eos%ionE) then
            if (allocated(tc_fl)) tc_fl%eint_from_T => eint_from_t_pi
            if (allocated(rc_fl)) then
                rc_fl%eint_from_T => eint_from_t_pi
                rc_fl%p2eint      => p2eint_pi
                rc_fl%T_from_eint => t_from_eint_pi
                rc_fl%y_from_eint => y_from_eint_pi
            end if
        end if
 
        if (allocated(fld_fl)) then
            !> Radiation (FLD) fluid: gas-EoS callbacks. get_temperature_from_pressure
            !> is the FI T=p/(R*rho) routine (the LTE variant is a later pass).
            fld_fl%gamma       =  eos%gamma
            fld_fl%get_tgas    => eos%get_temperature_from_pressure
            fld_fl%get_Rfactor => eos%get_Rfactor
        end if
    end subroutine bind_eos_to_source
 
    !> Transform primitive variables into conservative ones
    subroutine hd_to_conserved_origin(ixI^L, ixO^L, w, x)
        use mod_global_parameters
        use mod_dust, only: dust_to_conserved
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        integer :: ix^d
 
        timeeos0 = mpi_wtime() !> For monitoring cost of eos module
 
        {do ix^db=ixomin^db,ixomax^db\}
            if (hd_energy) then
                ! Calculate total energy from pressure and kinetic energy
                w(ix^d,e_)=w(ix^d,p_)*eos%inv_gamma_minus_1+&
                    half*(^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)
            end if
            ! Convert velocity to momentum
            ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)\
        {end do\}
 
        if (hd_dust) then
            call dust_to_conserved(ixi^l, ixo^l, w, x)
        end if
 
        timeeos_conv=timeeos_conv+(mpi_wtime()-timeeos0)
 
    end subroutine hd_to_conserved_origin
 
    !> Transform conservative variables into primitive ones
    subroutine hd_to_primitive_origin(ixI^L, ixO^L, w, x)
        use mod_global_parameters
        use mod_dust, only: dust_to_primitive
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision                :: inv_rho
        integer :: ix^d
 
        timeeos0 = mpi_wtime() !> For monitoring cost of eos module
 
        if (fix_small_values) then
            call hd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'hd_to_primitive')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
            ! Calculate pressure = (gamma-1) * (e-ek)
            if(hd_energy) then
                ! Compute pressure
                w(ix^d,p_)=(eos%gamma_minus_1)*(w(ix^d,e_)&
                    -half*w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+))
            end if
        {end do\}
 
        ! Convert dust momentum to dust velocity
        if (hd_dust) then
            call dust_to_primitive(ixi^l, ixo^l, w, x)
        end if
 
        timeeos_conv=timeeos_conv+(mpi_wtime()-timeeos0)
 
    end subroutine hd_to_primitive_origin
 
    !> LTE primitive -> conserved conversion.
    !>
    !> On entry: rho_ = density, m_ = velocity, p_ = pressure.
    !> On exit:  rho_ = density (unchanged), m_ = momentum, e_ = total energy.
    subroutine hd_to_conserved_lte(ixI^L, ixO^L, w, x)
        use mod_global_parameters
        use mod_dust, only: dust_to_conserved
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        timeeos0 = mpi_wtime()
 
        call p_to_e(ixi^l, ixo^l, w, x)
 
        ! Convert velocity to momentum
        ^c&w(ixo^s,m^c_)=w(ixo^s,rho_)*w(ixo^s,m^c_)\
 
        if (hd_dust) then
            call dust_to_conserved(ixi^l, ixo^l, w, x)
        end if
 
        timeeos_conv=timeeos_conv+(mpi_wtime()-timeeos0)
 
    end subroutine hd_to_conserved_lte
 
    subroutine p_to_e(ixI^L, ixO^L, w, x)
        !> Convert pressure to total energy: E = eint(rho, p) + KE.
        !>
        !> On entry: w(rho_) = density, w(m_) = velocity, w(p_) = pressure.
        !> On exit:  w(rho_) unchanged, w(m_) unchanged, w(e_) = total energy.
        !>
        !> Four paths for ionE (LTE with ionisation):
        !>   1. FI bypass (p/rho > threshold): analytic eint = p/(gamma-1) + eion*nH
        !>   2. Analytical Saha (eos_method='analytic'): direct Saha solve for T,y from p
        !>   3. WB mode (hd_well_balanced): cached bisection on forward tables
        !>   4. Standard: p2eint inverse table lookup (fast, ~0.01% round-trip error)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        integer :: ix^d
        double precision :: p_to_eint, p_over_rho
        double precision :: nh(ixi^s), nh_in(ixi^s), p_in(ixi^s)
        double precision :: log_eint_mid, eint_total
        double precision :: t_solve, y_solve, eint_nh_solve
 
        if (eos%ionE) then
            call eos%get_nH(w, x, ixi^l, ixo^l, nh)
            nh_in(ixo^s) = dlog10(nh(ixo^s))
            if (eos%p2eint_method /= 'bisect') then
                p_in(ixo^s) = dlog10(w(ixo^s,p_)) - nh_in(ixo^s)
            end if
        endif
 
        p_to_eint = eos%inv_gamma_minus_1
        {do ix^db=ixomin^db,ixomax^db\}
            if (hd_energy) then
                if (eos%ionE) then
                    p_over_rho = w(ix^d,p_) / w(ix^d,rho_)
                    if (p_over_rho > eos%p_rho_FI_threshold) then
                        !> FI bypass: exact inverse of to_primitive
                        p_to_eint = eos%inv_gamma_minus_1 &
                            + eos%eion_per_nH * nh(ix^d) / w(ix^d,p_)
                        w(ix^d,e_) = w(ix^d,p_)*p_to_eint + &
                            half*(^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)
                    else if (eos%method == 'analytic') then
                        !> Analytical Saha: bisect for T from p, return eint directly
                        call saha_state_from_nh_p(nh(ix^d), w(ix^d,p_), &
                            t_solve, y_solve, eint_nh_solve)
                        w(ix^d,e_) = eint_nh_solve * nh(ix^d) + &
                            half*(^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)
                    else if (eos%p2eint_method == 'bisect') then
                        !> Bisection on forward p table: finds eint such that
                        !> p(nH, eint) = p_target exactly. More accurate than
                        !> the p2eint table in the ionisation zone.
                        !> Entropy mode uses the pfwd/eintP pair: eos%log_p (used by
                        !> eint_from_p_bisect) is not loaded there.
                        if (eos%method == 'entropy') then
                            call entropy_eint_from_p_bisect( &
                                eos%pfwd, eos%pfwd_x, eos%pfwd_y, eos%pfwd_xy, &
                                eos%eintP, eos%eintP_x, eos%eintP_y, eos%eintP_xy, &
                                nh_in(ix^d), dlog10(w(ix^d,p_)) - nh_in(ix^d), &
                                log_eint_mid)
                            eint_total = nh(ix^d) * 10.0d0**log_eint_mid
                        else
                            call eint_from_p_bisect(nh_in(ix^d), &
                                dlog10(w(ix^d,p_)), log_eint_mid)
                            eint_total = nh(ix^d) * 10.0d0**log_eint_mid
                            eint_total = max(eint_total, &
                                nh(ix^d) * 10.0d0**eos%T%var2_min)
                        end if
                        w(ix^d,e_) = eint_total + &
                            half*(^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)
                    else
                        !> Standard table lookup (fast, default)
                        p_to_eint = p2eint_from_nh_p(nh_in(ix^d), p_in(ix^d))
                        w(ix^d,e_) = w(ix^d,p_)*p_to_eint + &
                            half*(^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)
                    end if
                else
                    w(ix^d,e_) = w(ix^d,p_)*p_to_eint + &
                        half*(^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)
                end if
            end if
        {end do\}
 
    end subroutine p_to_e
 
    !> LTE conserved -> primitive conversion.
    !>
    !> On entry: rho_ = density, m_ = momentum, e_ = total energy.
    !> On exit:  rho_ = density (unchanged), m_ = velocity, p_ = pressure.
    !>
    !> Pressure is computed energy-consistently from actual eint via EoS
    !> table lookups (T and ne/nH). Cannot use stored Ne_/Te_ because
    !> they may be stale after AMR prolongation/coarsening.
    subroutine hd_to_primitive_lte(ixI^L, ixO^L, w, x)
        use mod_global_parameters
        use mod_dust, only: dust_to_primitive
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision                :: inv_rho
        double precision                :: nh(ixi^s)
        double precision                :: log_nh(ixi^s)
        double precision                :: eint_val, eint_in, t_loc, y_loc
        integer :: ix^d
 
        timeeos0 = mpi_wtime()
 
        if (fix_small_values) then
            call hd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'hd_to_primitive_LTE')
        end if
 
        call eos%get_nH(w, x, ixi^l, ixo^l, nh)
 
        ! Cache log10(nH) for all cells (used by table lookups for IonE)
        if (eos%ionE) then
            log_nh(ixo^s) = dlog10(nh(ixo^s))
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
            ! Calculate pressure
            if(hd_energy) then
                if (eos%ionE) then
                    ! Energy-consistent pressure from actual eint via EoS tables.
                    ! Cannot use stored Ne_/Te_ because they may be stale after
                    ! AMR prolongation/coarsening (nonlinear EoS breaks averaging).
                    eint_val = w(ix^d,e_) - half*w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+)
                    ! Floor eint to prevent unphysical values
                    if (eos%method /= 'analytic') then
                        eint_val = max(eint_val, nh(ix^d) * 10.0d0**eos%T%var2_min)
                    end if
                    eint_val = max(eint_val, smalldouble)
                    if (eint_val * inv_rho > eos%eint_rho_FI_threshold) then
                        ! FI bypass: p = (gamma-1)*(eint - eion*nH)
                        w(ix^d,p_) = eos%gamma_minus_1 &
                            * (eint_val - eos%eion_per_nH * nh(ix^d))
                    else
                        ! Ionisation zone: single p/nH lookup
                        eint_in = dlog10(eint_val) - log_nh(ix^d)
                        w(ix^d,p_) = nh(ix^d) &
                            * p_nh_from_eint(log_nh(ix^d), eint_in)
                    end if
                else
                    w(ix^d,p_)=(eos%gamma_minus_1)*(w(ix^d,e_)&
                        -half*w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+))
                end if
            end if
        {end do\}
 
        ! Convert dust momentum to dust velocity
        if (hd_dust) then
            call dust_to_primitive(ixi^l, ixo^l, w, x)
        end if
 
        timeeos_conv=timeeos_conv+(mpi_wtime()-timeeos0)
 
    end subroutine hd_to_primitive_lte
 
    !> Calculate thermal pressure within ixO^L.
    !> For energy runs delegates to eos%get_thermal_pressure; for no-energy
    !> (isothermal) uses the adiabatic relation p = adiab * rho^gamma.
    subroutine hd_get_pthermal(w, x, ixI^L, ixO^L, pth)
        use mod_global_parameters
        use mod_physics, only: phys_get_ei
        use mod_usr_methods, only: usr_set_pthermal
        use mod_small_values, only: trace_small_values
 
        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):: pth(ixi^s)
        integer                      :: iw, ix^d
        double precision :: nh(ixi^s), ei(ixi^s), tpi, rpi
 
        if (hd_energy) then
            if (eos%eos_type == 'LTE') then
                ! LTE: p = nH * (1 + He + ne/nH) * T from stored state
                call eos%get_nH(w, x, ixi^l, ixo^l, nh)
                pth(ixo^s) = nh(ixo^s) * (1.0d0 + eos%He_abundance &
                    + (w(ixo^s,iw_ne) / nh(ixo^s))) * w(ixo^s,iw_te)
            else if (eos%eos_type == 'PI' .and. eos%ionE) then
                ! PI energy: invert eint (= e - KE, via phys_get_ei) -> p with
                ! the ionisation backend (eint carries the ionisation energy).
                ei(ixo^s) = phys_get_ei(w, ixi^l, ixo^l)
                {do ix^db=ixomin^db,ixomax^db\}
                    call state_from_eint_pi(w(ix^d,rho_), ei(ix^d), &
                        tpi, pth(ix^d), rpi)
                {end do\}
            else
                ! FI: p = (gamma-1) * (e - KE)
                pth(ixo^s) = eos%gamma_minus_1 * phys_get_ei(w, ixi^l, ixo^l)
            end if
        else
            if (.not. associated(usr_set_pthermal)) then
                pth(ixo^s) = hd_adiab * w(ixo^s, rho_)**eos%gamma
            else
                call usr_set_pthermal(w,x,ixi^l,ixo^l,pth)
            end if
        end if
 
        if (fix_small_values) then
            {do ix^db= ixo^lim^db\}
                if(pth(ix^d)<small_pressure) then
                    pth(ix^d)=small_pressure
                endif
            {enddo^d&\}
        else if (check_small_values) then
            {do ix^db= ixo^lim^db\}
                if(pth(ix^d)<small_pressure) then
                    write(*,*) "Error: small value of gas pressure",pth(ix^d),&
                        " encountered when call hd_get_pthermal"
                    write(*,*) "Iteration: ", it, " Time: ", global_time
                    write(*,*) "Location: ", x(ix^d,:)
                    write(*,*) "Cell number: ", ix^d
                    do iw=1,nw
                        write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
                    end do
                    if(trace_small_values) write(*,*) dsqrt(pth(ix^d)-bigdouble)
                    write(*,*) "Saving status at the previous time step"
                    crash=.true.
                end if
            {enddo^d&\}
        end if
 
    end subroutine hd_get_pthermal
 
    !> Sound speed squared for FI (fully ionized / constant gamma) EoS.
    !> Expects w in primitive form: w(p_) = pressure, w(rho_) = density.
    subroutine hd_get_csound2_fi(w, x, ixI^L, ixO^L, cs2)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(in)    :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
        double precision, intent(out)   :: cs2(ixi^s)
 
        timeeos0 = mpi_wtime()
 
        cs2(ixo^s) = eos%gamma * w(ixo^s, p_) / w(ixo^s, rho_)
 
        timeeos_csound = timeeos_csound + (mpi_wtime()-timeeos0)
 
    end subroutine hd_get_csound2_fi
 
    !> Sound speed squared for LTE+IonE EoS.
    !> Delegates Gamma_1 computation to hd_get_gamma1_LTE, then cs2 = Gamma_1 * p/rho.
    subroutine hd_get_csound2_lte(w, x, ixI^L, ixO^L, cs2)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(in)    :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
        double precision, intent(out)   :: cs2(ixi^s)
        integer :: ix^d
 
        timeeos0 = mpi_wtime()
 
        call hd_get_gamma1_lte(w, x, ixi^l, ixo^l, cs2)
        {do ix^db=ixomin^db,ixomax^db\}
            cs2(ix^d) = cs2(ix^d) * w(ix^d, p_) / w(ix^d, rho_)
        {end do\}
 
        timeeos_csound = timeeos_csound + (mpi_wtime()-timeeos0)
 
    end subroutine hd_get_csound2_lte
 
    !> PI energy-mode routines (eos_type='PI', ionE=.true.)
    !> Mirror the HD FI/LTE family; the eint<->p relation is delegated to the
    !> portable scalar backend (mod_eos_PI). Same eq_state_units /
    !> RR=1 normalisation as FI. Origin (total-energy) only, as in HD. eint
    !> carries the ionisation-energy term, so p=(gamma-1)*eint no longer holds.
 
    !> Primitive pressure -> total energy via backend (m_ is velocity here).
    subroutine p_to_e_pi(ixI^L, ixO^L, w, x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
        double precision :: eint
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            if (hd_energy) then
                call eint_from_rho_p_pi(w(ix^d,rho_), w(ix^d,p_), eint)
                w(ix^d,e_) = eint + half*(^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)
            end if
        {end do\}
    end subroutine p_to_e_pi
 
    !> Primitive -> conserved (PI energy)
    subroutine hd_to_conserved_pi(ixI^L, ixO^L, w, x)
        use mod_global_parameters
        use mod_dust, only: dust_to_conserved
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        timeeos0 = mpi_wtime()
        call p_to_e_pi(ixi^l, ixo^l, w, x)
        ! Convert velocity to momentum
        ^c&w(ixo^s,m^c_)=w(ixo^s,rho_)*w(ixo^s,m^c_)\
        if (hd_dust) call dust_to_conserved(ixi^l, ixo^l, w, x)
        timeeos_conv=timeeos_conv+(mpi_wtime()-timeeos0)
    end subroutine hd_to_conserved_pi
 
    !> Conserved -> primitive (PI energy): KE removed -> eint, backend -> p.
    subroutine hd_to_primitive_pi(ixI^L, ixO^L, w, x)
        use mod_global_parameters
        use mod_dust, only: dust_to_primitive
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
        double precision :: inv_rho, eint_val, t, rfac
        integer :: ix^d
 
        timeeos0 = mpi_wtime()
 
        if (fix_small_values) then
            call hd_handle_small_values(.false., w, x, ixi^l, ixo^l, &
                'hd_to_primitive_PI')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
            if (hd_energy) then
                eint_val = w(ix^d,e_) &
                    - half*w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+)
                eint_val = max(eint_val, smalldouble)
                call state_from_eint_pi(w(ix^d,rho_), eint_val, &
                    t, w(ix^d,p_), rfac)
            end if
        {end do\}
 
        if (hd_dust) call dust_to_primitive(ixi^l, ixo^l, w, x)
        timeeos_conv=timeeos_conv+(mpi_wtime()-timeeos0)
    end subroutine hd_to_primitive_pi
 
    !> Adiabatic sound speed squared from primitive (rho, p).
 
    !> Effective Gamma1 = cs2 * rho / p.
 
    !> PI energy mode: refresh Te_ from gas internal energy via the backend
    !> eint->T inversion. Uses HD's Te_ index (PI registers Te via
    !> var_set_auxvar, so the generic iw_te is unset). phys_get_ei supplies
    !> eint (KE removed).
 
    !> PI no-energy R-factor from the stored Te_ (HD's index). Mirrors MHD's
    !> Rfactor_from_temperature_ionization. The generic mod_eos version uses
    !> the global iw_te, which PI leaves unset (var_set_auxvar) -- so HD needs
    !> its own Te_-addressed routine, exactly as MHD does.
 
    !> PI no-energy Te_ update (pth/(rho*R) with lagged iz_H from wCT(Te_)).
    !> HD's Te_-addressed mirror of MHD's mhd_update_temperature.
 
    !> Prominence (T,p) R-factor: recompute from (rho, pth) via the module's
    !> prominence inversion (no-energy only -> pth=(gamma-1)*eint, R-indep).
 
    !> Prominence (T,p) Te_ update from (rho, pth).
 
    !> Return effective adiabatic index for LTE+IonE EoS.
    !> Dispatches on gamma1_method and eos%method.
    subroutine hd_get_gamma1_lte(w, x, ixI^L, ixO^L, gamma1)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(in)    :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
        double precision, intent(out)   :: gamma1(ixi^s)
 
        double precision :: nh_val, p_over_rho
        integer :: ix^d
 
        if (eos%gamma1_method == 'constant') then
            gamma1(ixo^s) = eos%gamma
            return
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            p_over_rho = w(ix^d, p_) / w(ix^d, rho_)
            if (p_over_rho > eos%p_rho_FI_threshold) then
                gamma1(ix^d) = eos%gamma
            else
                nh_val = w(ix^d, rho_) / eos%nH2rhoFactor
                if (eos%method == 'analytic') then
                    if (iw_te > 0 .and. w(ix^d,iw_te) > 0.0d0) then
                        gamma1(ix^d) = saha_gamma1_from_nh_t(nh_val, w(ix^d,iw_te))
                    else
                        gamma1(ix^d) = eos%gamma
                    end if
                else
                    gamma1(ix^d) = gamma1_from_nh_p(dlog10(nh_val), &
                        dlog10(w(ix^d, p_) / nh_val))
                end if
            end if
        {end do\}
 
    end subroutine hd_get_gamma1_lte
 
    !> Rfactor = p/(rho*T) for constant ionisation degree (FI/PI no-energy).
    !> Stays in the seam: RR is the physics module's gas-constant factor, not
    !> visible to mod_eos. Rfactor_from_LTE lives in mod_eos (no RR).
    subroutine rfactor_from_constant_ionization(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)
 
        rfactor(ixo^s)=rr
 
    end subroutine rfactor_from_constant_ionization
 
    !> Convert conserved (rho, rho*v, E) to prolong form (rho, v, T).
    !> T is stored in the e_ slot. Interpolation in this space avoids
    !> Jensen's inequality across the ionisation plateau.
    subroutine hd_to_prolong_lte(ixI^L, ixO^L, w, x)
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: inv_rho, eint_val, nh_val, log_nh, t_loc, y_loc
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0 / w(ix^d, rho_)
            nh_val = w(ix^d, rho_) / eos%nH2rhoFactor
            log_nh = dlog10(nh_val)
 
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
 
            ! Compute eint = E - 0.5*rho*v^2
            eint_val = w(ix^d, e_) &
                - half * w(ix^d, rho_) * (^c&w(ix^d,m^c_)**2+)
            ! Floor eint to prevent unphysical values
            if (eos%method /= 'analytic') then
                eint_val = max(eint_val, nh_val * 10.0d0**eos%T%var2_min)
            end if
            eint_val = max(eint_val, smalldouble)
 
            ! Convert eint to T via EoS
            if (eint_val * inv_rho > eos%eint_rho_FI_threshold) then
                ! FI bypass: T = (gamma-1)*(eint - eion*nH) / (nH * n_per_nH_FI)
                w(ix^d, e_) = eos%gamma_minus_1 &
                    * (eint_val - eos%eion_per_nH * nh_val) &
                    / (nh_val * eos%n_per_nH_FI)
            else if (eos%method == 'analytic') then
                ! Analytical Saha: solve for T from eint
                call saha_t_from_nh_eint(nh_val, &
                    eint_val / nh_val, t_loc, y_loc)
                w(ix^d, e_) = t_loc
            else
                ! Ionisation zone: T from table
                w(ix^d, e_) = t_from_nh_eint(log_nh, &
                    dlog10(eint_val) - log_nh)
            end if
        {end do\}
 
    end subroutine hd_to_prolong_lte
 
    !> Convert prolong form (rho, v, T) back to conserved (rho, rho*v, E).
    !> T is read from the e_ slot. Uses eint_nH_from_T for back-conversion.
    subroutine hd_from_prolong_lte(ixI^L, ixO^L, w, x)
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: t_val, eint_val, t_fi, nh_val, log_nh, log_t_min
        integer :: ix^d
 
        t_fi = (eos%eint_rho_FI_threshold &
            * eos%nH2rhoFactor - eos%eion_per_nH) &
            * eos%gamma_minus_1 / eos%n_per_nH_FI
 
        ! Floor for log_T when calling the (rho, T) inverse table.
        ! Legacy 'tables' method populates eos%eint_from_T; entropy method
        ! populates eos%eintT. Picking the wrong container leaves var2_min = 0
        ! (uninitialised), which floors T at 10^0 = 1 code unit (= 10^6 K) —
        ! that clobbers any cold cell going through AMR prolongation.
        if (eos%method == 'entropy') then
            log_t_min = eos%eintT%var2_min
        else
            log_t_min = eos%eint_from_T%var2_min
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            t_val = w(ix^d, e_)  ! T stored in e_ slot
            nh_val = w(ix^d, rho_) / eos%nH2rhoFactor
            log_nh = dlog10(nh_val)
 
            if (t_val > t_fi) then
                ! FI: eint = nH*(n_per_nH*T/(gamma-1) + eion)
                eint_val = nh_val &
                    * (eos%n_per_nH_FI * t_val * eos%inv_gamma_minus_1 &
                    + eos%eion_per_nH)
            else if (eos%method == 'analytic') then
                ! Analytical Saha: eint from T directly
                eint_val = saha_eint_from_nh_t(nh_val, t_val) * nh_val
            else
                ! Ionisation zone: eint/nH from T table
                eint_val = eint_nh_from_t(log_nh, &
                    dlog10(max(t_val, 10.0d0**log_t_min))) &
                    * nh_val
            end if
 
            ! E = eint + 0.5*rho*v^2
            w(ix^d, e_) = eint_val &
                + half * w(ix^d, rho_) * (^c&w(ix^d,m^c_)**2+)
            ! Convert velocity to momentum
            ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)\
        {end do\}
 
    end subroutine hd_from_prolong_lte
 
end module mod_hd_eos
!> Needs a line after to pass the preprocesor