! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! Maintainers : support@fluidnumerics.com ! Official Repository : https://github.com/FluidNumerics/self/ ! ! Copyright © 2024 Fluid Numerics LLC ! ! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in ! the documentation and/or other materials provided with the distribution. ! ! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from ! this software without specific prior written permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT ! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT ! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY ! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF ! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! module SELF_ECDGModel3D_t !! Entropy-conserving DGSEM base class for 3-D conservation laws. !! See SELF_ECDGModel2D_t for design documentation. !! !! On anisotropic (2-subface) mortar faces the surface flux is replaced by the !! entropy-stable two-point coupling of Friedrich et al. (2018); see ECMortarFlux. !! Meshes with 2x2 mortar faces are not supported yet. use SELF_DGModel3D use SELF_MappedTwoPointVector_3D implicit none type,extends(DGModel3D) :: ECDGModel3D_t type(MappedTwoPointVector3D) :: twoPointFlux ! Scaled outward face normals (nHat*nScale), staged through the anisotropic mortar ! exchange so that the rank of a small face has the big face's normals type(MappedScalar3D) :: mortarNormal contains procedure :: Init => Init_ECDGModel3D_t procedure :: Free => Free_ECDGModel3D_t procedure :: twopointflux3d => twopointflux3d_ECDGModel3D_t procedure :: TwoPointFluxMethod => TwoPointFluxMethod_ECDGModel3D_t procedure :: CalculateTendency => CalculateTendency_ECDGModel3D_t procedure :: ECMortarFlux => ECMortarFlux_ECDGModel3D_t endtype ECDGModel3D_t contains subroutine Init_ECDGModel3D_t(this,mesh,geometry) implicit none class(ECDGModel3D_t),intent(out) :: this type(Mesh3D),intent(in),target :: mesh type(SEMHex),intent(in),target :: geometry ! The entropy-conserving split form requires entropy-stable mortar operators at ! nonconforming interfaces. They are implemented for anisotropic (2-subface) faces ! only; the plain L2 projection of a 2x2 mortar would break the provable entropy ! estimate. if(mesh%nMortars > 0) then print*,__FILE__,' : Error : EC-DGSEM models do not support 2x2 nonconforming (mortar) faces.' error stop 1 endif #ifdef ENABLE_GPU ! The entropy-stable mortar coupling (ECMortarFlux) runs on the host and uses the ! host two-point flux; the GPU tendency has no device version of it yet. if(mesh%nAnisoMortars > 0) then print*,__FILE__,' : Error : EC-DGSEM models on GPU do not support anisotropic mortar faces yet.' error stop 1 endif #endif call Init_DGModel3D_t(this,mesh,geometry) call this%twoPointFlux%Init(geometry%x%interp,this%nvar,mesh%nElem) call this%twoPointFlux%AssociateGeometry(geometry) if(mesh%nAnisoMortars > 0) then call this%mortarNormal%Init(geometry%x%interp,3,mesh%nElem) call this%mortarNormal%AssociateGeometry(geometry) endif endsubroutine Init_ECDGModel3D_t subroutine Free_ECDGModel3D_t(this) implicit none class(ECDGModel3D_t),intent(inout) :: this call this%twoPointFlux%DissociateGeometry() call this%twoPointFlux%Free() if(this%mesh%nAnisoMortars > 0) then call this%mortarNormal%DissociateGeometry() call this%mortarNormal%Free() endif call Free_DGModel3D_t(this) endsubroutine Free_ECDGModel3D_t pure function twopointflux3d_ECDGModel3D_t(this,sL,sR) result(flux) !! Entropy-conserving two-point flux function (3-D). !! flux(ivar, d) is the d-th physical component (d=1:x, d=2:y, d=3:z). !! Override in concrete models. Stub returns zero. implicit none class(ECDGModel3D_t),intent(in) :: this real(prec),intent(in) :: sL(1:this%nvar) real(prec),intent(in) :: sR(1:this%nvar) real(prec) :: flux(1:this%nvar,1:3) flux = 0.0_prec if(.false.) flux(1,1) = flux(1,1)+sL(1)+sR(1) ! suppress unused-argument warning endfunction twopointflux3d_ECDGModel3D_t subroutine TwoPointFluxMethod_ECDGModel3D_t(this) !! Computes pre-projected SCALAR contravariant two-point fluxes for all !! node pairs, following Trixi.jl for curved meshes. Each direction r !! uses the correct partner AND the correct averaged metric Ja^r. implicit none class(ECDGModel3D_t),intent(inout) :: this ! Local integer :: nn,i,j,k,d,iEl,iVar real(prec) :: sL(1:this%nvar),sR(1:this%nvar) real(prec) :: Fphys(1:this%nvar,1:3) real(prec) :: Fc do concurrent(nn=1:this%solution%N+1,i=1:this%solution%N+1, & j=1:this%solution%N+1,k=1:this%solution%N+1, & iEl=1:this%mesh%nElem) sL = this%solution%interior(i,j,k,iEl,1:this%nvar) ! xi^1: pair (i,j,k)-(nn,j,k), project onto avg(Ja^1) sR = this%solution%interior(nn,j,k,iEl,1:this%nvar) Fphys = this%twopointflux3d(sL,sR) do iVar = 1,this%nvar Fc = 0.0_prec do d = 1,3 Fc = Fc+0.5_prec*( & this%geometry%dsdx%interior(i,j,k,iEl,1,d,1)+ & this%geometry%dsdx%interior(nn,j,k,iEl,1,d,1))* & Fphys(iVar,d) enddo this%twoPointFlux%interior(nn,i,j,k,iEl,iVar,1) = Fc enddo ! xi^2: pair (i,j,k)-(i,nn,k), project onto avg(Ja^2) sR = this%solution%interior(i,nn,k,iEl,1:this%nvar) Fphys = this%twopointflux3d(sL,sR) do iVar = 1,this%nvar Fc = 0.0_prec do d = 1,3 Fc = Fc+0.5_prec*( & this%geometry%dsdx%interior(i,j,k,iEl,1,d,2)+ & this%geometry%dsdx%interior(i,nn,k,iEl,1,d,2))* & Fphys(iVar,d) enddo this%twoPointFlux%interior(nn,i,j,k,iEl,iVar,2) = Fc enddo ! xi^3: pair (i,j,k)-(i,j,nn), project onto avg(Ja^3) sR = this%solution%interior(i,j,nn,iEl,1:this%nvar) Fphys = this%twopointflux3d(sL,sR) do iVar = 1,this%nvar Fc = 0.0_prec do d = 1,3 Fc = Fc+0.5_prec*( & this%geometry%dsdx%interior(i,j,k,iEl,1,d,3)+ & this%geometry%dsdx%interior(i,j,nn,iEl,1,d,3))* & Fphys(iVar,d) enddo this%twoPointFlux%interior(nn,i,j,k,iEl,iVar,3) = Fc enddo enddo endsubroutine TwoPointFluxMethod_ECDGModel3D_t 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 subroutine CalculateTendency_ECDGModel3D_t(this) implicit none class(ECDGModel3D_t),intent(inout) :: this ! Local integer :: i,j,k,iEl,iVar real(prec) :: bMi1,bMi2,bMj1,bMj2,bMk1,bMk2,qwi,qwj,qwk,jac call this%solution%BoundaryInterp() call this%solution%SideExchange(this%mesh) if(this%mesh%nAnisoMortars > 0) then call this%solution%AnisoMortarExchange(this%mesh) endif call this%PreTendencyHook() call this%SetBoundaryCondition() if(this%gradient_enabled) then call this%CalculateSolutionGradient() call this%SetGradientBoundaryCondition() call this%solutionGradient%AverageSides() endif call this%SourceMethod() call this%BoundaryFlux() if(this%mesh%nAnisoMortars > 0) then call this%ECMortarFlux() ! entropy-stable coupling on anisotropic mortar faces endif call this%TwoPointFluxMethod() call this%twoPointFlux%UpdateDevice() call this%twoPointFlux%MappedDivergence(this%fluxDivergence%interior) ! Add (1/J) * M^{-1} B^T f_Riemann ! 3D boundary side ordering: 1=Bottom, 2=South, 3=East, 4=North, 5=West, 6=Top do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & k=1:this%solution%N+1,iEl=1:this%mesh%nElem, & iVar=1:this%solution%nVar) bMi1 = this%solution%interp%bMatrix(i,2) bMi2 = this%solution%interp%bMatrix(i,1) bMj1 = this%solution%interp%bMatrix(j,2) bMj2 = this%solution%interp%bMatrix(j,1) bMk1 = this%solution%interp%bMatrix(k,2) bMk2 = this%solution%interp%bMatrix(k,1) qwi = this%solution%interp%qWeights(i) qwj = this%solution%interp%qWeights(j) qwk = this%solution%interp%qWeights(k) jac = this%geometry%J%interior(i,j,k,iEl,1) this%fluxDivergence%interior(i,j,k,iEl,iVar) = & this%fluxDivergence%interior(i,j,k,iEl,iVar)+ & (bMk1*this%flux%boundaryNormal(i,j,6,iEl,iVar)+ & ! top bMk2*this%flux%boundaryNormal(i,j,1,iEl,iVar))/(qwk*jac)+ & ! bottom (bMi1*this%flux%boundaryNormal(j,k,3,iEl,iVar)+ & ! east bMi2*this%flux%boundaryNormal(j,k,5,iEl,iVar))/(qwi*jac)+ & ! west (bMj1*this%flux%boundaryNormal(i,k,4,iEl,iVar)+ & ! north bMj2*this%flux%boundaryNormal(i,k,2,iEl,iVar))/(qwj*jac) ! south enddo do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & k=1:this%solution%N+1,iEl=1:this%mesh%nElem, & iVar=1:this%solution%nVar) this%dSdt%interior(i,j,k,iEl,iVar) = & this%source%interior(i,j,k,iEl,iVar)- & this%fluxDivergence%interior(i,j,k,iEl,iVar) enddo endsubroutine CalculateTendency_ECDGModel3D_t endmodule SELF_ECDGModel3D_t