LMARS (Low-Mach Approximate Riemann Solver, Chen et al. 2013) interface flux. No hydrostatic pressure split: gravity is folded into SourceMethod via the Souza non-conservative form using the geopotential carried in the state vector (variable index 5).
un = 0.5(unL + unR) - (pR - pL) / (2 * rho_bar * c_bar) p = 0.5(pL + pR) - 0.5 * rho_bar * c_bar * (unR - unL)
Upwind on sign of un* selects which side supplies the conserved state for the advective part. Pressure adds to normal momentum. Geopotential (var 5) has zero flux.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(ESAtmo2D_t), | intent(inout) | :: | this |
subroutine BoundaryFlux_ESAtmo2D_t(this)
!! LMARS (Low-Mach Approximate Riemann Solver, Chen et al. 2013)
!! interface flux. No hydrostatic pressure split: gravity is folded
!! into SourceMethod via the Souza non-conservative form using the
!! geopotential carried in the state vector (variable index 5).
!!
!! un* = 0.5*(unL + unR) - (pR - pL) / (2 * rho_bar * c_bar)
!! p* = 0.5*(pL + pR) - 0.5 * rho_bar * c_bar * (unR - unL)
!!
!! Upwind on sign of un* selects which side supplies the conserved
!! state for the advective part. Pressure adds to normal momentum.
!! Geopotential (var 5) has zero flux.
implicit none
class(ESAtmo2D_t),intent(inout) :: this
! Local
integer :: i,k,iel
real(prec) :: nhat(1:2),nmag
real(prec) :: rhoL,uL,vL,thetaL,pL,unL,cL
real(prec) :: rhoR,uR,vR,thetaR,pR,unR,cR
real(prec) :: rho_bar,c_bar,rc,un_star,p_star,gamma
real(prec) :: sL(1:4),sR(1:4),s_up(1:4)
gamma = this%cp/this%cv
do concurrent(i=1:this%solution%N+1,k=1:4,iel=1:this%mesh%nElem)
nhat = this%geometry%nHat%boundary(i,k,iEl,1,1:2)
nmag = this%geometry%nScale%boundary(i,k,iEl,1)
sL = this%solution%boundary(i,k,iel,1:4)
sR = this%solution%extBoundary(i,k,iel,1:4)
rhoL = sL(1)
uL = sL(2)/rhoL
vL = sL(3)/rhoL
thetaL = sL(4)/rhoL
pL = this%p0*(rhoL*this%Rd*thetaL/this%p0)**gamma
unL = uL*nhat(1)+vL*nhat(2)
cL = sqrt(gamma*pL/rhoL)
rhoR = sR(1)
uR = sR(2)/rhoR
vR = sR(3)/rhoR
thetaR = sR(4)/rhoR
pR = this%p0*(rhoR*this%Rd*thetaR/this%p0)**gamma
unR = uR*nhat(1)+vR*nhat(2)
cR = sqrt(gamma*pR/rhoR)
rho_bar = 0.5_prec*(rhoL+rhoR)
c_bar = 0.5_prec*(cL+cR)
rc = rho_bar*c_bar
un_star = 0.5_prec*(unL+unR)-(pR-pL)/(2.0_prec*rc)
p_star = 0.5_prec*(pL+pR)-0.5_prec*rc*(unR-unL)
! Upwind on un_star (advective part of LMARS)
if(un_star >= 0.0_prec) then
s_up = sL
else
s_up = sR
endif
this%flux%boundaryNormal(i,k,iEl,1) = (s_up(1)*un_star)*nmag
this%flux%boundaryNormal(i,k,iEl,2) = (s_up(2)*un_star+p_star*nhat(1))*nmag
this%flux%boundaryNormal(i,k,iEl,3) = (s_up(3)*un_star+p_star*nhat(2))*nmag
this%flux%boundaryNormal(i,k,iEl,4) = (s_up(4)*un_star)*nmag
this%flux%boundaryNormal(i,k,iEl,5) = 0.0_prec
enddo
endsubroutine BoundaryFlux_ESAtmo2D_t