mod_eos_pi.t

Go to the documentation of this file.

!=============================================================================
!> PI (partial-ionisation, eos_type='PI') arm of the eos% family.
!>
!> The adapter between the ionisation backend (mod_eos_pi_tables) and the eos%
!> authority. PI rides the FI fully-ionised ideal-gas closure (RR=1
!> normalisation, p=(gamma-1) eint when ionE=F) and differs only through the
!> variable mean-mass that enters via eos%get_Rfactor and -- in energy mode --
!> the ionisation energy folded into eint. mod_eos_pi_tables has no mod_eos
!> dependency, so eos% -> ionisation is acyclic.
!>
!> Structure mirrors mod_eos_fi / mod_eos_lte:
!>   A. init / finalise               -- wire the eos% pointer targets
!>   B. generic block getters         -- the eos%-interface routines bound for PI
!>                                       (csound2, gamma1, Rfactor, Te update);
!>                                       module-agnostic via iw_rho/iw_e/iw_te,
!>                                       so the hd/mhd/ffhd seams share ONE copy
!>                                       instead of three. The seam keeps only
!>                                       the B/KE-specific conversions.
!>   C. scalar combined-solve backend -- one ionisation solve yields T,p,R(,iz)
!>                                       together; called per-cell by the seam
!>                                       conversions and by section D.
!>   D. fl-port scalar shims          -- f(log_nH,log_x) callbacks for the
!>                                       cooling / conduction port objects.
!> mod_eos re-exports the public names so the seams reach them through the single
!> `use mod_eos` facade, exactly as for the FI/LTE kernels.
!=============================================================================
module mod_eos_pi
    use mod_global_parameters
    use mod_eos_container, only: eos, eos_type_pi
    use mod_eos_fi, only: get_temperature_from_eint_fi, &
         get_temperature_from_pressure_fi
    use mod_eos_pi_tables
    use mod_comm_lib, only: mpistop
    implicit none
    private
 
    !> Lifecycle (PI arms of the eos_init / eos_finalise dispatchers). The eos%
    !> block getters (update_eos_PI, get_Te_PI, get_csound2_PI, get_Rfactor_*_PI)
    !> are PRIVATE: they are bound to eos% pointers inside eos_finalise_PI, so no
    !> external module needs to name them.
    public :: eos_init_pi, eos_finalise_pi
 
    !> get_gamma1_PI is bound to the phys_get_gamma1 pointer in the seam (the
    !> phys layer owns that pointer), so it must stay public.
    public :: get_gamma1_pi
 
    !> Scalar combined-solve backend, called per-cell by the seam conversions
    !> (rho = nH = 1 per H under the FI normalisation; eint is the GAS internal
    !> energy, carrying the ionisation energy in energy mode).
    public :: state_from_eint_pi   !> rho,eint -> T, p, R   (energy mode)
    public :: p_eint_from_rho_t_pi !> rho,T    -> p, eint, R
    public :: eint_from_rho_p_pi   !> rho,p    -> eint      (prim -> cons)
    public :: csound2_prim_pi      !> rho,p    -> csound2
 
    !> fl-port scalar callbacks (signature f(log_nH, log_x) -> scalar) for the
    !> cooling / conduction port objects. PI's T-only ionisation is
    !> nH-independent, so log_nH is ignored and quantities are evaluated per H.
    public :: eint_from_t_pi       !> (log_nH, log_T)      -> eint/nH
    public :: p2eint_pi            !> (log_nH, log_p_nH)   -> eint/p factor
    public :: t_from_eint_pi       !> (log_nH, log_eint_nH)-> T
    public :: y_from_eint_pi       !> (log_nH, log_eint_nH)-> ne/nH
 
contains
 
    !=========================================================================
    ! A. Lifecycle
    !=========================================================================
 
    !> PI arm of eos_init (before units are known): PI shares the FI
    !> temperature-from-pressure getter -- T=p/(R*rho) routed through
    !> eos%get_Rfactor, which the seam points at the ionisation R-factor.
    subroutine eos_init_pi()
        eos%get_temperature_from_pressure => get_temperature_from_pressure_fi
    end subroutine eos_init_pi
 
    !> PI arm of eos_finalise (after unit_* are set): ride FI's ideal-gas getter
    !> set and RR=1 normalisation (the variable mean-mass enters only via
    !> eos%get_Rfactor, wired by the seam link arm), then initialise the
    !> ionisation backend. In energy mode eint carries the ionisation energy, so
    !> swap in the backend eint->T inversion.
    subroutine eos_finalise_pi()
        use mod_global_parameters, only: mype
        double precision :: rfactor_norm
 
        eos%update_eos => update_eos_pi   !> direct Te cache refresh each substep
        eos%get_temperature_from_eint => get_temperature_from_eint_fi
        eos%get_Te => get_te_pi           !> read the cached Te (refreshed above)
 
        !> R-factor: dispatch on pi_table. Physics-independent, so bound here
        !> (was in the hd/mhd/ffhd link arms) -> these targets stay private.
        if (eos%pi_table == 'prominence') then
            eos%get_Rfactor => get_rfactor_prominence_pi
        else
            eos%get_Rfactor => get_rfactor_tonly_pi
        end if
        !> Fully-ionised particle counts (2 + 3*A_He per H; ne/nH = 1 + 2*A_He)
        eos%n_per_nH_FI = 2.0d0 + 3.0d0 * eos%He_abundance
        eos%neOnH_FI    = 1.0d0 + 2.0d0 * eos%He_abundance
 
        !> Initialise the ionisation backend. Under the FI normalisation the
        !> (a,b) factors are absorbed into unit_*, RR=1, and b=(2+3 A_He) is the
        !> full-ionisation reference in unit_pressure; the module's
        !> R = (1+iz_H+A_He(...))/Rfactor_norm must reduce to 1 at full ionisation,
        !> so Rfactor_norm = 2+3 A_He (the LTE a=b=1 value 1+4 A_He would be wrong).
        rfactor_norm = 2.0d0 + 3.0d0 * eos%He_abundance
        call ionization_degree_init(eos%He_abundance, rfactor_norm, &
             trim(eos%pi_table), include_energy=eos%ionE)
        if (mype == 0) write(*,*) "[eos PI] ionisation backend initialised: table=", &
             trim(eos%pi_table), " ionE=", eos%ionE, " Rfactor_norm=", rfactor_norm
 
        !> Energy mode: eint carries the ionisation energy, so swap in the
        !> backend eint->T inversion and the energy-mode csound2 (both
        !> physics-independent; get_csound2_PI was in the link arms).
        if (eos%ionE) then
            eos%get_temperature_from_eint => get_temperature_from_eint_pi
            eos%get_csound2 => get_csound2_pi
        end if
    end subroutine eos_finalise_pi
 
    !=========================================================================
    ! B. Generic block getters (eos%-interface shape).
    !    Module-agnostic: the primitive pressure is w(iw_e) (e_ aliases p_ in
    !    primitive state) and the cached temperature is w(iw_te), published by
    !    the phys module. One copy serves hd, mhd and ffhd.
    !=========================================================================
 
    !> Adiabatic sound speed squared from primitive (rho, p). PI energy mode.
    subroutine get_csound2_pi(w, x, ixI^L, ixO^L, cs2)
        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
 
        {do ix^db=ixomin^db,ixomax^db\}
            call csound2_prim_pi(w(ix^d,iw_rho), w(ix^d,iw_e), cs2(ix^d))
        {end do\}
    end subroutine get_csound2_pi
 
    !> Effective Gamma1 = cs2 * rho / p for the same primitive state.
    subroutine get_gamma1_pi(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)
 
        double precision :: cs2
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            call csound2_prim_pi(w(ix^d,iw_rho), w(ix^d,iw_e), cs2)
            gamma1(ix^d) = cs2 * w(ix^d,iw_rho) / w(ix^d,iw_e)
        {end do\}
    end subroutine get_gamma1_pi
 
    !> No-energy R-factor (chromosphere/flare T-only tables) from the stored Te_.
    !> Denominator (2+3 He) = b (full-ion reference) absorbed into unit_pressure
    !> under the FI normalisation, RR=1 -> R->1 at full ionisation.
    subroutine get_rfactor_tonly_pi(w, x, ixI^L, ixO^L, Rfactor)
        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 :: iz_h(ixo^s), iz_he(ixo^s)
 
        call ionization_degree_from_temperature(ixi^l,ixo^l,w(ixi^s,iw_te),iz_h,iz_he)
        rfactor(ixo^s) = (1.d0 + iz_h(ixo^s) + eos%He_abundance &
             * (1.d0 + iz_he(ixo^s)*(1.d0+iz_he(ixo^s)))) / (2.d0 + 3.d0*eos%He_abundance)
    end subroutine get_rfactor_tonly_pi
 
    !> Prominence (T,p) R-factor: recompute from (rho, pth) via the backend's
    !> prominence inversion (no-energy only -> pth=(gamma-1)*eint, R-independent).
    subroutine get_rfactor_prominence_pi(w, x, ixI^L, ixO^L, Rfactor)
        use mod_physics, only: phys_get_ei
        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 :: ei(ixi^s), pth(ixi^s), rho(ixi^s), t(ixi^s)
 
        ei(ixo^s)  = phys_get_ei(w, ixi^l, ixo^l)
        pth(ixo^s) = eos%gamma_minus_1 * ei(ixo^s)
        rho(ixo^s) = w(ixo^s,iw_rho)
        call ionization_get_state(ixi^l, ixo^l, rho, pth, t, rfactor)
    end subroutine get_rfactor_prominence_pi
 
    !> Direct (no-lag) temperature refresh: recompute the cached w(iw_te) from the
    !> CURRENT conserved state, mirroring update_eos_LTE. Energy mode inverts the
    !> gas internal energy (incl. ionisation energy) by Newton; no-energy modes
    !> form pth=(gamma-1)*eint and invert (rho,pth)->T via the backend
    !> (y-bisection for chromosphere/flare, p-bisection for prominence). Shared by
    !> update_eos_PI (eos%update_eos) and the gated source/STS call sites, refresh Te
    subroutine update_eos_pi(ixI^L, ixO^L, w, x)
        use mod_physics, only: phys_get_ei
        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 :: ei(ixi^s), pth(ixi^s), rho(ixi^s), t(ixi^s), rf(ixi^s)
        double precision :: tsc, psc, rsc
        integer :: ix^d
 
        ei(ixo^s) = phys_get_ei(w, ixi^l, ixo^l)
        if (eos%ionE) then
            !> energy mode (chromosphere/flare): eint carries the ionisation
            !> energy, invert eint->T by Newton per cell.
            {do ix^db=ixomin^db,ixomax^db\}
                call state_from_eint_pi(w(ix^d,iw_rho), ei(ix^d), tsc, psc, rsc)
                w(ix^d,iw_te) = tsc
            {end do\}
        else
            !> no-energy: pth=(gamma-1)*eint; ionization_get_state inverts
            !> (rho,pth)->T (T-only y-bisection or prominence p-bisection).
            pth(ixo^s) = eos%gamma_minus_1 * ei(ixo^s)
            rho(ixo^s) = w(ixo^s,iw_rho)
            call ionization_get_state(ixi^l, ixo^l, rho, pth, t, rf)
            w(ixo^s,iw_te) = t(ixo^s)
        end if
    end subroutine update_eos_pi
 
    !> eos%get_Te arm: read the cached temperature (refreshed by update_eos_PI /
    !> the gated source/STS sites), mirroring get_Te_LTE.
    subroutine get_te_pi(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)
        t(ixo^s) = w(ixo^s,iw_te)
    end subroutine get_te_pi
 
    !> PI energy mode: temperature from GAS internal energy via the backend
    !> eint->T inversion (eint carries the ionisation-energy term, so the FI
    !> p=(gamma-1)*eint relation does not hold). w(iw_e) is assumed to be eint
    !> already (the etot hook subtracts KE/ME via phys_e_to_ei before calling).
    subroutine get_temperature_from_eint_pi(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 :: pth, rfac
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            call state_from_eint_pi(w(ix^d,iw_rho), w(ix^d,iw_e), &
                 res(ix^d), pth, rfac)
        {end do\}
    end subroutine get_temperature_from_eint_pi
 
    !=========================================================================
    ! C. Scalar combined-solve backend (eos%inv_gamma_minus_1 throughout,
    !    matching the FI ideal-gas split; the ionisation-energy term is folded
    !    in by the ionisation module when ionE=.true.).
    !=========================================================================
 
    !> Conserved gas internal energy -> temperature, thermal pressure, R.
    !> Inverse of p_eint_from_rho_T_PI. eint is the GAS internal energy
    !> (mechanical energy already removed by the caller).
    subroutine state_from_eint_pi(rho, eint, T, p, Rfactor, iz_H, iz_He)
        double precision, intent(in)  :: rho, eint
        double precision, intent(out) :: t, p, rfactor
        double precision, intent(out), optional :: iz_h, iz_he
 
        if (eos%type_id /= eos_type_pi) call mpistop( &
            "state_from_eint_PI called outside eos_type='PI'")
        call ionization_get_state_from_eint(rho, eint, eos%inv_gamma_minus_1, &
             t, p, rfactor, iz_h, iz_he)
    end subroutine state_from_eint_pi
 
    !> Temperature -> thermal pressure and GAS internal energy (incl. ion energy)
    subroutine p_eint_from_rho_t_pi(rho, T, p, eint, Rfactor, iz_H, iz_He)
        double precision, intent(in)  :: rho, t
        double precision, intent(out) :: p, eint, rfactor
        double precision, intent(out), optional :: iz_h, iz_he
 
        if (eos%type_id /= eos_type_pi) call mpistop( &
            "p_eint_from_rho_T_PI called outside eos_type='PI'")
        call ionization_get_p_eint_from_rho_t(rho, t, eos%inv_gamma_minus_1, &
             p, eint, rfactor, iz_h, iz_he)
    end subroutine p_eint_from_rho_t_pi
 
    !> Primitive pressure -> GAS internal energy (prim -> conserved direction):
    !> invert (rho,p)->T then map T->eint, so the forward/inverse pair stay consistent.
    subroutine eint_from_rho_p_pi(rho, p, eint)
        double precision, intent(in)  :: rho, p
        double precision, intent(out) :: eint
        double precision :: t, pcheck, rfactor
 
        if (eos%type_id /= eos_type_pi) call mpistop( &
            "eint_from_rho_p_PI called outside eos_type='PI'")
        call ionization_get_state_scalar(rho, p, t, rfactor)
        call ionization_get_p_eint_from_rho_t(rho, t, eos%inv_gamma_minus_1, &
             pcheck, eint, rfactor)
    end subroutine eint_from_rho_p_pi
 
    !> Adiabatic sound speed squared from primitive (rho, p).
    subroutine csound2_prim_pi(rho, p, csound2)
        double precision, intent(in)  :: rho, p
        double precision, intent(out) :: csound2
        double precision :: t, rfactor
 
        if (eos%type_id /= eos_type_pi) call mpistop( &
            "csound2_prim_PI called outside eos_type='PI'")
        call ionization_get_state_scalar(rho, p, t, rfactor)
        call ionization_get_csound2_t(t, eos%inv_gamma_minus_1, csound2)
    end subroutine csound2_prim_pi
 
    !=========================================================================
    ! D. fl-port scalar callbacks (rho = nH = 1 per H). Thin shims over the
    !    section-C backend matching the cooling/conduction port signatures.
    !=========================================================================
 
    !> Internal energy per H from temperature: eint/nH(T).
    double precision function eint_from_t_pi(log_nH, log_T) result(eint)
        double precision, intent(in) :: log_nh, log_t
        double precision :: p, rfactor
        call p_eint_from_rho_t_pi(1.d0, 10.d0**log_t, p, eint, rfactor)
    end function eint_from_t_pi
 
    !> eint/p factor from pressure per H: maps p -> eint = p * (this).
    double precision function p2eint_pi(log_nH, log_p_nH) result(factor)
        double precision, intent(in) :: log_nh, log_p_nh
        double precision :: p, eint
        p = 10.d0**log_p_nh
        call eint_from_rho_p_pi(1.d0, p, eint)
        factor = eint / p
    end function p2eint_pi
 
    !> Temperature from internal energy per H.
    double precision function t_from_eint_pi(log_nH, log_eint_nH) result(T)
        double precision, intent(in) :: log_nh, log_eint_nh
        double precision :: p, rfactor
        call state_from_eint_pi(1.d0, 10.d0**log_eint_nh, t, p, rfactor)
    end function t_from_eint_pi
 
    !> Electron-to-hydrogen ratio ne/nH from internal energy per H.
    !> ne/nH = iz_H + A_He*iz_He*(1+iz_He) (matches the R-factor numerator's
    !> electron count; second He ionisation assumed equal to the first).
    double precision function y_from_eint_pi(log_nH, log_eint_nH) result(y)
        double precision, intent(in) :: log_nh, log_eint_nh
        double precision :: t, p, rfactor, iz_h, iz_he
        call state_from_eint_pi(1.d0, 10.d0**log_eint_nh, t, p, rfactor, iz_h, iz_he)
        y = iz_h + eos%He_abundance*iz_he*(1.d0+iz_he)
    end function y_from_eint_pi
 
end module mod_eos_pi