riemann_hllc_prim Subroutine

private subroutine riemann_hllc_prim(uL, uR, xC, flux_dim, F)

One-face HLLC numerical flux from primitive L/R states (HD only). Uses scalar pressure assumption; not valid for MHD (use HLLD instead). Reference: Toro (2010), chapter 10 (Variant 2) does not yet take flux_type(flux_dim,nw_flux) into account

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(inout) :: uL(nw_phys)
real(kind=dp), intent(inout) :: uR(nw_phys)
real(kind=dp), intent(in) :: xC(ndim)
integer, intent(in) :: flux_dim
real(kind=dp), intent(out) :: F(nw_flux)

Calls

proc~~riemann_hllc_prim~2~~CallsGraph proc~riemann_hllc_prim~2 riemann_hllc_prim proc~estimate_speeds_toro_pvrs~2 estimate_speeds_toro_pvrs proc~riemann_hllc_prim~2->proc~estimate_speeds_toro_pvrs~2 proc~get_flux~2 get_flux proc~riemann_hllc_prim~2->proc~get_flux~2 proc~to_conservative~6 to_conservative proc~riemann_hllc_prim~2->proc~to_conservative~6 tracer tracer proc~get_flux~2->tracer

Called by

proc~~riemann_hllc_prim~2~~CalledByGraph proc~riemann_hllc_prim~2 riemann_hllc_prim proc~reconflux_muscl_hllc_prim~2 reconflux_muscl_hllc_prim proc~reconflux_muscl_hllc_prim~2->proc~riemann_hllc_prim~2 proc~finite_volume_local_muscl_hllc~2 finite_volume_local_muscl_hllc proc~finite_volume_local_muscl_hllc~2->proc~reconflux_muscl_hllc_prim~2 proc~finite_volume_local~2 finite_volume_local proc~finite_volume_local~2->proc~finite_volume_local_muscl_hllc~2

Source Code

  subroutine riemann_hllc_prim(uL, uR, xC, flux_dim, F)
    
    real(dp), intent(inout) :: uL(nw_phys), uR(nw_phys)
    real(dp), intent(in)    :: xC(ndim)
    integer,  intent(in)    :: flux_dim
    real(dp), intent(out)   :: F(nw_flux)

    real(dp) :: flux_l(nw_flux), flux_r(nw_flux)
    real(dp) :: sL, sR, s_star
    real(dp) :: dL, dR, denom
    real(dp) :: rhoL, rhoR, unL, unR, pL, pR
    real(dp) :: PLR
    real(dp) :: Dstar
    real(dp) :: wmax
    real(dp), parameter :: eps = 1e-14_dp
    integer :: i

    call get_flux(uL, xC, flux_dim, flux_l)
    call get_flux(uR, xC, flux_dim, flux_r)

    call estimate_speeds_toro_pvrs(uL, uR, xC, flux_dim, sL, sR)

    ! cache primitive scalars before uL/uR are overwritten to conservative
    rhoL = uL(iw_rho)
    rhoR = uR(iw_rho)
    unL  = uL(iw_mom(flux_dim))
    unR  = uR(iw_mom(flux_dim))
    pL   = uL(iw_e)
    pR   = uR(iw_e)

    ! we re-use these later (Eq. 10.70, Eq. 10.76)
    dL    = rhoL*(sL - unL)
    dR    = rhoR*(sR - unR)
    denom = dL - dR

    ! needed because (Eq. 10.75) uses U_K
    call to_conservative(uL)
    call to_conservative(uR)

    ! fall back to LLF/Rusanov if denom is ~0
    if (abs(denom) <= eps*(abs(dL) + abs(dR) + 1._dp)) then
      wmax = max(abs(sL), abs(sR))
      F = 0.5_dp * ((flux_l + flux_r) - wmax * (uR(1:nw_flux) - uL(1:nw_flux)))
      return
    end if

    ! contact speed S* (Eq. 10.70)
    s_star = ((pR - pL) + unL*dL - unR*dR) / denom

    ! pressure-like term (Eq. 10.76)
    PLR = 0.5_dp * (pL + pR + dL*(s_star - unL) + dR*(s_star - unR))

    ! HLLC flux (Eq. 10.71) + variant 2 (Eq. 10.75)
    if (0._dp .le. sL) then
      F = flux_l
      return
    elseif (0._dp .ge. sR) then
      F = flux_r
      return
    end if

    if (0._dp .le. s_star) then
      ! F_*L
      do i = 1, nw_flux
        Dstar = 0._dp
        if (i .eq. iw_mom(flux_dim)) Dstar = 1._dp
        if (i .eq. iw_e)             Dstar = s_star
        F(i) = (s_star*(sL*uL(i) - flux_l(i)) + sL*PLR*Dstar) / (sL - s_star)
      end do
    else
      ! F_*R
      do i = 1, nw_flux
        Dstar = 0._dp
        if (i .eq. iw_mom(flux_dim)) Dstar = 1._dp
        if (i .eq. iw_e)             Dstar = s_star
        F(i) = (s_star*(sR*uR(i) - flux_r(i)) + sR*PLR*Dstar) / (sR - s_star)
      end do
    end if

  end subroutine riemann_hllc_prim