mod_eos_fi.t

Go to the documentation of this file.

!=============================================================================
!> FI (fully-ionised, constant-gamma) EoS kernels for the eos% family.
!>
!> Carved out of mod_eos.t (the per-type split). The fully-ionised ideal gas:
!> p = (gamma-1) eint, constant Gamma_1 = gamma, temperature T = p/(R rho) with
!> R the constant ionisation R-factor. update_eos is a no-op (nothing cached).
!> mod_eos re-exports the public names so existing `use mod_eos` callers are
!> unaffected; eos_finalise wires these as the FI pointer targets.
!=============================================================================
module mod_eos_fi
    use mod_global_parameters
    use mod_eos_container
    use mod_eos_shared_functions, only: get_rho
 
    implicit none
    private
 
    !> get_Te_FI is PRIVATE: bound to eos%get_Te inside eos_finalise_FI only.
    !> The rest stay public: update_eos_FI (ffhd_phys), get_gamma1_FI (the seam
    !> phys_get_gamma1), get_temperature_from_{eint,pressure}_FI (reused by
    !> mod_eos_pi), and the init/finalise dispatcher arms.
    public :: update_eos_fi
    public :: get_gamma1_fi
    public :: get_temperature_from_eint_fi
    public :: get_temperature_from_pressure_fi
    public :: eos_init_fi
    public :: eos_finalise_fi
 
contains
 
    !> FI arm of eos_init (before units are known): wire the FI
    !> temperature-from-pressure getter, which must be live before
    !> (m)hd_link_eos captures it into the radiation fluid object.
    subroutine eos_init_fi()
        eos%get_temperature_from_pressure => get_temperature_from_pressure_fi
    end subroutine eos_init_fi
 
    !> FI arm of eos_finalise: wire the ideal-gas getter set and the fully-ionised
    !> particle counts (2 + 3*A_He per H; ne/nH = 1 + 2*A_He). No tables to load.
    subroutine eos_finalise_fi()
        eos%update_eos => update_eos_fi   !> nothing cached for ideal gas
        eos%get_temperature_from_eint => get_temperature_from_eint_fi
        eos%get_Te => get_te_fi
        eos%n_per_nH_FI = 2.0d0 + 3.0d0 * eos%He_abundance
        eos%neOnH_FI    = 1.0d0 + 2.0d0 * eos%He_abundance
    end subroutine eos_finalise_fi
 
    !> FI update_eos: nothing to cache for an ideal gas. Called each RK substep.
    subroutine update_eos_fi(ixI^L, ixO^L, w, x)
        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)
    end subroutine update_eos_fi
 
    !> Gamma_1 for the fully-ionised ideal gas: constant eos%gamma.
    subroutine get_gamma1_fi(w, x, ixI^L, ixO^L, gamma1)
        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)
 
        gamma1(ixo^s) = eos%gamma
 
    end subroutine get_gamma1_fi
 
    !> Cached-Te accessor for FI: T = pth/(rho*R) (no stored Te_; recompute).
    subroutine get_te_fi(w,x,ixI^L,ixO^L,T)
        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)
        double precision :: rfactor(ixi^s), pth(ixi^s), rho(ixi^s)
 
        !> No timing here: get_thermal_pressure is already timed.
        !> The division and Rfactor call are trivial.
        !> get_thermal_pressure already returns the total pressure under
        !> equilibrium splitting, so only the density needs the background.
        call eos%get_thermal_pressure(w, x, ixi^l, ixo^l, pth)
        call eos%get_Rfactor(w,x,ixi^l,ixo^l,rfactor)
        call get_rho(w,x,ixi^l,ixo^l,rho)
        t(ixo^s) = pth(ixo^s) / (rho(ixo^s) * rfactor(ixo^s))
 
    end subroutine get_te_fi
 
    !> Temperature from internal energy, fully-ionised ideal gas: T = pth/(rho*R).
    subroutine get_temperature_from_eint_fi(w, x, ixI^L, ixO^L, res)
        !> Assumes input energy is internal energy
        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 :: rfactor(ixi^s), rho(ixi^s), pth(ixi^s)
 
        call eos%get_Rfactor(w,x,ixi^l,ixo^l,rfactor)
        call get_rho(w,x,ixi^l,ixo^l,rho)
        !> Under equilibrium splitting the e slot holds the perturbation internal
        !> energy, so the background pressure is added back to recover the total.
        pth(ixo^s) = eos%gamma_minus_1 * w(ixo^s,iw_e)
        if (iw_equi_p > 0) pth(ixo^s) = pth(ixo^s) &
             + block%equi_vars(ixo^s,iw_equi_p,b0i)
        res(ixo^s) = pth(ixo^s) / (rfactor(ixo^s) * rho(ixo^s))
 
    end subroutine get_temperature_from_eint_fi
 
    !> FI temperature from primitive pressure: T = p / (R * rho).
    !> w is PRIMITIVE here, so iw_e holds the thermal pressure. Lives in the EoS
    !> so hd and mhd share one routine (replaces the orphaned, never-called
    !> hd_get_temperature_from_prim that used to sit in mod_hd_phys).
    subroutine get_temperature_from_pressure_fi(w, x, ixI^L, ixO^L, res)
        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 :: rfactor(ixi^s)
 
        call eos%get_Rfactor(w, x, ixi^l, ixo^l, rfactor)
        res(ixo^s) = w(ixo^s,iw_e) / (rfactor(ixo^s) * w(ixo^s,iw_rho))
    end subroutine get_temperature_from_pressure_fi
 
end module mod_eos_fi