riemannflux2d_LinearEuler2D_t Function

  pure function riemannflux2d_LinearEuler2D_t(this,sL,sR,dsdx,nhat) result(flux)
    !! Characteristic-decomposition (impedance-matched) interface flux for
    !! linear acoustics with possibly discontinuous sound speed.
    !!
    !! The normal-flux Jacobian has eigenstructure
    !!   +c : right-going acoustic mode, W_+ = rho0*c*u_n + p
    !!   -c : left-going  acoustic mode, W_- = -rho0*c*u_n + p
    !!    0 : entropy density mode      , W_0 = rho' - p/c^2
    !!    0 : tangential vorticity mode , u_t
    !! Upwinding each mode by its characteristic direction at the face,
    !!   W_+|* = W_+|_L   (transmitted from left,  at speed +c_L)
    !!   W_-|* = W_-|_R   (transmitted from right, at speed -c_R)
    !! yields the impedance-matched interface state
    !!   u_n* = (Z_L u_n,L + Z_R u_n,R + (p_L - p_R)) / (Z_L + Z_R)
    !!   p*   = (Z_R p_L + Z_L p_R + Z_L Z_R (u_n,L - u_n,R)) / (Z_L + Z_R)
    !! with the per-side acoustic impedance Z = rho0*c (each side using its own
    !! background density rho0 and sound speed c). This is exact upwind /
    !! Godunov for the linearised acoustic system and reduces correctly to
    !! Fresnel reflection / transmission across an impedance jump (Z_L .ne. Z_R,
    !! whether from a density jump, a sound-speed jump, or both). LLF with
    !! cmax = max(c_L, c_R) over-dissipates tangential and entropy modes
    !! and at high polynomial order fails to stably handle the
    !! impedance mismatch (aliasing instability at material interfaces).
    !!
    !! The reconstructed density/momentum/pressure fluxes need a single rho0
    !! (and c^2) at the face, but rho0 and c are two-valued across a material
    !! interface. We use the arithmetic averages rho0_avg and c2_avg; this is a
    !! pragmatic treatment of the non-conservative products p/rho0 and
    !! rho0*c^2*div(v) at a face where the coefficients jump. A fully
    !! path-conservative (Castro-Pares) treatment would use each side's own
    !! coefficients in its surface integral. For piecewise-constant material
    !! regions the interior is exactly entropy-conservative and the impedance
    !! solver above carries the (entropy-stable) interface dissipation.
    class(LinearEuler2D_t),intent(in) :: this
    real(prec),intent(in) :: sL(1:this%nvar)
    real(prec),intent(in) :: sR(1:this%nvar)
    real(prec),intent(in) :: dsdx(1:this%nvar,1:2)
    real(prec),intent(in) :: nhat(1:2)
    real(prec) :: flux(1:this%nvar)
    ! Local
    real(prec) :: rho0L,rho0R,rho0_avg,cL,cR,ZL,ZR,unL,unR,pL,pR,un_star,p_star,c2_avg

    rho0L = sL(6)
    rho0R = sR(6)
    rho0_avg = 0.5_prec*(rho0L+rho0R)
    cL = sL(5)
    cR = sR(5)
    ZL = rho0L*cL
    ZR = rho0R*cR

    unL = sL(2)*nhat(1)+sL(3)*nhat(2)
    unR = sR(2)*nhat(1)+sR(3)*nhat(2)
    pL = sL(4)
    pR = sR(4)

    un_star = (ZL*unL+ZR*unR+(pL-pR))/(ZL+ZR)
    p_star = (ZR*pL+ZL*pR+ZL*ZR*(unL-unR))/(ZL+ZR)
    c2_avg = 0.5_prec*(cL*cL+cR*cR)

    flux(1) = rho0_avg*un_star
    flux(2) = p_star*nhat(1)/rho0_avg
    flux(3) = p_star*nhat(2)/rho0_avg
    flux(4) = rho0_avg*c2_avg*un_star
    flux(5) = 0.0_prec
    flux(6) = 0.0_prec
    if(.false.) flux(1) = flux(1)+dsdx(1,1) ! suppress unused-dummy-argument warning

  endfunction riemannflux2d_LinearEuler2D_t