! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! 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_SetInterface_3D_t !! Interface map and trace exchange between element sets. !! !! An interface face is a face of the global mesh whose two elements belong to different !! element sets (SELF_ElementSetMesh_3D). In the sub-mesh of each set it is tagged !! SELF_BC_SET_INTERFACE, and its exterior state must come from the other set. For every !! interface face owned by this rank, this module receives the partner element's face !! trace, !! !! ghost(i,j,v,f) = b_p(i2,j2,s2,e2,v), v = 1..nvar(p), !! !! where p is the partner set, (e2,s2) the partner element and side, b_p the partner's !! boundary trace, and (i2,j2) the receiver's flip mapping of (i,j). The mapping is the !! one SideExchange applies for a face with flip = mod(sideInfo(4),10), so a ghost value !! equals the exterior state the unsplit mesh would have had on that face. !! !! Faces. Local interface faces are numbered f = 1..nFaces, grouped by receiving set: !! the faces of set k are setFaceStart(k) .. setFaceStart(k+1)-1, in global element !! order. Each face both sends its own trace to the partner's rank and receives the !! partner's trace. !! !! Messages. Faces whose partner is on another rank are grouped by partner rank. Within !! each group the faces are ordered by global side id, which both ranks know, and each !! face carries Np*Np*nvar values of the sending set (Np = N+1). The receiver knows the !! partner set's nvar, so both ends agree on the layout without negotiation. The exchange !! uses persistent point-to-point requests with a constant tag on a duplicate of the !! global mesh communicator, so it cannot match any halo message of the sub-meshes. !! Faces whose partner is on the same rank are copied directly. !! !! Each exchange runs in two parts. GatherOutbox copies every local interface face's own !! trace, unflipped, into outbox(:,:,v,f); a GPU backend overrides it to gather on the !! device and copy the outbox to the host. ExchangeOutbox then sends, receives and flips !! on the host, so both backends fill ghost with the same values. !! !! Collectives. BuildMap duplicates the communicator and Free releases it; both are !! collective over the global mesh communicator. Exchange posts point-to-point messages !! only, and only between ranks that share interface faces. use SELF_Constants use SELF_Mesh_3D use iso_c_binding use mpi implicit none type SetTraceRef !! Face trace of one set: trace(1:Np,1:Np,1:6,1:nElem,1:nvar), for example !! solution%boundary of the set's model. real(prec),pointer,contiguous :: trace(:,:,:,:,:) => null() !! Device copy of the same trace (GPU backend); left null on the host path. type(c_ptr) :: trace_gpu = c_null_ptr endtype SetTraceRef type SetInterface3D_t integer :: nSets = 0 integer :: N = 0 integer,allocatable :: nvarOfSet(:) integer :: nvarMax = 0 !! Local interface faces, f = 1..nFaces integer :: nFaces = 0 integer,allocatable :: setFaceStart(:) !! (1:nSets+1) integer,allocatable :: faceSet(:) !! receiving set k integer,allocatable :: faceElem(:) !! element in set k's sub-mesh integer,allocatable :: faceSide(:) !! local side 1..6 integer,allocatable :: faceFlip(:) !! mod(sideInfo(4),10) of the receiving side integer,allocatable :: facePartnerSet(:) !! partner set p integer,allocatable :: facePartnerRank(:) !! rank owning the partner element integer,allocatable :: facePartnerElem(:) !! partner element in set p's sub-mesh; 0 if remote integer,allocatable :: facePartnerSide(:) !! partner local side 1..6 integer,allocatable :: faceSideId(:) !! |global side id| integer,allocatable :: facePartnerFace(:) !! local face index of the partner; 0 if remote !! Own and partner traces, (1:Np,1:Np,1:nvarMax,1:nFaces). Pointers so a GPU backend !! can copy them to and from device memory with c_loc. real(prec),pointer,contiguous :: outbox(:,:,:,:) => null() real(prec),pointer,contiguous :: ghost(:,:,:,:) => null() !! Neighbour ranks and their faces integer :: nNbr = 0 integer,allocatable :: nbrRank(:) integer,allocatable :: nbrFaceStart(:) !! (1:nNbr+1) integer,allocatable :: nbrFace(:) !! local face index, ordered by neighbour then side id integer,allocatable :: sendStart(:) !! (1:nNbr+1), offsets into sendBuf integer,allocatable :: recvStart(:) !! (1:nNbr+1), offsets into recvBuf real(prec),allocatable :: sendBuf(:) real(prec),allocatable :: recvBuf(:) integer,allocatable :: requests(:) integer :: comm = MPI_COMM_NULL integer :: mpiPrec integer :: rankId = 0 contains procedure :: BuildMap => BuildMap_SetInterface3D_t procedure :: Exchange => Exchange_SetInterface3D_t procedure :: GatherOutbox => GatherOutbox_SetInterface3D_t procedure :: ExchangeOutbox => ExchangeOutbox_SetInterface3D_t procedure :: Free => Free_SetInterface3D_t endtype SetInterface3D_t integer,parameter,private :: SET_INTERFACE_TAG = 0 contains pure subroutine FlipIndex(i,j,flip,N,i2,j2) !! Index (i2,j2) on the sending face that supplies node (i,j) of the receiving face, !! for the receiving side's flip in 0..7. Identical to the mapping in !! SideExchange_MappedScalar3D_t and ApplyFlip_MappedScalar3D_t. implicit none integer,intent(in) :: i,j,flip,N integer,intent(out) :: i2,j2 select case(flip) case(1) i2 = N+2-i; j2 = j case(2) i2 = N+2-i; j2 = N+2-j case(3) i2 = i; j2 = N+2-j case(4) i2 = j; j2 = i case(5) i2 = N+2-j; j2 = i case(6) i2 = N+2-j; j2 = N+2-i case(7) i2 = j; j2 = N+2-i case default i2 = i; j2 = j endselect endsubroutine FlipIndex subroutine HeapSortByKey(keys,vals) !! Sort keys(1:n) ascending in place, permuting vals(1:n) alongside. Heap sort: !! O(n log n), no extra storage. Keys are global side ids, distinct within one !! neighbour's face list, so stability does not matter. implicit none integer,intent(inout) :: keys(1:) integer,intent(inout) :: vals(1:) ! Local integer :: n,i,last n = size(keys) do i = n/2,1,-1 call SiftDown(i,n) enddo do last = n,2,-1 call SwapEntries(1,last) call SiftDown(1,last-1) enddo contains subroutine SiftDown(first,last) implicit none integer,intent(in) :: first,last integer :: root,child root = first do while(2*root <= last) child = 2*root if(child < last) then if(keys(child+1) > keys(child)) child = child+1 endif if(keys(root) >= keys(child)) return call SwapEntries(root,child) root = child enddo endsubroutine SiftDown subroutine SwapEntries(i1,i2) implicit none integer,intent(in) :: i1,i2 integer :: t t = keys(i1); keys(i1) = keys(i2); keys(i2) = t t = vals(i1); vals(i1) = vals(i2); vals(i2) = t endsubroutine SwapEntries endsubroutine HeapSortByKey subroutine BuildMap_SetInterface3D_t(this,globalMesh,setOfElem,setGlobalId, & setOffsetThisRank,nvarOfSet,N) !! Build the interface map from the global mesh and the set numbering, without !! communication beyond duplicating the communicator. Collective over !! globalMesh%decomp%mpiComm. !! !! Input !! - globalMesh : the decomposed global mesh !! - setOfElem(1:nGlobalElem) : set id of every global element !! - setGlobalId(1:nGlobalElem) : set-global id of every global element !! - setOffsetThisRank(1:nSets) : number of set-k elements on lower ranks, so the !! sub-mesh element of global element e is setGlobalId(e)-setOffsetThisRank(s(e)) !! - nvarOfSet(1:nSets) : number of variables in each set's trace !! - N : polynomial degree; faces hold (N+1)*(N+1) nodes implicit none class(SetInterface3D_t),intent(inout) :: this type(Mesh3D),intent(in) :: globalMesh integer,intent(in) :: setOfElem(1:) integer,intent(in) :: setGlobalId(1:) integer,intent(in) :: setOffsetThisRank(1:) integer,intent(in) :: nvarOfSet(1:) integer,intent(in) :: N ! Local integer :: nSets,iEl,s1,eg,e2g,k,p,f,c,r,a,b,nn,ierror,np2,offset integer,allocatable :: counts(:),order(:),keys(:),faceOfSide(:,:) logical,allocatable :: isNbr(:) nSets = size(nvarOfSet) this%nSets = nSets this%N = N np2 = (N+1)*(N+1) this%nvarOfSet = nvarOfSet this%nvarMax = maxval(nvarOfSet) this%rankId = globalMesh%decomp%rankId this%mpiPrec = globalMesh%decomp%mpiPrec call mpi_comm_dup(globalMesh%decomp%mpiComm,this%comm,ierror) ! Count local interface faces per receiving set allocate(counts(1:nSets)) counts = 0 do iEl = 1,globalMesh%nElem eg = globalMesh%decomp%offsetElem(this%rankId+1)+iEl do s1 = 1,6 e2g = globalMesh%sideInfo(3,s1,iEl) if(e2g <= 0) cycle if(setOfElem(e2g) == setOfElem(eg)) cycle counts(setOfElem(eg)) = counts(setOfElem(eg))+1 enddo enddo allocate(this%setFaceStart(1:nSets+1)) this%setFaceStart(1) = 1 do k = 1,nSets this%setFaceStart(k+1) = this%setFaceStart(k)+counts(k) enddo this%nFaces = this%setFaceStart(nSets+1)-1 nn = this%nFaces allocate(this%faceSet(1:nn),this%faceElem(1:nn),this%faceSide(1:nn), & this%faceFlip(1:nn),this%facePartnerSet(1:nn),this%facePartnerRank(1:nn), & this%facePartnerElem(1:nn),this%facePartnerSide(1:nn),this%faceSideId(1:nn), & this%facePartnerFace(1:nn)) allocate(this%outbox(1:N+1,1:N+1,1:this%nvarMax,1:nn)) allocate(this%ghost(1:N+1,1:N+1,1:this%nvarMax,1:nn)) this%outbox = 0.0_prec this%ghost = 0.0_prec offset = globalMesh%decomp%offsetElem(this%rankId+1) allocate(faceOfSide(1:6,1:globalMesh%nElem)) faceOfSide = 0 ! Fill the faces, grouped by set, in global element order within each set counts = 0 do iEl = 1,globalMesh%nElem eg = globalMesh%decomp%offsetElem(this%rankId+1)+iEl do s1 = 1,6 e2g = globalMesh%sideInfo(3,s1,iEl) if(e2g <= 0) cycle k = setOfElem(eg) p = setOfElem(e2g) if(p == k) cycle f = this%setFaceStart(k)+counts(k) counts(k) = counts(k)+1 faceOfSide(s1,iEl) = f this%faceSet(f) = k this%faceElem(f) = setGlobalId(eg)-setOffsetThisRank(k) this%faceSide(f) = s1 this%facePartnerSide(f) = globalMesh%sideInfo(4,s1,iEl)/10 this%faceFlip(f) = globalMesh%sideInfo(4,s1,iEl)-10*this%facePartnerSide(f) this%facePartnerSet(f) = p this%facePartnerRank(f) = globalMesh%decomp%elemToRank(e2g) this%faceSideId(f) = abs(globalMesh%sideInfo(2,s1,iEl)) if(this%facePartnerRank(f) == this%rankId) then this%facePartnerElem(f) = setGlobalId(e2g)-setOffsetThisRank(p) else this%facePartnerElem(f) = 0 endif enddo enddo deallocate(counts) ! Same-rank partner faces, through the partner's global element and side do f = 1,nn this%facePartnerFace(f) = 0 enddo do iEl = 1,globalMesh%nElem do s1 = 1,6 f = faceOfSide(s1,iEl) if(f == 0) cycle if(this%facePartnerRank(f) /= this%rankId) cycle e2g = globalMesh%sideInfo(3,s1,iEl) this%facePartnerFace(f) = faceOfSide(this%facePartnerSide(f),e2g-offset) enddo enddo deallocate(faceOfSide) ! Neighbour ranks, in increasing rank order allocate(isNbr(0:globalMesh%decomp%nRanks-1)) isNbr = .false. do f = 1,nn if(this%facePartnerRank(f) /= this%rankId) isNbr(this%facePartnerRank(f)) = .true. enddo this%nNbr = count(isNbr) allocate(this%nbrRank(1:this%nNbr)) c = 0 do r = 0,globalMesh%decomp%nRanks-1 if(isNbr(r)) then c = c+1 this%nbrRank(c) = r endif enddo deallocate(isNbr) ! Faces per neighbour, ordered by global side id allocate(this%nbrFaceStart(1:this%nNbr+1)) allocate(this%nbrFace(1:count(this%facePartnerRank /= this%rankId))) this%nbrFaceStart(1) = 1 do c = 1,this%nNbr a = this%nbrFaceStart(c) b = a-1 do f = 1,nn if(this%facePartnerRank(f) == this%nbrRank(c)) then b = b+1 this%nbrFace(b) = f endif enddo this%nbrFaceStart(c+1) = b+1 if(b > a) then allocate(order(1:b-a+1),keys(1:b-a+1)) order = this%nbrFace(a:b) keys = this%faceSideId(order) call HeapSortByKey(keys,order) this%nbrFace(a:b) = order deallocate(order,keys) endif enddo ! Buffer offsets: a face sends nvar of its own set and receives nvar of the partner's allocate(this%sendStart(1:this%nNbr+1),this%recvStart(1:this%nNbr+1)) this%sendStart(1) = 1 this%recvStart(1) = 1 do c = 1,this%nNbr this%sendStart(c+1) = this%sendStart(c) this%recvStart(c+1) = this%recvStart(c) do b = this%nbrFaceStart(c),this%nbrFaceStart(c+1)-1 f = this%nbrFace(b) this%sendStart(c+1) = this%sendStart(c+1)+np2*nvarOfSet(this%faceSet(f)) this%recvStart(c+1) = this%recvStart(c+1)+np2*nvarOfSet(this%facePartnerSet(f)) enddo enddo allocate(this%sendBuf(1:this%sendStart(this%nNbr+1)-1)) allocate(this%recvBuf(1:this%recvStart(this%nNbr+1)-1)) ! Persistent requests: receives first, then sends allocate(this%requests(1:2*this%nNbr)) do c = 1,this%nNbr call mpi_recv_init(this%recvBuf(this%recvStart(c)), & this%recvStart(c+1)-this%recvStart(c),this%mpiPrec, & this%nbrRank(c),SET_INTERFACE_TAG,this%comm, & this%requests(c),ierror) call mpi_send_init(this%sendBuf(this%sendStart(c)), & this%sendStart(c+1)-this%sendStart(c),this%mpiPrec, & this%nbrRank(c),SET_INTERFACE_TAG,this%comm, & this%requests(this%nNbr+c),ierror) enddo endsubroutine BuildMap_SetInterface3D_t subroutine Exchange_SetInterface3D_t(this,traces) !! Fill ghost(:,:,1:nvar(p),f) for every local interface face f with the partner's !! face trace, flipped into the receiving face's orientation. !! !! Input !! - traces(1:nSets) : each set's current face trace; read only for sets that own !! local interface faces implicit none class(SetInterface3D_t),intent(inout) :: this type(SetTraceRef),intent(in) :: traces(1:) call this%GatherOutbox(traces) call this%ExchangeOutbox() endsubroutine Exchange_SetInterface3D_t subroutine GatherOutbox_SetInterface3D_t(this,traces) !! outbox(:,:,v,f) = trace of set k on face f, v = 1..nvar(k), unflipped, for every !! local interface face f of every set k. Host traces. implicit none class(SetInterface3D_t),intent(inout) :: this type(SetTraceRef),intent(in) :: traces(1:) ! Local integer :: f,k,nv do f = 1,this%nFaces k = this%faceSet(f) nv = this%nvarOfSet(k) this%outbox(:,:,1:nv,f) = traces(k)%trace(:,:,this%faceSide(f),this%faceElem(f),1:nv) enddo endsubroutine GatherOutbox_SetInterface3D_t subroutine ExchangeOutbox_SetInterface3D_t(this) !! Send each remote face's outbox to the partner's rank, and fill !! !! ghost(i,j,v,f) = partner outbox(i2,j2,v), v = 1..nvar(p), !! !! with (i2,j2) the receiving side's flip mapping of (i,j). Same-rank partners are read !! from the local outbox. Point-to-point only, on this%comm. implicit none class(SetInterface3D_t),intent(inout) :: this ! Local integer :: c,b,f,g,k,p,v,i,j,i2,j2,pos,ierror integer :: np1 np1 = this%N+1 ! Pack this rank's side of every remote face do c = 1,this%nNbr pos = this%sendStart(c) do b = this%nbrFaceStart(c),this%nbrFaceStart(c+1)-1 f = this%nbrFace(b) k = this%faceSet(f) do v = 1,this%nvarOfSet(k) do j = 1,np1 do i = 1,np1 this%sendBuf(pos) = this%outbox(i,j,v,f) pos = pos+1 enddo enddo enddo enddo enddo if(this%nNbr > 0) call mpi_startall(2*this%nNbr,this%requests,ierror) ! Same-rank partners do f = 1,this%nFaces g = this%facePartnerFace(f) if(g == 0) cycle p = this%facePartnerSet(f) do v = 1,this%nvarOfSet(p) do j = 1,np1 do i = 1,np1 call FlipIndex(i,j,this%faceFlip(f),this%N,i2,j2) this%ghost(i,j,v,f) = this%outbox(i2,j2,v,g) enddo enddo enddo enddo if(this%nNbr == 0) return call mpi_waitall(2*this%nNbr,this%requests,MPI_STATUSES_IGNORE,ierror) ! Unpack remote partners' traces; the sender packed (i,j) unflipped do c = 1,this%nNbr pos = this%recvStart(c) do b = this%nbrFaceStart(c),this%nbrFaceStart(c+1)-1 f = this%nbrFace(b) p = this%facePartnerSet(f) do v = 1,this%nvarOfSet(p) do j = 1,np1 do i = 1,np1 call FlipIndex(i,j,this%faceFlip(f),this%N,i2,j2) this%ghost(i,j,v,f) = this%recvBuf(pos+(i2-1)+np1*(j2-1)) enddo enddo pos = pos+np1*np1 enddo enddo enddo endsubroutine ExchangeOutbox_SetInterface3D_t subroutine Free_SetInterface3D_t(this) !! Release the persistent requests, buffers and the duplicate communicator. !! Collective over the communicator the map was built on. implicit none class(SetInterface3D_t),intent(inout) :: this ! Local integer :: c,ierror if(allocated(this%requests)) then do c = 1,size(this%requests) call mpi_request_free(this%requests(c),ierror) enddo deallocate(this%requests) endif if(allocated(this%nvarOfSet)) deallocate(this%nvarOfSet) if(allocated(this%setFaceStart)) deallocate(this%setFaceStart) if(allocated(this%faceSet)) deallocate(this%faceSet) if(allocated(this%faceElem)) deallocate(this%faceElem) if(allocated(this%faceSide)) deallocate(this%faceSide) if(allocated(this%faceFlip)) deallocate(this%faceFlip) if(allocated(this%facePartnerSet)) deallocate(this%facePartnerSet) if(allocated(this%facePartnerRank)) deallocate(this%facePartnerRank) if(allocated(this%facePartnerElem)) deallocate(this%facePartnerElem) if(allocated(this%facePartnerSide)) deallocate(this%facePartnerSide) if(allocated(this%faceSideId)) deallocate(this%faceSideId) if(allocated(this%facePartnerFace)) deallocate(this%facePartnerFace) if(associated(this%outbox)) deallocate(this%outbox) if(associated(this%ghost)) deallocate(this%ghost) if(allocated(this%nbrRank)) deallocate(this%nbrRank) if(allocated(this%nbrFaceStart)) deallocate(this%nbrFaceStart) if(allocated(this%nbrFace)) deallocate(this%nbrFace) if(allocated(this%sendStart)) deallocate(this%sendStart) if(allocated(this%recvStart)) deallocate(this%recvStart) if(allocated(this%sendBuf)) deallocate(this%sendBuf) if(allocated(this%recvBuf)) deallocate(this%recvBuf) if(this%comm /= MPI_COMM_NULL) call mpi_comm_free(this%comm,ierror) this%nSets = 0 this%nFaces = 0 this%nNbr = 0 endsubroutine Free_SetInterface3D_t endmodule SELF_SetInterface_3D_t