mod_ffhd_eos.t

Go to the documentation of this file.

!=============================================================================
!> FFHD <-> EoS seam: binds the eos% authority into force-free hydrodynamics.
!>
!> Mirrors mod_hd_eos.t, adapted to FFHD's structure: a single field-aligned
!> momentum component mom(1) and NO magnetic energy (force-free, frozen B0), so
!> every kinetic-energy term is half*rho*mom(1)**2. The B0 projection lives in
!> the flux/cmax/velocity layer (mod_ffhd_phys.t) and is orthogonal to the EoS
!> seam, so the thermodynamic conversions are HD-like with one momentum
!> component.
!>
!> ffhd_link_eos        -- wires eos%/phys_ pointers per eos_type (FI now;
!>                         LTE/PI guarded until milestones 2/3).
!> ffhd_bind_eos_to_source -- wires the thermal-conduction (tc_fl), radiative-
!>                         cooling (rc_fl) and thermal-emission (te_fl_ffhd)
!>                         fluid-port callbacks from eos%. This is what the
!>                         pre-eos% hand-wiring was missing (rc_fl%get_Te /
!>                         get_ne_nH), which caused the cooling SIGSEGV.
!=============================================================================
module mod_ffhd_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_saha
    use mod_eos_pi
    use mod_eos_container
    use mod_ffhd_phys
    use mod_radiative_cooling, only: build_y_mod_table
 
    use mod_comm_lib, only: mpistop
 
    implicit none
    private
 
    public :: ffhd_link_eos
 
contains
 
    !> Link the appropriate EoS conversion/closure routines per eos_type.
    !> Called from ffhd_activate after ffhd_phys_init (mirrors hd_activate).
    subroutine ffhd_link_eos()
        use mod_usr_methods
 
        ! Isothermal FFHD (no energy equation) has no thermodynamic state for
        ! the EoS to close -- only FI is meaningful there.
        if (.not. ffhd_energy .and. eos%eos_type /= 'FI') &
            call mpistop("FFHD: ffhd_energy=.false. requires eos_type='FI'")
 
        ! FI / PI base conversions = ffhd's origin routines (single-component
        ! KE). PI energy mode overrides them below (milestone 3). LTE routes
        ! primitive<->conserved through ffhd_p_to_e / table inversion.
        if (eos%eos_type == 'FI' .or. eos%eos_type == 'PI') then
            eos%to_conserved => ffhd_to_conserved_origin
            eos%to_primitive => ffhd_to_primitive_origin
        else if (eos%eos_type == 'LTE') then
            eos%to_conserved => ffhd_to_conserved_lte
            eos%to_primitive => ffhd_to_primitive_lte
        else
            call mpistop("Unknown FFHD EOS type: " // trim(eos%eos_type))
        end if
 
        phys_to_primitive => eos%to_primitive
        phys_to_conserved => eos%to_conserved
        ffhd_to_primitive => eos%to_primitive
        ffhd_to_conserved => eos%to_conserved
        phys_get_rho      => eos%get_rho
 
        ! suitable for both FI (origin inlines it) and LTE (conversions call it)
        eos%p_to_e => ffhd_p_to_e
 
        ! pthermal: the isothermal (no-energy) pointer is wired in phys_init and
        ! left untouched; for the energy equation route through eos%. LTE reads
        ! the stored Te_/Ne_ aux state; FI/PI use (gamma-1)*(e-KE).
        if (ffhd_energy) then
            if (eos%eos_type == 'LTE') then
                phys_get_pthermal        => ffhd_get_pthermal_lte
                ffhd_get_pthermal        => ffhd_get_pthermal_lte
                eos%get_thermal_pressure => ffhd_get_pthermal_lte
            else
                phys_get_pthermal        => ffhd_get_pthermal_origin
                ffhd_get_pthermal        => ffhd_get_pthermal_origin
                eos%get_thermal_pressure => ffhd_get_pthermal_origin
            end if
        end if
 
        phys_bind_eos_to_source => ffhd_bind_eos_to_source
 
        ! Energy<->internal heads: ffhd previously handed these only to the TC
        ! STS via set_conversion_methods_to_head; the LTE fast-T kernel and
        ! phys_update_temperature also need the module-level phys_ pointers.
        phys_e_to_ei => ffhd_e_to_ei
        phys_ei_to_e => ffhd_ei_to_e
        phys_get_ei  => ffhd_get_ei   ! LTE+ionE cooling recovers eint via this
 
        ! sound speed / gamma1. LTE+ionE: Gamma_1(nH,p) from the EoS tables and
        ! EoS-aware AMR prolongation (interpolate in (rho,v,T) to avoid Jensen's
        ! inequality across the ionisation plateau). Otherwise FI: cs2=gamma*p/rho.
        if (eos%eos_type == 'LTE' .and. eos%ionE) then
            eos%get_csound2   => ffhd_get_csound2_lte
            phys_get_gamma1   => ffhd_get_gamma1_lte
            phys_to_prolong   => ffhd_to_prolong_lte
            phys_from_prolong => ffhd_from_prolong_lte
        else
            eos%get_csound2 => ffhd_get_csound2
            phys_get_gamma1 => get_gamma1_fi
        end if
 
        ! ffhd_get_temperature (phys pointer) by eos_type: PI reads the cached
        ! Te_, the others recompute from etot.
        if (eos%eos_type == 'PI') then
            ffhd_get_temperature => ffhd_get_temperature_from_te
        else
            ffhd_get_temperature => ffhd_get_temperature_from_etot
        end if
        ! Rfactor: only the FI case is physics-dependent (usr_Rfactor). LTE and PI
        ! bind pure eos% routines in eos_finalise_{LTE,PI} (run after this, before
        ! ffhd_bind_eos_to_source -> targets stay private). usr_Rfactor thus applies
        ! to FI only, consistent with mhd/hd.
        if (associated(usr_rfactor)) then
            eos%get_Rfactor => usr_rfactor
        else
            eos%get_Rfactor => rfactor_from_constant_ionization
        end if
        ! ffhd_get_Rfactor (runtime alias) captured in ffhd_bind_eos_to_source,
        ! AFTER eos_finalise has set eos%get_Rfactor for LTE/PI.
 
        !> PI energy-mode overrides (eos_type='PI', ionE=.true.). Energy mode
        !> differs from no-energy PI only in the eint<->p relation (eint carries
        !> the ionisation energy), so override exactly the routines touching it,
        !> all via the portable mod_eos_PI backend. get_Rfactor stays the
        !> Te_-addressed routine above (Te_ refreshed each step by
        !> ffhd_update_temperature_PI).
        if (eos%eos_type == 'PI' .and. eos%ionE) then
            if (.not. ffhd_energy) &
                call mpistop('FFHD PI energy EoS requires ffhd_energy=.true.')
            eos%to_conserved  => ffhd_to_conserved_pi
            eos%to_primitive  => ffhd_to_primitive_pi
            phys_to_primitive => eos%to_primitive
            phys_to_conserved => eos%to_conserved
            ffhd_to_primitive => eos%to_primitive
            ffhd_to_conserved => eos%to_conserved
            eos%p_to_e        => ffhd_p_to_e_pi
            phys_get_pthermal        => ffhd_get_pthermal_pi
            ffhd_get_pthermal        => ffhd_get_pthermal_pi
            eos%get_thermal_pressure => ffhd_get_pthermal_pi
            !> eos%get_csound2 => get_csound2_PI set in eos_finalise_PI (private target)
            phys_get_gamma1   => get_gamma1_pi
        end if
 
    end subroutine ffhd_link_eos
 
    !> Wire the source-term fluid ports from eos%. Called from eos_finalise via
    !> phys_bind_eos_to_source, AFTER eos_finalise has set eos%get_Te,
    !> get_ne_nH, eion_per_nH, n_per_nH_FI/neOnH_FI, get_temperature_from_*.
    subroutine ffhd_bind_eos_to_source()
 
        !> Runtime R-factor alias: eos%get_Rfactor is fully set by now
        !> (eos_finalise ran before phys_bind_eos_to_source), incl. the LTE/PI
        !> bindings done in eos_finalise_{LTE,PI}.
        ffhd_get_rfactor => eos%get_Rfactor
 
        if (allocated(tc_fl)) then
            tc_fl%get_temperature_from_conserved => eos%get_temperature_from_etot
            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
            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 eint<->p kernels at the portable backend.
        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
 
        {^ifthreed
        if (allocated(te_fl_ffhd)) then
            te_fl_ffhd%get_rho         => eos%get_rho
            te_fl_ffhd%get_pthermal    => eos%get_thermal_pressure
            te_fl_ffhd%get_var_Rfactor => eos%get_Rfactor
            te_fl_ffhd%get_ne_nH       => eos%get_ne_nH
        end if
        }
 
    end subroutine ffhd_bind_eos_to_source
 
    !=========================================================================
    !> LTE (Saha table) thermodynamics for FFHD. Mirrors mod_hd_eos.t adapted
    !> to FFHD's single field-aligned momentum mom(1) and total-energy-only
    !> ("origin") formulation; no dust, no magnetic energy.
    !=========================================================================
 
    !> LTE primitive -> conserved. On entry rho/v/p; on exit rho/mom/E.
    subroutine ffhd_to_conserved_lte(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)
 
        call ffhd_p_to_e(ixi^l, ixo^l, w, x)
        w(ixo^s,mom(1)) = w(ixo^s,rho_)*w(ixo^s,mom(1))
    end subroutine ffhd_to_conserved_lte
 
    !> Convert pressure to total energy: E = eint(rho, p) + KE.
    !> On entry rho/v/p; on exit w(e_) = total energy (rho, mom unchanged).
    !> Four ionE paths mirror hd p_to_e: FI bypass, analytic Saha, p-bisect,
    !> p2eint table. Non-ionE LTE uses the constant inverse-gamma relation.
    subroutine ffhd_p_to_e(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)
 
        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 (.not. ffhd_energy) return
 
        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 (eos%ionE) then
                p_over_rho = w(ix^d,p_) / w(ix^d,rho_)
                if (p_over_rho > eos%p_rho_FI_threshold) then
                    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*w(ix^d,mom(1))**2*w(ix^d,rho_)
                else if (eos%method == 'analytic') then
                    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*w(ix^d,mom(1))**2*w(ix^d,rho_)
                else if (eos%p2eint_method == 'bisect') then
                    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)
                    w(ix^d,e_) = eint_total + &
                        half*w(ix^d,mom(1))**2*w(ix^d,rho_)
                else
                    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*w(ix^d,mom(1))**2*w(ix^d,rho_)
                end if
            else
                w(ix^d,e_) = w(ix^d,p_)*p_to_eint + &
                    half*w(ix^d,mom(1))**2*w(ix^d,rho_)
            end if
        {end do\}
    end subroutine ffhd_p_to_e
 
    !> LTE conserved -> primitive. On entry rho/mom/E; on exit rho/v/p.
    !> Pressure is energy-consistent from actual eint (stored Ne_/Te_ may be
    !> stale after AMR prolong/coarsen on the nonlinear EoS).
    subroutine ffhd_to_primitive_lte(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                :: inv_rho
        double precision                :: nh(ixi^s), log_nh(ixi^s)
        double precision                :: eint_val, eint_in
        integer :: ix^d
 
        if (fix_small_values) then
            call ffhd_handle_small_values(.false., w, x, ixi^l, ixo^l, &
                'ffhd_to_primitive_LTE')
        end if
 
        call eos%get_nH(w, x, ixi^l, ixo^l, nh)
        if (eos%ionE) log_nh(ixo^s) = dlog10(nh(ixo^s))
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            w(ix^d,mom(1)) = w(ix^d,mom(1))*inv_rho
            if (ffhd_energy) then
                if (eos%ionE) then
                    eint_val = w(ix^d,e_) - half*w(ix^d,rho_)*w(ix^d,mom(1))**2
                    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
                        w(ix^d,p_) = eos%gamma_minus_1 &
                            * (eint_val - eos%eion_per_nH * nh(ix^d))
                    else
                        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_)*w(ix^d,mom(1))**2)
                end if
            end if
        {end do\}
    end subroutine ffhd_to_primitive_lte
 
    !> LTE thermal pressure from the stored aux state: p = nH(1+He+ne/nH)*Te.
    !> Mirrors the LTE branch of hd_get_pthermal (small-value guards from
    !> ffhd_get_pthermal_origin).
    subroutine ffhd_get_pthermal_lte(w, x, ixI^L, ixO^L, pth)
        use mod_global_parameters
        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)
 
        if (ffhd_energy) then
            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
            pth(ixo^s) = ffhd_adiab * w(ixo^s, rho_)**eos%gamma
        end if
 
        if (fix_small_values) then
            {do ix^db= ixo^lim^db\}
               if(pth(ix^d)<small_pressure) pth(ix^d)=small_pressure
            {end do^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 ffhd_get_pthermal_LTE"
                 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
            {end do^D&\}
        end if
    end subroutine ffhd_get_pthermal_lte
 
    !> cs2 = Gamma_1(LTE) * p/rho.
    subroutine ffhd_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
 
        call ffhd_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\}
    end subroutine ffhd_get_csound2_lte
 
    !> Effective Gamma_1 for LTE+IonE (dispatch on gamma1_method / eos%method).
    subroutine ffhd_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 ffhd_get_gamma1_lte
 
    !> Conserved (rho, rho*v, E) -> prolong form (rho, v, T); T in the e_ slot.
    !> Interpolating in (rho, v, T) avoids Jensen's inequality across the
    !> ionisation plateau. Mirrors hd_to_prolong_LTE.
    subroutine ffhd_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)
 
            w(ix^d,mom(1)) = w(ix^d,mom(1))*inv_rho
 
            eint_val = w(ix^d, e_) - half*w(ix^d, rho_)*w(ix^d,mom(1))**2
            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)
 
            if (eint_val * inv_rho > eos%eint_rho_FI_threshold) then
                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
                call saha_t_from_nh_eint(nh_val, eint_val / nh_val, t_loc, y_loc)
                w(ix^d, e_) = t_loc
            else
                w(ix^d, e_) = t_from_nh_eint(log_nh, &
                    dlog10(eint_val) - log_nh)
            end if
        {end do\}
    end subroutine ffhd_to_prolong_lte
 
    !> Prolong form (rho, v, T) -> conserved (rho, rho*v, E); T read from e_.
    !> Mirrors hd_from_prolong_LTE.
    subroutine ffhd_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
 
        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_)
            nh_val = w(ix^d, rho_) / eos%nH2rhoFactor
            log_nh = dlog10(nh_val)
 
            if (t_val > t_fi) then
                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
                eint_val = saha_eint_from_nh_t(nh_val, t_val) * nh_val
            else
                eint_val = eint_nh_from_t(log_nh, &
                    dlog10(max(t_val, 10.0d0**log_t_min))) * nh_val
            end if
 
            w(ix^d, e_) = eint_val + half*w(ix^d, rho_)*w(ix^d,mom(1))**2
            w(ix^d,mom(1)) = w(ix^d,rho_)*w(ix^d,mom(1))
        {end do\}
    end subroutine ffhd_from_prolong_lte
 
    !=========================================================================
    !> PI energy-mode (eos_type='PI', ionE=.true.) thermodynamics for FFHD.
    !> eint carries the ionisation-energy term, so p=(gamma-1)*eint no longer
    !> holds; the eint<->p relation is delegated to the portable scalar backend
    !> (mod_eos_PI). Single momentum, total-energy only. Mirrors the
    !> hd PI-energy family.
    !=========================================================================
 
    !> Primitive pressure -> total energy via backend (mom(1) is velocity here).
    subroutine ffhd_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 (ffhd_energy) then
                call eint_from_rho_p_pi(w(ix^d,rho_), w(ix^d,p_), eint)
                w(ix^d,e_) = eint + half*w(ix^d,mom(1))**2*w(ix^d,rho_)
            end if
        {end do\}
    end subroutine ffhd_p_to_e_pi
 
    !> Primitive -> conserved (PI energy).
    subroutine ffhd_to_conserved_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)
 
        call ffhd_p_to_e_pi(ixi^l, ixo^l, w, x)
        w(ixo^s,mom(1)) = w(ixo^s,rho_)*w(ixo^s,mom(1))
    end subroutine ffhd_to_conserved_pi
 
    !> Conserved -> primitive (PI energy): KE removed -> eint, backend -> p.
    subroutine ffhd_to_primitive_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 :: inv_rho, eint_val, t, rfac
        integer :: ix^d
 
        if (fix_small_values) then
            call ffhd_handle_small_values(.false., w, x, ixi^l, ixo^l, &
                'ffhd_to_primitive_PI')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            w(ix^d,mom(1)) = w(ix^d,mom(1))*inv_rho
            if (ffhd_energy) then
                eint_val = w(ix^d,e_) - half*w(ix^d,rho_)*w(ix^d,mom(1))**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\}
    end subroutine ffhd_to_primitive_pi
 
    !> PI-energy thermal pressure: eint (= e-KE) -> p via the backend.
    subroutine ffhd_get_pthermal_pi(w, x, ixI^L, ixO^L, pth)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:ndim)
        double precision, intent(out):: pth(ixi^s)
        double precision :: ei(ixi^s), tpi, rpi
        integer :: ix^d
 
        ei(ixo^s) = ffhd_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\}
    end subroutine ffhd_get_pthermal_pi
 
    !> Adiabatic sound speed^2 from primitive (rho, p) via the backend.
 
    !> Effective Gamma_1 = cs2 * rho / p (PI energy).
 
    !> PI-energy Te_ refresh: eint (= e-KE) -> T via the backend.
 
    !> Prominence (T,p table) R-factor from (rho, pth) -- no-energy only, so
    !> pth=(gamma-1)*eint and R is independent of it. Mirrors hd_Rfactor_prominence.
 
    !> Prominence (T,p) Te_ update from (rho, pth).
 
end module mod_ffhd_eos