fix_prim_state Subroutine

private pure subroutine fix_prim_state(u)

Make a reconstructed primitive face state admissible before it is handed to the Riemann solver.

Two things happen here, and only the second is optional. The MUSCL reconstruction limits every slot of the state independently, including the auxiliaries xi and lfac -- but those are not independent variables, they are functions of (rho, u^i, p). A reconstructed face state therefore does not in general satisfy lfac = sqrt(1+|u|^2) or xi = lfac^2rhoh, while get_flux, get_cmax and estimate_speeds_minmax all read them. Left alone, v2 = 1-1/lfac^2 and v_n = u_n/lfac come from different sources, v_n can exceed v (even exceed 1), the radicand of the relativistic characteristic speeds goes negative and the HLL wave speeds collapse onto each other -- which removes the solver's dissipation exactly where it is needed. So the auxiliaries are always recomputed from the reconstructed primitives, by the same expressions to_conservative uses.

The floors are the optional part. Both default to zero, in which case the max() calls only clamp states that were already unphysical; set srhd_small_pressure (and srhd_small_density) positive in &srhd_list to get an actual floor, which is the useful configuration -- a state floored to exactly zero still gives csound2 = gammap/(rhoh) a zero denominator.

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(inout) :: u(nw_phys)

Called by

proc~~fix_prim_state~~CalledByGraph proc~fix_prim_state fix_prim_state proc~muscl_reconstruct_prim muscl_reconstruct_prim proc~muscl_reconstruct_prim->proc~fix_prim_state proc~reconflux_muscl_hll_prim reconflux_muscl_hll_prim proc~reconflux_muscl_hll_prim->proc~muscl_reconstruct_prim proc~reconflux_muscl_hllc_prim reconflux_muscl_hllc_prim proc~reconflux_muscl_hllc_prim->proc~muscl_reconstruct_prim proc~reconflux_muscl_llf_prim reconflux_muscl_llf_prim proc~reconflux_muscl_llf_prim->proc~muscl_reconstruct_prim proc~finite_volume_local_muscl_hll finite_volume_local_muscl_hll proc~finite_volume_local_muscl_hll->proc~reconflux_muscl_hll_prim proc~finite_volume_local_muscl_hllc finite_volume_local_muscl_hllc proc~finite_volume_local_muscl_hllc->proc~reconflux_muscl_hllc_prim proc~finite_volume_local_muscl_llf finite_volume_local_muscl_llf proc~finite_volume_local_muscl_llf->proc~reconflux_muscl_llf_prim proc~finite_volume_local finite_volume_local proc~finite_volume_local->proc~finite_volume_local_muscl_hll proc~finite_volume_local->proc~finite_volume_local_muscl_hllc proc~finite_volume_local->proc~finite_volume_local_muscl_llf

Source Code

  pure subroutine fix_prim_state(u)
    
    real(dp), intent(inout) :: u(nw_phys)

    real(dp) :: rho, pth, rhoh, E, E_th

    rho = max(u(iw_rho), srhd_small_density)
    pth = max(u(iw_e),   srhd_small_pressure)
    u(iw_rho) = rho
    u(iw_e)   = pth

    u(lfac_) = dsqrt(1.0d0 + u(iw_mom(1))**2 + u(iw_mom(2))**2 &
         + u(iw_mom(3))**2)

    !! begin: call srhd_get_enthalpy_eos(rho,pth,rhoh)
    if (srhd_eos) then
       E_th = pth*inv_gamma_1
       E    = E_th + dsqrt(E_th**2 + rho**2)
       ! writing rho/E on purpose, for numerics
       rhoh = 0.5_dp*((srhd_gamma + 1.0_dp)*E - gamma_1*rho*(rho/E))
    else
       rhoh = rho + gamma_to_gamma_1*pth
    end if
    !! end: call srhd_get_enthalpy_eos(rho,pth,rhoh)

    u(xi_) = u(lfac_)**2 * rhoh

  end subroutine fix_prim_state