! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! 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_ECDGModel2D_t !! Entropy-conserving DGSEM base class for 2-D conservation laws. !! !! Extends DGModel2D_t by replacing the standard volume flux divergence with !! the split-form (two-point) EC volume term following !! Gassner, Winters, Kopriva (2016) !! Winters, Kopriva, Gassner, Hindenlang (2021) !! !! Semidiscretization: !! du/dt = source - (1/J) * [ 2 sum_n D_{n,i} F~^1(n,i,j) !! + 2 sum_n D_{n,j} F~^2(n,i,j) !! + M^{-1} B^T f_Riemann ] !! !! where F~^r is the Jacobian-weighted contravariant EC two-point flux in the !! r-th computational direction (computed by TwoPointFluxMethod via twopointflux2d), !! and f_Riemann is a standard Riemann solver flux at element faces (computed !! by BoundaryFlux via riemannflux2d). !! !! On 2:1 nonconforming (mortar) edges the surface flux is replaced by the !! entropy-stable two-point coupling of Friedrich et al. (2018); see !! ECMortarFlux. use SELF_DGModel2D use SELF_MappedTwoPointVector_2D implicit none type,extends(DGModel2D) :: ECDGModel2D_t type(MappedTwoPointVector2D) :: twoPointFlux ! Scaled outward edge normals (nHat*nScale), staged through the mortar exchange so ! that the rank of a small edge has the big edge's normals (mortar meshes only) type(MappedScalar2D) :: mortarNormal logical :: mortarNormalReady = .false. contains procedure :: Init => Init_ECDGModel2D_t procedure :: Free => Free_ECDGModel2D_t procedure :: Regrid => Regrid_ECDGModel2D_t procedure :: twopointflux2d => twopointflux2d_ECDGModel2D_t procedure :: TwoPointFluxMethod => TwoPointFluxMethod_ECDGModel2D_t procedure :: CalculateTendency => CalculateTendency_ECDGModel2D_t procedure :: ECMortarFlux => ECMortarFlux_ECDGModel2D_t endtype ECDGModel2D_t contains subroutine Init_ECDGModel2D_t(this,mesh,geometry) implicit none class(ECDGModel2D_t),intent(out) :: this type(Mesh2D),intent(in),target :: mesh type(SEMQuad),intent(in),target :: geometry ! Initialise all parent fields (solution, flux, source, ...) #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%nMortars > 0) then print*,__FILE__,' : Error : EC-DGSEM models on GPU do not support mortar edges yet.' error stop 1 endif #endif call Init_DGModel2D_t(this,mesh,geometry) ! Additional two-point flux field call this%twoPointFlux%Init(geometry%x%interp,this%nvar,mesh%nElem) call this%twoPointFlux%AssociateGeometry(geometry) if(mesh%nMortars > 0) then call this%mortarNormal%Init(geometry%x%interp,2,mesh%nElem) call this%mortarNormal%AssociateGeometry(geometry) this%mortarNormalReady = .true. endif endsubroutine Init_ECDGModel2D_t subroutine Free_ECDGModel2D_t(this) implicit none class(ECDGModel2D_t),intent(inout) :: this call this%twoPointFlux%DissociateGeometry() call this%twoPointFlux%Free() if(this%mortarNormalReady) then call this%mortarNormal%DissociateGeometry() call this%mortarNormal%Free() this%mortarNormalReady = .false. endif call Free_DGModel2D_t(this) endsubroutine Free_ECDGModel2D_t subroutine Regrid_ECDGModel2D_t(this,mesh,geometry) !! Rebind the EC model to a new mesh/geometry pair (AMR regrid). The parent Regrid !! resizes the DG fields and rebuilds the boundary conditions; this adds the mesh-sized !! state of the split form, which the parent does not know about: the two-point flux, and !! the mortar-normal staging field, which exists only while the mesh has mortars. The !! backend's Regrid runs, so a GPU backend's device-buffer regrid is included. !! The same mortar restrictions as Init apply to the new mesh. implicit none class(ECDGModel2D_t),intent(inout) :: this type(Mesh2D),intent(in),target :: mesh type(SEMQuad),intent(in),target :: geometry #ifdef ENABLE_GPU if(mesh%nMortars > 0) then print*,__FILE__,' : Error : EC-DGSEM models on GPU do not support mortar edges yet.' error stop 1 endif #endif call this%twoPointFlux%DissociateGeometry() call this%twoPointFlux%Free() if(this%mortarNormalReady) then call this%mortarNormal%DissociateGeometry() call this%mortarNormal%Free() this%mortarNormalReady = .false. endif ! The backend's Regrid is called by its specific name with the polymorphic this. A call ! through the parent component (this%DGModel2D%Regrid) is not polymorphic, so inside it ! AdditionalFree, AdditionalInit and SetMetadata would resolve to the DGModel2D defaults ! and the concrete model's hooks (its boundary conditions, for one) would not run. #ifdef ENABLE_GPU call Regrid_DGModel2D(this,mesh,geometry) #else call Regrid_DGModel2D_t(this,mesh,geometry) #endif call this%twoPointFlux%Init(geometry%x%interp,this%nvar,mesh%nElem) call this%twoPointFlux%AssociateGeometry(geometry) if(mesh%nMortars > 0) then call this%mortarNormal%Init(geometry%x%interp,2,mesh%nElem) call this%mortarNormal%AssociateGeometry(geometry) this%mortarNormalReady = .true. endif endsubroutine Regrid_ECDGModel2D_t pure function twopointflux2d_ECDGModel2D_t(this,sL,sR) result(flux) !! Entropy-conserving two-point flux function. !! !! Returns the physical-space two-point EC flux for the state pair (sL, sR). !! flux(ivar, d) is the d-th physical component (d=1: x, d=2: y) for variable ivar. !! !! Override this procedure in a concrete model to supply an EC flux. !! The stub returns zero (no flux), which is safe but trivial. implicit none class(ECDGModel2D_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:2) flux = 0.0_prec if(.false.) flux(1,1) = flux(1,1)+sL(1)+sR(1) ! suppress unused-argument warning endfunction twopointflux2d_ECDGModel2D_t subroutine TwoPointFluxMethod_ECDGModel2D_t(this) !! Computes pre-projected SCALAR contravariant two-point fluxes for all !! node pairs and stores them in twoPointFlux%interior(n,i,j,iel,ivar,r). !! !! Following Trixi.jl for curved meshes, each computational direction r !! uses the correct partner node AND the correct averaged metric Ja^r: !! !! r=1: pair (i,j)-(nn,j), Ja^1_avg = 0.5*(Ja^1(i,j) + Ja^1(nn,j)) !! interior(nn,i,j,...,1) = sum_d Ja^1_avg(d) * F_EC_d(sL, sR) !! !! r=2: pair (i,j)-(i,nn), Ja^2_avg = 0.5*(Ja^2(i,j) + Ja^2(i,nn)) !! interior(nn,i,j,...,2) = sum_d Ja^2_avg(d) * F_EC_d(sL, sR) !! !! The result is a SCALAR per (nn,i,j,iel,ivar,r) - no cross-contamination !! between directions. implicit none class(ECDGModel2D_t),intent(inout) :: this ! Local integer :: nn,i,j,d,iEl,iVar real(prec) :: sL(1:this%nvar),sR(1:this%nvar) real(prec) :: Fphys(1:this%nvar,1:2) real(prec) :: Fc do concurrent(nn=1:this%solution%N+1,i=1:this%solution%N+1, & j=1:this%solution%N+1,iEl=1:this%mesh%nElem) sL = this%solution%interior(i,j,iEl,1:this%nvar) ! xi^1: pair (i,j)-(nn,j), project onto avg(Ja^1) sR = this%solution%interior(nn,j,iEl,1:this%nvar) Fphys = this%twopointflux2d(sL,sR) do iVar = 1,this%nvar Fc = 0.0_prec do d = 1,2 Fc = Fc+0.5_prec*( & this%geometry%dsdx%interior(i,j,iEl,1,d,1)+ & this%geometry%dsdx%interior(nn,j,iEl,1,d,1))* & Fphys(iVar,d) enddo this%twoPointFlux%interior(nn,i,j,iEl,iVar,1) = Fc enddo ! xi^2: pair (i,j)-(i,nn), project onto avg(Ja^2) sR = this%solution%interior(i,nn,iEl,1:this%nvar) Fphys = this%twopointflux2d(sL,sR) do iVar = 1,this%nvar Fc = 0.0_prec do d = 1,2 Fc = Fc+0.5_prec*( & this%geometry%dsdx%interior(i,j,iEl,1,d,2)+ & this%geometry%dsdx%interior(i,nn,iEl,1,d,2))* & Fphys(iVar,d) enddo this%twoPointFlux%interior(nn,i,j,iEl,iVar,2) = Fc enddo enddo endsubroutine TwoPointFluxMethod_ECDGModel2D_t 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 subroutine CalculateTendency_ECDGModel2D_t(this) !! Computes du/dt = source - EC-DG flux divergence. !! !! Volume term : MappedDivergence of twoPointFlux (EC split-form sum, /J) !! Surface term: (1/J) * M^{-1} B^T f_Riemann (weak-form Riemann flux) implicit none class(ECDGModel2D_t),intent(inout) :: this ! Local integer :: i,j,iEl,iVar real(prec) :: bM1,bM2,qwi,qwj,jac call this%solution%BoundaryInterp() call this%solution%SideExchange(this%mesh) if(this%mesh%nMortars > 0) then call this%solution%MortarExchange(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() ! fills flux%boundaryNormal with Riemann fluxes if(this%mesh%nMortars > 0) then call this%ECMortarFlux() ! entropy-stable coupling on mortar edges endif ! EC volume divergence (includes J weighting via MappedDivergence) call this%TwoPointFluxMethod() call this%twoPointFlux%UpdateDevice() call this%twoPointFlux%MappedDivergence(this%fluxDivergence%interior) ! Add weak-form DG surface term: (1/J) * M^{-1} B^T f_Riemann ! Boundary side ordering: 1=South, 2=East, 3=North, 4=West do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar) bM1 = this%solution%interp%bMatrix(i,2) ! right/north boundary value at node i bM2 = this%solution%interp%bMatrix(i,1) ! left/south boundary value at node i qwi = this%solution%interp%qWeights(i) qwj = this%solution%interp%qWeights(j) jac = this%geometry%J%interior(i,j,iEl,1) ! East (side 2) and West (side 4) contributions - vary in xi^1 (i-direction) this%fluxDivergence%interior(i,j,iEl,iVar) = & this%fluxDivergence%interior(i,j,iEl,iVar)+ & (bM1*this%flux%boundaryNormal(j,2,iEl,iVar)+ & bM2*this%flux%boundaryNormal(j,4,iEl,iVar))/(qwi*jac)+ & (this%solution%interp%bMatrix(j,2)*this%flux%boundaryNormal(i,3,iEl,iVar)+ & this%solution%interp%bMatrix(j,1)*this%flux%boundaryNormal(i,1,iEl,iVar))/(qwj*jac) enddo do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar) this%dSdt%interior(i,j,iEl,iVar) = & this%source%interior(i,j,iEl,iVar)- & this%fluxDivergence%interior(i,j,iEl,iVar) enddo endsubroutine CalculateTendency_ECDGModel2D_t endmodule SELF_ECDGModel2D_t