! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! 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_MappedScalar_2D_t use SELF_Constants use SELF_Lagrange use SELF_Scalar_2D use SELF_Tensor_2D use SELF_Mesh_2D use SELF_Geometry_2D use SELF_DomainDecomposition use FEQParse use iso_c_binding implicit none type,extends(Scalar2D),public :: MappedScalar2D_t logical :: geometry_associated = .false. type(SEMQuad),pointer :: geometry => null() ! Mortar exchange work array, allocated on first use for meshes with 2:1 ! nonconforming interfaces. mortarBuff(i,slot,m,ivar) holds edge traces in the ! big side's edge orientation for mortar m : ! slot 1 - big-side trace, used to compute sub-edge 1's external state ! slot 2 - big-side trace, used to compute sub-edge 2's external state ! slot 3 - small-side trace on sub-edge 1 ! slot 4 - small-side trace on sub-edge 2 ! Slots 1 and 2 carry the same data when the big element is rank-local; they are ! distinct so that MPI receives of the big-side trace on the small sides' ranks ! never alias. real(prec),allocatable,dimension(:,:,:,:) :: mortarBuff contains procedure,public :: AssociateGeometry => AssociateGeometry_MappedScalar2D_t procedure,public :: DissociateGeometry => DissociateGeometry_MappedScalar2D_t procedure,public :: SideExchange => SideExchange_MappedScalar2D_t procedure,public :: MortarExchange => MortarExchange_MappedScalar2D_t procedure,private :: MPIMortarExchangeAsync => MPIMortarExchangeAsync_MappedScalar2D_t generic,public :: MappedGradient => MappedGradient_MappedScalar2D_t procedure,private :: MappedGradient_MappedScalar2D_t generic,public :: MappedDGGradient => MappedDGGradient_MappedScalar2D_t procedure,private :: MappedDGGradient_MappedScalar2D_t procedure,private :: MPIExchangeAsync => MPIExchangeAsync_MappedScalar2D_t procedure,private :: ApplyFlip => ApplyFlip_MappedScalar2D_t procedure,public :: SetInteriorFromEquation => SetInteriorFromEquation_MappedScalar2D_t endtype MappedScalar2D_t contains subroutine AssociateGeometry_MappedScalar2D_t(this,geometry) implicit none class(MappedScalar2D_t),intent(inout) :: this type(SEMQuad),target,intent(in) :: geometry if(.not. associated(this%geometry)) then this%geometry => geometry this%geometry_associated = .true. endif endsubroutine AssociateGeometry_MappedScalar2D_t subroutine DissociateGeometry_MappedScalar2D_t(this) implicit none class(MappedScalar2D_t),intent(inout) :: this if(associated(this%geometry)) then this%geometry => null() this%geometry_associated = .false. endif endsubroutine DissociateGeometry_MappedScalar2D_t subroutine SetInteriorFromEquation_MappedScalar2D_t(this,geometry,time) !! Sets the this % interior attribute using the eqn attribute, !! geometry (for physical positions), and provided simulation time. implicit none class(MappedScalar2D_t),intent(inout) :: this type(SEMQuad),intent(in) :: geometry real(prec),intent(in) :: time ! Local integer :: i,j,iEl,iVar real(prec) :: x real(prec) :: y do iVar = 1,this%nVar do iEl = 1,this%nElem do j = 1,this%interp%N+1 do i = 1,this%interp%N+1 ! Get the mesh positions x = geometry%x%interior(i,j,iEl,1,1) y = geometry%x%interior(i,j,iEl,1,2) this%interior(i,j,iEl,iVar) = & this%eqn(iVar)%Evaluate((/x,y,0.0_prec,time/)) enddo enddo enddo enddo endsubroutine SetInteriorFromEquation_MappedScalar2D_t subroutine MPIExchangeAsync_MappedScalar2D_t(this,mesh) implicit none class(MappedScalar2D_t),intent(inout) :: this type(Mesh2D),intent(inout) :: mesh ! Local integer :: e1,s1,e2,s2,ivar integer :: globalSideId,r2,tag integer :: iError integer :: msgCount msgCount = 0 do ivar = 1,this%nvar do e1 = 1,this%nElem do s1 = 1,4 e2 = mesh%sideInfo(3,s1,e1) ! Neighbor Element if(e2 > 0) then r2 = mesh%decomp%elemToRank(e2) ! Neighbor Rank if(r2 /= mesh%decomp%rankId) then s2 = mesh%sideInfo(4,s1,e1)/10 globalSideId = abs(mesh%sideInfo(2,s1,e1)) ! create unique tag for each side and each variable tag = globalsideid+mesh%nUniqueSides*(ivar-1) msgCount = msgCount+1 call MPI_IRECV(this%extBoundary(:,s1,e1,ivar), & (this%interp%N+1), & mesh%decomp%mpiPrec, & r2,tag, & mesh%decomp%mpiComm, & mesh%decomp%requests(msgCount),iError) msgCount = msgCount+1 call MPI_ISEND(this%boundary(:,s1,e1,ivar), & (this%interp%N+1), & mesh%decomp%mpiPrec, & r2,tag, & mesh%decomp%mpiComm, & mesh%decomp%requests(msgCount),iError) endif endif enddo enddo enddo mesh%decomp%msgCount = msgCount endsubroutine MPIExchangeAsync_MappedScalar2D_t subroutine ApplyFlip_MappedScalar2D_t(this,mesh) ! Apply side flips to sides where MPI exchanges took place. implicit none class(MappedScalar2D_t),intent(inout) :: this type(Mesh2D),intent(in) :: mesh ! Local integer :: e1,s1,e2,s2 integer :: i,i2 integer :: r2,flip,ivar real(prec) :: extBuff(1:this%interp%N+1) do ivar = 1,this%nvar do e1 = 1,this%nElem do s1 = 1,4 e2 = mesh%sideInfo(3,s1,e1) ! Neighbor Element (global id) if(e2 > 0) then ! Interior Element r2 = mesh%decomp%elemToRank(e2) ! Neighbor Rank if(r2 /= mesh%decomp%rankId) then s2 = mesh%sideInfo(4,s1,e1)/10 flip = mesh%sideInfo(4,s1,e1)-s2*10 ! Need to update extBoundary with flip applied if(flip == 1) then do i = 1,this%interp%N+1 i2 = this%interp%N+2-i extBuff(i) = this%extBoundary(i2,s1,e1,ivar) enddo do i = 1,this%interp%N+1 this%extBoundary(i,s1,e1,ivar) = extBuff(i) enddo endif endif endif enddo enddo enddo endsubroutine ApplyFlip_MappedScalar2D_t subroutine SideExchange_MappedScalar2D_t(this,mesh) implicit none class(MappedScalar2D_t),intent(inout) :: this type(Mesh2D),intent(inout) :: mesh ! Local integer :: e1,e2,s1,s2,e2Global integer :: flip integer :: i1,i2,ivar integer :: r2 integer :: rankId,offset,N integer,pointer :: elemtorank(:) ! This mapping is needed to resolve a build error with ! amdflang that appears to be caused by referencing ! the elemToRank attribute within the do concurrent ! https://github.com/FluidNumerics/SELF/issues/54 elemtorank => mesh%decomp%elemToRank(:) rankId = mesh%decomp%rankId offset = mesh%decomp%offsetElem(rankId+1) N = this%interp%N if(mesh%decomp%mpiEnabled) then call this%MPIExchangeAsync(mesh) endif do concurrent(s1=1:4,e1=1:mesh%nElem,ivar=1:this%nvar) e2Global = mesh%sideInfo(3,s1,e1) e2 = e2Global-offset s2 = mesh%sideInfo(4,s1,e1)/10 flip = mesh%sideInfo(4,s1,e1)-s2*10 if(e2Global > 0) then r2 = elemToRank(e2Global) if(r2 == rankId) then if(flip == 0) then do i1 = 1,N+1 this%extBoundary(i1,s1,e1,ivar) = & this%boundary(i1,s2,e2,ivar) enddo elseif(flip == 1) then do i1 = 1,N+1 i2 = N+2-i1 this%extBoundary(i1,s1,e1,ivar) = & this%boundary(i2,s2,e2,ivar) enddo endif endif endif enddo if(mesh%decomp%mpiEnabled) then call mesh%decomp%FinalizeMPIExchangeAsync() ! Apply side flips for data exchanged with MPI call this%ApplyFlip(mesh) endif endsubroutine SideExchange_MappedScalar2D_t subroutine MPIMortarExchangeAsync_MappedScalar2D_t(this,mesh) !! Posts the point-to-point messages required for mortar interfaces whose big and !! small elements reside on different ranks. The big-side rank sends its edge trace !! to each remote small-side rank and receives each remote small-side trace; tags !! are built from the sub-edge global side ids, following the conforming SideExchange !! tag convention. implicit none class(MappedScalar2D_t),intent(inout) :: this type(Mesh2D),intent(inout) :: mesh ! Local integer :: m,k,ivar integer :: eB,sB,rB,eS,sS,rS integer :: globalSideId,tag integer :: offset integer :: iError integer :: msgCount msgCount = 0 offset = mesh%decomp%offsetElem(mesh%decomp%rankId+1) do ivar = 1,this%nvar do m = 1,mesh%nMortars eB = mesh%mortarInfo(1,m) sB = mesh%mortarInfo(2,m) rB = mesh%decomp%elemToRank(eB) do k = 1,2 eS = mesh%mortarInfo(2*k+1,m) sS = mesh%mortarInfo(2*k+2,m)/10 rS = mesh%decomp%elemToRank(eS) globalSideId = mesh%mortarInfo(6+k,m) tag = globalSideId+mesh%nUniqueSides*(ivar-1) if(rB == mesh%decomp%rankId .and. rS /= mesh%decomp%rankId) then msgCount = msgCount+1 call MPI_IRECV(this%mortarBuff(:,2+k,m,ivar), & (this%interp%N+1), & mesh%decomp%mpiPrec, & rS,tag, & mesh%decomp%mpiComm, & mesh%decomp%requests(msgCount),iError) msgCount = msgCount+1 call MPI_ISEND(this%boundary(:,sB,eB-offset,ivar), & (this%interp%N+1), & mesh%decomp%mpiPrec, & rS,tag, & mesh%decomp%mpiComm, & mesh%decomp%requests(msgCount),iError) elseif(rS == mesh%decomp%rankId .and. rB /= mesh%decomp%rankId) then msgCount = msgCount+1 call MPI_IRECV(this%mortarBuff(:,k,m,ivar), & (this%interp%N+1), & mesh%decomp%mpiPrec, & rB,tag, & mesh%decomp%mpiComm, & mesh%decomp%requests(msgCount),iError) msgCount = msgCount+1 call MPI_ISEND(this%boundary(:,sS,eS-offset,ivar), & (this%interp%N+1), & mesh%decomp%mpiPrec, & rB,tag, & mesh%decomp%mpiComm, & mesh%decomp%requests(msgCount),iError) endif enddo enddo enddo mesh%decomp%msgCount = msgCount endsubroutine MPIMortarExchangeAsync_MappedScalar2D_t subroutine MortarExchange_MappedScalar2D_t(this,mesh) !! Fills the extBoundary attribute on all sides that participate in a 2:1 !! nonconforming (mortar) interface. !! !! Each small side receives the big side's trace restricted (exactly) to its !! sub-edge with the interp's mortarR matrices. The big side receives the L2 !! projection of the two small-side traces onto its trace space with the interp's !! mortarP matrices. Traces are staged in mortarBuff in the big side's edge !! orientation; the flip recorded in mesh%mortarInfo translates between the small !! sides' orientations and the big side's orientation. Must be called after !! BoundaryInterp and may be called before, after, or without SideExchange (the !! side sets touched are disjoint from conforming sides). implicit none class(MappedScalar2D_t),intent(inout) :: this type(Mesh2D),intent(inout) :: mesh ! Local integer :: m,k,ivar,i,ii integer :: eB,sB,eS,sS,flip integer :: rankId,offset,N integer,pointer :: elemtorank(:) real(prec) :: fm real(prec) :: extBuff(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:4,1:mesh%nMortars,1:this%nvar)) this%mortarBuff = 0.0_prec endif if(mesh%decomp%mpiEnabled) then call this%MPIMortarExchangeAsync(mesh) endif ! Stage rank-local traces in the big side's edge orientation do concurrent(m=1:mesh%nMortars,ivar=1:this%nvar) eB = mesh%mortarInfo(1,m) if(elemtorank(eB) == rankId) then sB = mesh%mortarInfo(2,m) do i = 1,N+1 this%mortarBuff(i,1,m,ivar) = this%boundary(i,sB,eB-offset,ivar) this%mortarBuff(i,2,m,ivar) = this%boundary(i,sB,eB-offset,ivar) enddo endif do k = 1,2 eS = mesh%mortarInfo(2*k+1,m) if(elemtorank(eS) == rankId) then sS = mesh%mortarInfo(2*k+2,m)/10 flip = mesh%mortarInfo(2*k+2,m)-10*sS if(flip == 0) then do i = 1,N+1 this%mortarBuff(i,2+k,m,ivar) = this%boundary(i,sS,eS-offset,ivar) enddo else do i = 1,N+1 this%mortarBuff(i,2+k,m,ivar) = this%boundary(N+2-i,sS,eS-offset,ivar) enddo endif endif enddo enddo if(mesh%decomp%mpiEnabled) then call mesh%decomp%FinalizeMPIExchangeAsync() ! Reorient small-side traces received over MPI into the big side's orientation do ivar = 1,this%nvar do m = 1,mesh%nMortars eB = mesh%mortarInfo(1,m) if(elemtorank(eB) == rankId) then do k = 1,2 eS = mesh%mortarInfo(2*k+1,m) sS = mesh%mortarInfo(2*k+2,m)/10 flip = mesh%mortarInfo(2*k+2,m)-10*sS if(elemtorank(eS) /= rankId .and. flip == 1) then do i = 1,N+1 extBuff(i) = this%mortarBuff(N+2-i,2+k,m,ivar) enddo do i = 1,N+1 this%mortarBuff(i,2+k,m,ivar) = extBuff(i) enddo endif enddo endif enddo enddo endif ! Compute external states : ! small sides get the restricted big-side trace (exact), ! the big side gets the L2 projection of the small-side traces do concurrent(m=1:mesh%nMortars,ivar=1:this%nvar) do k = 1,2 eS = mesh%mortarInfo(2*k+1,m) if(elemtorank(eS) == rankId) then sS = mesh%mortarInfo(2*k+2,m)/10 flip = mesh%mortarInfo(2*k+2,m)-10*sS do i = 1,N+1 fm = 0.0_prec do ii = 1,N+1 fm = fm+this%interp%mortarR(ii,i,k)*this%mortarBuff(ii,k,m,ivar) enddo if(flip == 0) then this%extBoundary(i,sS,eS-offset,ivar) = fm else this%extBoundary(N+2-i,sS,eS-offset,ivar) = fm endif enddo endif enddo eB = mesh%mortarInfo(1,m) if(elemtorank(eB) == rankId) then sB = mesh%mortarInfo(2,m) do i = 1,N+1 fm = 0.0_prec do k = 1,2 do ii = 1,N+1 fm = fm+this%interp%mortarP(ii,i,k)*this%mortarBuff(ii,2+k,m,ivar) enddo enddo this%extBoundary(i,sB,eB-offset,ivar) = fm enddo endif enddo endsubroutine MortarExchange_MappedScalar2D_t subroutine MappedGradient_MappedScalar2D_t(this,df) !! Calculates the gradient of a function using the strong form of the gradient !! in mapped coordinates. implicit none class(MappedScalar2D_t),intent(in) :: this real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%nelem,1:this%nvar,1:2) ! Local integer :: iEl,iVar,i,j,ii,idir real(prec) :: dfdx,ja do concurrent(i=1:this%N+1,j=1:this%N+1,iel=1:this%nElem,ivar=1:this%nVar,idir=1:2) dfdx = 0.0_prec do ii = 1,this%N+1 ! dsdx(j,i) is contravariant vector i, component j ja = this%geometry%dsdx%interior(ii,j,iel,1,idir,1) dfdx = dfdx+this%interp%dMatrix(ii,i)*this%interior(ii,j,iel,ivar)*ja enddo df(i,j,iel,ivar,idir) = dfdx enddo do concurrent(i=1:this%N+1,j=1:this%N+1,iel=1:this%nElem,ivar=1:this%nVar,idir=1:2) dfdx = 0.0_prec do ii = 1,this%N+1 ja = this%geometry%dsdx%interior(i,ii,iel,1,idir,2) dfdx = dfdx+this%interp%dMatrix(ii,j)*this%interior(i,ii,iel,ivar)*ja enddo df(i,j,iel,ivar,idir) = (df(i,j,iel,ivar,idir)+dfdx)/this%geometry%J%interior(i,j,iEl,1) enddo endsubroutine MappedGradient_MappedScalar2D_t subroutine MappedDGGradient_MappedScalar2D_t(this,df) !! implicit none class(MappedScalar2D_t),intent(in) :: this real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%nelem,1:this%nvar,1:2) ! Local integer :: iEl,iVar,i,j,ii,idir real(prec) :: dfdx,dfdxb,ja,bfl,bfr do concurrent(i=1:this%N+1,j=1:this%N+1,iel=1:this%nElem,ivar=1:this%nVar,idir=1:2) dfdx = 0.0_prec do ii = 1,this%N+1 ja = this%geometry%dsdx%interior(ii,j,iel,1,idir,1) dfdx = dfdx+this%interp%dgMatrix(ii,i)*this%interior(ii,j,iel,ivar)*ja enddo bfl = this%avgboundary(j,4,iel,ivar)*this%geometry%dsdx%boundary(j,4,iel,1,idir,1) ! west bfr = this%avgboundary(j,2,iel,ivar)*this%geometry%dsdx%boundary(j,2,iel,1,idir,1) ! east dfdxb = (this%interp%bMatrix(i,1)*bfl+this%interp%bMatrix(i,2)*bfr)/this%interp%qweights(i) df(i,j,iel,ivar,idir) = dfdx+dfdxb enddo do concurrent(i=1:this%N+1,j=1:this%N+1,iel=1:this%nElem,ivar=1:this%nVar,idir=1:2) dfdx = 0.0_prec do ii = 1,this%N+1 ja = this%geometry%dsdx%interior(i,ii,iel,1,idir,2) dfdx = dfdx+this%interp%dgMatrix(ii,j)*this%interior(i,ii,iel,ivar)*ja enddo bfl = this%avgboundary(i,1,iel,ivar)*this%geometry%dsdx%boundary(i,1,iel,1,idir,2) ! south bfr = this%avgboundary(i,3,iel,ivar)*this%geometry%dsdx%boundary(i,3,iel,1,idir,2) ! north dfdxb = (this%interp%bMatrix(j,1)*bfl+this%interp%bMatrix(j,2)*bfr)/this%interp%qweights(j) df(i,j,iel,ivar,idir) = (df(i,j,iel,ivar,idir)+dfdx+dfdxb)/this%geometry%J%interior(i,j,iEl,1) enddo endsubroutine MappedDGGradient_MappedScalar2D_t endmodule SELF_MappedScalar_2D_t