Entropy-stable surface fluxes on anisotropic (2-subface) mortar faces, the 3-D counterpart of ECMortarFlux_ECDGModel2D_t (Friedrich et al., J. Sci. Comput. 77, 2018). Overwrites flux%boundaryNormal on every anisotropic mortar face, after BoundaryFlux.
A big face split along trace coordinate s couples only along s: big node (i',l) and node (i,l) of half h share the line index l of the other coordinate. With R_ii' = mortarR(i',i,h) along s, the 1-D weights w along s, the big face's scaled outward normal N = nHat*nScale and c = 1/2 (the half to big-face area ratio),
half h, node (i,l) : F = sum_i' R_ii' f#(u_h(i,l), u(i',l)) . (-c N(i',l)) big node (i',l) : F = (1/w_i') sum_h sum_i R_ii' w_i [ c f#(u(i',l), u_h(i,l)) . N(i',l) - d_h(i,l) ]
evaluated, as in 2-D, as the local flux plus a correction (the small face's own normal in its local term), so a constant state is preserved exactly; and with the dissipation d_h of ECMortarFlux_ECDGModel2D_t evaluated against the restricted big trace. Conservation and the entropy estimate follow as in 2-D; the weights along l cancel because R is the identity along l. Traces come from solution%anisoMortarBuff (slot h: big trace, 2+h: small trace in the big face's coordinates) and the big faces' normals from mortarNormal%anisoMortarBuff.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(ECDGModel3D_t), | intent(inout) | :: | this |
subroutine ECMortarFlux_ECDGModel3D_t(this)
!! Entropy-stable surface fluxes on anisotropic (2-subface) mortar faces, the 3-D
!! counterpart of ECMortarFlux_ECDGModel2D_t (Friedrich et al., J. Sci. Comput. 77,
!! 2018). Overwrites flux%boundaryNormal on every anisotropic mortar face, after
!! BoundaryFlux.
!!
!! A big face split along trace coordinate s couples only along s: big node (i',l)
!! and node (i,l) of half h share the line index l of the other coordinate. With
!! R_ii' = mortarR(i',i,h) along s, the 1-D weights w along s, the big face's
!! scaled outward normal N = nHat*nScale and c = 1/2 (the half to big-face area
!! ratio),
!!
!! half h, node (i,l) : F = sum_i' R_ii' f#(u_h(i,l), u(i',l)) . (-c N(i',l))
!! big node (i',l) : F = (1/w_i') sum_h sum_i R_ii' w_i [ c f#(u(i',l), u_h(i,l)) . N(i',l)
!! - d_h(i,l) ]
!!
!! evaluated, as in 2-D, as the local flux plus a correction (the small face's own
!! normal in its local term), so a constant state is preserved exactly; and
!! with the dissipation d_h of ECMortarFlux_ECDGModel2D_t evaluated against the
!! restricted big trace. Conservation and the entropy estimate follow as in 2-D; the
!! weights along l cancel because R is the identity along l. Traces come from
!! solution%anisoMortarBuff (slot h: big trace, 2+h: small trace in the big face's
!! coordinates) and the big faces' normals from mortarNormal%anisoMortarBuff.
implicit none
class(ECDGModel3D_t),intent(inout) :: this
! Local
real(prec),parameter :: c = 0.5_prec
integer :: N,nv,m,h,a,ap,l,i,j,ib,jb,i2,j2,iel,s,eB,sB,eS,sS,flip,split,offset,rankId,iv
real(prec) :: uB(1:this%nvar),uS(1:this%nvar),uT(1:this%nvar)
real(prec) :: fp(1:this%nvar,1:3),fs(1:this%nvar,1:3),fR(1:this%nvar),acc(1:this%nvar)
real(prec) :: dsdx0(1:this%nvar,1:3)
real(prec) :: NB(1:3),Nt(1:3),nS,nh(1:3),r
real(prec) :: d(1:this%nvar,1:this%solution%N+1,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 faces.
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
do iel = 1,this%mesh%nElem
do s = 1,6
do j = 1,N+1
do i = 1,N+1
this%mortarNormal%boundary(i,j,s,iel,1:3) = this%geometry%nHat%boundary(i,j,s,iel,1,1:3)* &
this%geometry%nScale%boundary(i,j,s,iel,1)
enddo
enddo
enddo
enddo
call this%mortarNormal%AnisoMortarExchange(this%mesh)
do m = 1,this%mesh%nAnisoMortars
split = this%mesh%anisoMortarInfo(9,m)
! ---- Small faces ----
do h = 1,2
eS = this%mesh%anisoMortarInfo(2*h+1,m)
if(this%mesh%decomp%elemToRank(eS) /= rankId) cycle
sS = this%mesh%anisoMortarInfo(2*h+2,m)/10
flip = this%mesh%anisoMortarInfo(2*h+2,m)-10*sS
do l = 1,N+1
do a = 1,N+1
call SplitIndex(split,a,l,i,j)
uS = this%solution%anisoMortarBuff(i,j,2+h,m,1:nv)
call MortarFaceMap(i,j,N,flip,i2,j2)
fs = this%twopointflux3d(uS,uS)
NB = this%mortarNormal%boundary(i2,j2,sS,eS-offset,1:3)
acc = fs(:,1)*NB(1)+fs(:,2)*NB(2)+fs(:,3)*NB(3)
Nt = 0.0_prec
uT = 0.0_prec
do ap = 1,N+1
call SplitIndex(split,ap,l,ib,jb)
r = this%solution%interp%mortarR(ap,a,h)
uB = this%solution%anisoMortarBuff(ib,jb,h,m,1:nv)
NB = this%mortarNormal%anisoMortarBuff(ib,jb,h,m,1:3)
fp = this%twopointflux3d(uS,uB)-fs
acc = acc-c*r*(fp(:,1)*NB(1)+fp(:,2)*NB(2)+fp(:,3)*NB(3))
Nt = Nt-c*r*NB
uT = uT+r*uB
enddo
nS = sqrt(Nt(1)*Nt(1)+Nt(2)*Nt(2)+Nt(3)*Nt(3))
nh = Nt/nS
fR = this%riemannflux3d(uS,uT,dsdx0,nh)
fp = this%twopointflux3d(uS,uT)
do iv = 1,nv
this%flux%boundaryNormal(i2,j2,sS,eS-offset,iv) = &
acc(iv)+(fR(iv)-(fp(iv,1)*nh(1)+fp(iv,2)*nh(2)+fp(iv,3)*nh(3)))*nS
enddo
enddo
enddo
enddo
! ---- Big face ----
eB = this%mesh%anisoMortarInfo(1,m)
if(this%mesh%decomp%elemToRank(eB) /= rankId) cycle
sB = this%mesh%anisoMortarInfo(2,m)
do h = 1,2
do l = 1,N+1
do a = 1,N+1
call SplitIndex(split,a,l,i,j)
uS = this%solution%anisoMortarBuff(i,j,2+h,m,1:nv)
Nt = 0.0_prec
uT = 0.0_prec
do ap = 1,N+1
call SplitIndex(split,ap,l,ib,jb)
r = this%solution%interp%mortarR(ap,a,h)
Nt = Nt-c*r*this%mortarNormal%boundary(ib,jb,sB,eB-offset,1:3)
uT = uT+r*this%solution%boundary(ib,jb,sB,eB-offset,1:nv)
enddo
nS = sqrt(Nt(1)*Nt(1)+Nt(2)*Nt(2)+Nt(3)*Nt(3))
nh = Nt/nS
fR = this%riemannflux3d(uS,uT,dsdx0,nh)
fp = this%twopointflux3d(uS,uT)
d(1:nv,i,j,h) = (fR-(fp(:,1)*nh(1)+fp(:,2)*nh(2)+fp(:,3)*nh(3)))*nS
enddo
enddo
enddo
do l = 1,N+1
do ap = 1,N+1
call SplitIndex(split,ap,l,ib,jb)
uB = this%solution%boundary(ib,jb,sB,eB-offset,1:nv)
NB = this%mortarNormal%boundary(ib,jb,sB,eB-offset,1:3)
fs = this%twopointflux3d(uB,uB)
acc = 0.0_prec
do h = 1,2
do a = 1,N+1
call SplitIndex(split,a,l,i,j)
r = this%solution%interp%mortarR(ap,a,h)*this%solution%interp%qWeights(a)
uS = this%solution%anisoMortarBuff(i,j,2+h,m,1:nv)
fp = this%twopointflux3d(uB,uS)-fs
acc = acc+r*(c*(fp(:,1)*NB(1)+fp(:,2)*NB(2)+fp(:,3)*NB(3))-d(1:nv,i,j,h))
enddo
enddo
this%flux%boundaryNormal(ib,jb,sB,eB-offset,1:nv) = fs(:,1)*NB(1)+fs(:,2)*NB(2)+fs(:,3)*NB(3)+ &
acc/this%solution%interp%qWeights(ap)
enddo
enddo
enddo
contains
pure subroutine SplitIndex(split,a,l,i,j)
!! Face indices (i,j) of the node at position a along the split coordinate and l
!! along the other.
implicit none
integer,intent(in) :: split,a,l
integer,intent(out) :: i,j
if(split == 1) then
i = a
j = l
else
i = l
j = a
endif
endsubroutine SplitIndex
endsubroutine ECMortarFlux_ECDGModel3D_t