Entropy-stable surface fluxes on 2:1 nonconforming (mortar) edges, following Friedrich, Winters, Del Rey Fernandez, Gassner, Parsani & Carpenter, J. Sci. Comput. 77 (2018), "An entropy stable h/p non-conforming discontinuous Galerkin method with the summation-by-parts property".
Overwrites flux%boundaryNormal on every mortar edge, after BoundaryFlux. With the big edge's nodes i, the nodes j of sub-edge k (in the big edge's orientation), the restriction R_ji = mortarR(i,j,k) (exact interpolation of the big trace to sub-edge node j), the quadrature weights w, the big edge's scaled outward normal N_i = nHat*nScale and c = 1/2 (the sub-edge to big-edge length ratio), the entropy-conserving part is
small edge k, node j : F_kj = sum_i R_ji f#(u_kj, u_i) . (-c N_i) big edge, node i : F_i = (c/w_i) sum_k sum_j R_ji w_j f#(u_i, u_kj) . N_i
where f# is the model's two-point flux (twopointflux2d). They are evaluated as the local flux plus a correction,
F_kj = f#(u_kj,u_kj) . N_kj + sum_i R_ji [f#(u_kj,u_i) - f#(u_kj,u_kj)] . (-c N_i) F_i = f#(u_i,u_i) . N_i + (c/w_i) sum_k sum_j R_ji w_j [f#(u_i,u_kj) - f#(u_i,u_i)] . N_i
with N_kj the small edge's own scaled normal. This equals the form above in exact arithmetic, by the two identities below, and makes the correction vanish identically for a constant state, so free-stream preservation is that of a conforming edge instead of resting on the quadrature identity in floating point. Because sum_i R_ji = 1, the big edge's surface integral equals minus the sum of the small edges', so the coupling is conservative. For a symmetric entropy-conservative f# the entropy flux also telescopes, using sum_k c sum_j w_j R_ji = w_i (exact: the composite quadrature of a degree-N polynomial) and the restriction identity N_kj = -c sum_i R_ji N_i of watertight geometry (the small edge's scaled normal is the restriction of the big edge's). The pair normal is built from the big edge's normals, so both identities hold to round-off whenever the normals are polynomials of degree <= N in the edge coordinate.
Dissipation: on sub-edge node j the model's Riemann flux minus its EC part, d_kj = [ f_R(u_kj, u~_kj, n) - f#(u_kj, u~_kj) . n ] |N~_kj| with u~ = R u (the restricted big trace) and N~_kj = -c sum_i R_ji N_i, is added to the small edge and projected conservatively onto the big edge, -(1/w_i) sum_kj R_ji w_j d_kj. It is entropy stable whenever the entropy variables are linear in the state (a quadratic entropy, e.g. advection or linear shallow water); for a nonlinear entropy its sign is not guaranteed.
Needs solution%MortarExchange after the side exchange (it stages the traces in solution%mortarBuff, slots k: big trace, 2+k: small trace) and exchanges the big edges' normals through mortarNormal so that a small edge's rank has them.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(ECDGModel2D_t), | intent(inout) | :: | this |
subroutine ECMortarFlux_ECDGModel2D_t(this)
!! Entropy-stable surface fluxes on 2:1 nonconforming (mortar) edges, following
!! Friedrich, Winters, Del Rey Fernandez, Gassner, Parsani & Carpenter, J. Sci.
!! Comput. 77 (2018), "An entropy stable h/p non-conforming discontinuous Galerkin
!! method with the summation-by-parts property".
!!
!! Overwrites flux%boundaryNormal on every mortar edge, after BoundaryFlux. With
!! the big edge's nodes i, the nodes j of sub-edge k (in the big edge's
!! orientation), the restriction R_ji = mortarR(i,j,k) (exact interpolation of the
!! big trace to sub-edge node j), the quadrature weights w, the big edge's scaled
!! outward normal N_i = nHat*nScale and c = 1/2 (the sub-edge to big-edge length
!! ratio), the entropy-conserving part is
!!
!! small edge k, node j : F_kj = sum_i R_ji f#(u_kj, u_i) . (-c N_i)
!! big edge, node i : F_i = (c/w_i) sum_k sum_j R_ji w_j f#(u_i, u_kj) . N_i
!!
!! where f# is the model's two-point flux (twopointflux2d). They are evaluated as
!! the local flux plus a correction,
!!
!! F_kj = f#(u_kj,u_kj) . N_kj + sum_i R_ji [f#(u_kj,u_i) - f#(u_kj,u_kj)] . (-c N_i)
!! F_i = f#(u_i,u_i) . N_i + (c/w_i) sum_k sum_j R_ji w_j [f#(u_i,u_kj) - f#(u_i,u_i)] . N_i
!!
!! with N_kj the small edge's own scaled normal. This equals the form above in exact
!! arithmetic, by the two identities below, and makes the correction vanish
!! identically for a constant state, so free-stream preservation is that of a
!! conforming edge instead of resting on the quadrature identity in floating point.
!! Because sum_i R_ji = 1,
!! the big edge's surface integral equals minus the sum of the small edges', so
!! the coupling is conservative. For a symmetric entropy-conservative f# the
!! entropy flux also telescopes, using sum_k c sum_j w_j R_ji = w_i (exact: the
!! composite quadrature of a degree-N polynomial) and the restriction identity
!! N_kj = -c sum_i R_ji N_i of watertight geometry (the small edge's scaled normal
!! is the restriction of the big edge's). The pair normal is built from the big
!! edge's normals, so both identities hold to round-off whenever the normals are
!! polynomials of degree <= N in the edge coordinate.
!!
!! Dissipation: on sub-edge node j the model's Riemann flux minus its EC part,
!! d_kj = [ f_R(u_kj, u~_kj, n) - f#(u_kj, u~_kj) . n ] |N~_kj| with u~ = R u (the
!! restricted big trace) and N~_kj = -c sum_i R_ji N_i, is added to the small edge
!! and projected conservatively onto the big edge, -(1/w_i) sum_kj R_ji w_j d_kj.
!! It is entropy stable whenever the entropy variables are linear in the state
!! (a quadratic entropy, e.g. advection or linear shallow water); for a nonlinear
!! entropy its sign is not guaranteed.
!!
!! Needs solution%MortarExchange after the side exchange (it stages the traces in
!! solution%mortarBuff, slots k: big trace, 2+k: small trace) and exchanges the big
!! edges' normals through mortarNormal so that a small edge's rank has them.
implicit none
class(ECDGModel2D_t),intent(inout) :: this
! Local
real(prec),parameter :: c = 0.5_prec
integer :: N,nv,m,k,i,j,j2,iel,s,eB,sB,eS,sS,flip,offset,rankId,iv
real(prec) :: uB(1:this%nvar),uS(1:this%nvar),uT(1:this%nvar)
real(prec) :: fp(1:this%nvar,1:2),fs(1:this%nvar,1:2),fR(1:this%nvar),acc(1:this%nvar)
real(prec) :: dsdx0(1:this%nvar,1:2)
real(prec) :: NB(1:2),Nt(1:2),nS,nh(1:2),r
real(prec) :: d(1:this%nvar,1:this%solution%N+1,1:2)
if(this%gradient_enabled) then
! The coupling below has no viscous (gradient) contribution, so it would drop
! the viscous flux on mortar edges.
print*,__FILE__,' : Error : EC-DGSEM models with gradients do not support mortar meshes yet.'
error stop 1
endif
N = this%solution%N
nv = this%nvar
rankId = this%mesh%decomp%rankId
offset = this%mesh%decomp%offsetElem(rankId+1)
dsdx0 = 0.0_prec
! Scaled outward normals on every edge, then the mortar staging of the big edges'
! normals (slot k of mortarNormal%mortarBuff on the rank of sub-edge k)
do iel = 1,this%mesh%nElem
do s = 1,4
do i = 1,N+1
this%mortarNormal%boundary(i,s,iel,1:2) = this%geometry%nHat%boundary(i,s,iel,1,1:2)* &
this%geometry%nScale%boundary(i,s,iel,1)
enddo
enddo
enddo
call this%mortarNormal%MortarExchange(this%mesh)
do m = 1,this%mesh%nMortars
! ---- Small edges ----
do k = 1,2
eS = this%mesh%mortarInfo(2*k+1,m)
if(this%mesh%decomp%elemToRank(eS) /= rankId) cycle
sS = this%mesh%mortarInfo(2*k+2,m)/10
flip = this%mesh%mortarInfo(2*k+2,m)-10*sS
do j = 1,N+1
uS = this%solution%mortarBuff(j,2+k,m,1:nv)
j2 = j
if(flip == 1) j2 = N+2-j
fs = this%twopointflux2d(uS,uS)
NB = this%mortarNormal%boundary(j2,sS,eS-offset,1:2)
acc = fs(:,1)*NB(1)+fs(:,2)*NB(2)
Nt = 0.0_prec
uT = 0.0_prec
do i = 1,N+1
r = this%solution%interp%mortarR(i,j,k)
uB = this%solution%mortarBuff(i,k,m,1:nv)
NB = this%mortarNormal%mortarBuff(i,k,m,1:2)
fp = this%twopointflux2d(uS,uB)-fs
acc = acc-c*r*(fp(:,1)*NB(1)+fp(:,2)*NB(2))
Nt = Nt-c*r*NB
uT = uT+r*uB
enddo
nS = sqrt(Nt(1)*Nt(1)+Nt(2)*Nt(2))
nh = Nt/nS
fR = this%riemannflux2d(uS,uT,dsdx0,nh)
fp = this%twopointflux2d(uS,uT)
do iv = 1,nv
this%flux%boundaryNormal(j2,sS,eS-offset,iv) = acc(iv)+ &
(fR(iv)-(fp(iv,1)*nh(1)+fp(iv,2)*nh(2)))*nS
enddo
enddo
enddo
! ---- Big edge ----
eB = this%mesh%mortarInfo(1,m)
if(this%mesh%decomp%elemToRank(eB) /= rankId) cycle
sB = this%mesh%mortarInfo(2,m)
! Dissipation on the sub-edge nodes, from the restricted big trace and normals
do k = 1,2
do j = 1,N+1
uS = this%solution%mortarBuff(j,2+k,m,1:nv)
Nt = 0.0_prec
uT = 0.0_prec
do i = 1,N+1
r = this%solution%interp%mortarR(i,j,k)
Nt = Nt-c*r*this%mortarNormal%boundary(i,sB,eB-offset,1:2)
uT = uT+r*this%solution%boundary(i,sB,eB-offset,1:nv)
enddo
nS = sqrt(Nt(1)*Nt(1)+Nt(2)*Nt(2))
nh = Nt/nS
fR = this%riemannflux2d(uS,uT,dsdx0,nh)
fp = this%twopointflux2d(uS,uT)
d(1:nv,j,k) = (fR-(fp(:,1)*nh(1)+fp(:,2)*nh(2)))*nS
enddo
enddo
do i = 1,N+1
uB = this%solution%boundary(i,sB,eB-offset,1:nv)
NB = this%mortarNormal%boundary(i,sB,eB-offset,1:2)
fs = this%twopointflux2d(uB,uB)
acc = 0.0_prec
do k = 1,2
do j = 1,N+1
r = this%solution%interp%mortarR(i,j,k)*this%solution%interp%qWeights(j)
uS = this%solution%mortarBuff(j,2+k,m,1:nv)
fp = this%twopointflux2d(uB,uS)-fs
acc = acc+r*(c*(fp(:,1)*NB(1)+fp(:,2)*NB(2))-d(1:nv,j,k))
enddo
enddo
this%flux%boundaryNormal(i,sB,eB-offset,1:nv) = fs(:,1)*NB(1)+fs(:,2)*NB(2)+ &
acc/this%solution%interp%qWeights(i)
enddo
enddo
endsubroutine ECMortarFlux_ECDGModel2D_t