Replaces the big-face boundaryNormal trace on each mortar interface with the L2 projection of the four small faces' boundaryNormal traces.
boundaryNormal holds the Riemann-solved surface-flux integrand f* . nHat * nScale (see BoundaryFlux in the DG models). Because the small faces' nScale is one quarter of the big face's and the sub-face coordinate Jacobian is 1/4, the projected big-face integrand is -4 * sum_q (P_kx x P_ky) g_q, where g_q are the small-face integrands and the sign accounts for the opposing outward normals. With this choice, the discrete surface integral of the big face equals minus the sum of the small faces' discrete surface integrals to roundoff, so the mortar interface is discretely conservative. Must be called after the model's BoundaryFlux and before the flux divergence is computed.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(MappedVector3D_t), | intent(inout) | :: | this | |||
| type(Mesh3D), | intent(inout) | :: | mesh |
subroutine MortarFluxCollect_MappedVector3D_t(this,mesh)
!! Replaces the big-face boundaryNormal trace on each mortar interface with the L2
!! projection of the four small faces' boundaryNormal traces.
!!
!! boundaryNormal holds the Riemann-solved surface-flux integrand f* . nHat * nScale
!! (see BoundaryFlux in the DG models). Because the small faces' nScale is one
!! quarter of the big face's and the sub-face coordinate Jacobian is 1/4, the
!! projected big-face integrand is -4 * sum_q (P_kx x P_ky) g_q, where g_q are the
!! small-face integrands and the sign accounts for the opposing outward normals.
!! With this choice, the discrete surface integral of the big face equals minus the
!! sum of the small faces' discrete surface integrals to roundoff, so the mortar
!! interface is discretely conservative. Must be called after the model's
!! BoundaryFlux and before the flux divergence is computed.
implicit none
class(MappedVector3D_t),intent(inout) :: this
type(Mesh3D),intent(inout) :: mesh
! Local
integer :: m,q,ivar,i,j,i2,j2
integer :: eB,eS,sS,flip
integer :: rankId,offset,N
integer,pointer :: elemtorank(:)
real(prec) :: extBuff(1:this%interp%N+1,1:this%interp%N+1)
! See https://github.com/FluidNumerics/SELF/issues/54 for the reason behind
! this pointer alias
elemtorank => mesh%decomp%elemToRank(:)
rankId = mesh%decomp%rankId
offset = mesh%decomp%offsetElem(rankId+1)
N = this%interp%N
if(.not. allocated(this%mortarBuff)) then
allocate(this%mortarBuff(1:N+1,1:N+1,1:8,1:mesh%nMortars,1:this%nvar,1:3))
this%mortarBuff = 0.0_prec
endif
if(mesh%decomp%mpiEnabled) then
call this%MPIMortarFluxAsync(mesh)
endif
! Stage rank-local small-face integrands in the big face's coordinates
do concurrent(m=1:mesh%nMortars,ivar=1:this%nvar)
block
integer :: i,j,i2,j2,q
integer :: eS,sS,flip
do q = 1,4
eS = mesh%mortarInfo(2*q+1,m)
if(elemtorank(eS) == rankId) then
sS = mesh%mortarInfo(2*q+2,m)/10
flip = mesh%mortarInfo(2*q+2,m)-10*sS
do j = 1,N+1
do i = 1,N+1
call MortarFaceMap(i,j,N,flip,i2,j2)
this%mortarBuff(i,j,4+q,m,ivar,1) = &
this%boundaryNormal(i2,j2,sS,eS-offset,ivar)
enddo
enddo
endif
enddo
endblock
enddo
if(mesh%decomp%mpiEnabled) then
call mesh%decomp%FinalizeMPIExchangeAsync()
! Reorient small-face integrands received over MPI into the big face's coordinates
do ivar = 1,this%nvar
do m = 1,mesh%nMortars
eB = mesh%mortarInfo(1,m)
if(elemtorank(eB) == rankId) then
do q = 1,4
eS = mesh%mortarInfo(2*q+1,m)
sS = mesh%mortarInfo(2*q+2,m)/10
flip = mesh%mortarInfo(2*q+2,m)-10*sS
if(elemtorank(eS) /= rankId .and. flip /= 0) then
do j = 1,N+1
do i = 1,N+1
call MortarFaceMap(i,j,N,flip,i2,j2)
extBuff(i,j) = this%mortarBuff(i2,j2,4+q,m,ivar,1)
enddo
enddo
do j = 1,N+1
do i = 1,N+1
this%mortarBuff(i,j,4+q,m,ivar,1) = extBuff(i,j)
enddo
enddo
endif
enddo
endif
enddo
enddo
endif
! Project the small-face integrands onto the big face's trace space. The factor of
! four converts the solution-space projection (the tensor-product mortarP carries
! the 1/4 sub-face Jacobian) into the integrand-space projection; the sign accounts
! for the opposing outward normals.
do concurrent(m=1:mesh%nMortars,ivar=1:this%nvar)
block
integer :: i,j,ii,jj,kx,ky,q
integer :: eB,sB
real(prec) :: fm
real(prec) :: tmp(1:N+1,1:N+1)
real(prec) :: acc(1:N+1,1:N+1)
eB = mesh%mortarInfo(1,m)
if(elemtorank(eB) == rankId) then
sB = mesh%mortarInfo(2,m)
acc = 0.0_prec
do q = 1,4
kx = mortarQuadKx(q)
ky = mortarQuadKy(q)
do jj = 1,N+1
do i = 1,N+1
fm = 0.0_prec
do ii = 1,N+1
fm = fm+this%interp%mortarP(ii,i,kx)* &
this%mortarBuff(ii,jj,4+q,m,ivar,1)
enddo
tmp(i,jj) = fm
enddo
enddo
do j = 1,N+1
do i = 1,N+1
fm = 0.0_prec
do jj = 1,N+1
fm = fm+this%interp%mortarP(jj,j,ky)*tmp(i,jj)
enddo
acc(i,j) = acc(i,j)+fm
enddo
enddo
enddo
do j = 1,N+1
do i = 1,N+1
this%boundaryNormal(i,j,sB,eB-offset,ivar) = -4.0_prec*acc(i,j)
enddo
enddo
endif
endblock
enddo
endsubroutine MortarFluxCollect_MappedVector3D_t