! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! 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_AMRController_Extruded !! Adaptive-mesh-refinement controller for models on EXTRUDED meshes: a 2-D quadtree on the !! surface, with every column extruded through the same nz layers. One call to Adapt performs !! one adaptation epoch between time steps, carrying a 3-D model on the extruded mesh and, !! optionally, a 2-D model on the surface mesh it was extruded from: !! !! 1. estimate - the 3-D modal-decay indicator flags every element of the 3-D model, and !! each column takes the strongest request: refine if any of its layers asks, !! coarsen only if all of them agree, keep otherwise; !! 2. cap, halo - as in AMRController2D: refine flags at maxLevel are demoted, and refine !! flags spread to face-neighbouring columns for nHalo passes; !! 3. mutate - AdaptFromFlags + Balance2to1 update the surface forest; an unchanged leaf !! set ends the epoch with both models untouched; !! 4. rebuild - EmitMesh produces the new surface mesh (with 2:1 mortar edges), the !! caller's levels procedure fills the level heights z_k at its geometry !! nodes, and ExtrudedMesh builds the new 3-D mesh (whose lateral faces !! between columns of different levels become anisotropic mortars); !! 5. transfer - the 2-D transfer plan of the forest drives both models: the 2-D model !! through its own ApplyTransferPlan, and the 3-D model through !! ApplyExtrudedTransferPlan, which applies the plan to every layer of a !! column (the surface operator tensored with the identity in zeta). !! !! The levels procedure is called with each new surface mesh because the levels depend on the !! position of the surface nodes: with terrain-following levels, refinement samples the !! bathymetry at the new nodes (ExtrudedSigmaLevels and ExtrudedHybridLevels are the usual !! bodies). It must return levels that are continuous across lateral faces, as ExtrudedMesh !! requires. !! !! MPI: as in AMRController2D, the forest is rank-replicated, the column flags are allgathered !! once per epoch, and EmitMesh re-decomposes the new surface leaves into contiguous ranges. !! Columns never straddle ranks, so the 3-D decomposition is the surface one expanded by nz and !! the 2-D migration windows (PlanWindows) become 3-D windows by the same factor. The 3-D model !! migrates point-to-point with MigrateOldWindow_DGModel3D_t; the 2-D model with its own !! MigrateOldWindow. !! !! Backends: the 3-D model always stages and migrates through the portable host routines of !! DGModel3D_t (called by name), because the extruded apply reads host data; on a GPU build !! that means one device-to-host and one host-to-device copy of the 3-D field per adaptation !! epoch. The 2-D model uses its backend's own transfer. Nothing is added to the time-stepping !! loop. !! !! Ownership: the caller's base surface mesh, extruded mesh and geometries are never freed !! here; the base surface mesh must outlive the controller (it supplies the boundary-condition !! metadata and the communicator for every emitted mesh). Meshes and geometries emitted by an !! epoch are owned by the controller and freed when they are replaced or by Free. use iso_fortran_env,only:int64 use SELF_Constants use SELF_Lagrange use SELF_Mesh_2D use SELF_Mesh_3D use SELF_Geometry_2D use SELF_Geometry_3D use SELF_DGModel2D_t use SELF_DGModel3D_t,only:DGModel3D_t,StageSolutionForTransfer_DGModel3D_t, & MigrateOldWindow_DGModel3D_t,ApplyExtrudedTransferPlan_DGModel3D_t use SELF_QuadTreeMesh_2D use SELF_AdaptiveMesh_2D use SELF_RefinementIndicator_3D use SELF_TransferPlan_2D use SELF_AMRController_2D,only:PlanWindows,AllgatherPerElemInts,InitForestFromDecomposedMesh use mpi implicit none abstract interface subroutine SELF_ExtrudedLevels(surf,zLevels) !! Fill zLevels(1:nGeo+1,1:nGeo+1,0:nz,1:surf%nElem), the level heights z_0 < ... < !! z_nz = 0 at the geometry nodes of every rank-local surface element, for ExtrudedMesh. use SELF_Constants,only:prec use SELF_Mesh_2D,only:Mesh2D implicit none type(Mesh2D),intent(in) :: surf real(prec),intent(out) :: zLevels(1:,1:,0:,1:) endsubroutine SELF_ExtrudedLevels endinterface type :: ExtrudedAMRController type(QuadTreeMesh2D) :: forest type(RefinementIndicator3D) :: indicator procedure(SELF_ExtrudedLevels),pointer,nopass :: levels => null() type(Mesh2D),pointer :: baseSurf => null() !! caller-owned; metadata source for EmitMesh type(Mesh2D),pointer :: activeSurf => null() !! surface mesh of the current columns type(Mesh3D),pointer :: activeMesh => null() !! extruded mesh the 3-D model runs on type(SEMHex),pointer :: activeGeom => null() !! its geometry type(SEMQuad),pointer :: activeGeom2 => null() !! surface geometry (only with a 2-D model) type(Lagrange),pointer :: interp => null() !! the models' solution interpolant logical :: ownsActive = .false. !! whether the active meshes/geometries were emitted by us logical :: hasModel2D = .false. !! whether Init was given a 2-D model integer :: nz = 0 !! layers per column integer :: bottomBC = 0 !! boundary condition id of the column bottoms integer :: topBC = 0 !! boundary condition id of the column tops integer :: ivar = SELF_AMR_ALLVARS !! driving 3-D variable for the indicator integer :: maxLevel = 1 !! refinement-level cap integer :: nHalo = 1 !! refine-flag halo-expansion passes (in columns) real(prec) :: refineThreshold = 0.0_prec real(prec) :: coarsenThreshold = 0.0_prec real(prec) :: relativeEnergyFloor = SELF_AMR_DEFAULT_RELFLOOR real(prec) :: significantEnergyFloor = SELF_AMR_DEFAULT_RELFLOOR real(prec),allocatable :: energyWeights(:) !! Migration windows for every rank, in surface elements (columns) and in 3-D elements integer,allocatable :: winFirst(:),winLast(:) integer,allocatable :: winFirst3(:),winLast3(:) !! Cumulative migration accounting, both models together integer(int64) :: nMigrateBytesRecv = 0 integer(int64) :: nMigrateBytesSent = 0 integer(int64) :: nMigrateElemRemote = 0 contains procedure,public :: Init => Init_ExtrudedAMRController procedure,public :: Free => Free_ExtrudedAMRController procedure,public :: Adapt => Adapt_ExtrudedAMRController procedure,public :: RecommendedTimeStep => RecommendedTimeStep_ExtrudedAMRController procedure,private :: ApplyIndicatorSettings => ApplyIndicatorSettings_ExtrudedAMRController endtype ExtrudedAMRController contains subroutine Init_ExtrudedAMRController(this,model3d,surf,levels,bottomBC,topBC, & refineThreshold,coarsenThreshold,ivar,maxLevel,nHalo, & model2d,relativeEnergyFloor,energyWeights, & significantEnergyFloor) !! Attach the controller to an initialized 3-D model on an extruded mesh, and optionally to !! a 2-D model on the surface mesh. !! !! - model3d : runs on ExtrudedMesh(surf, ...) with nz = model3d%mesh%nLayers layers. !! - surf : the surface mesh the 3-D mesh was extruded from; it becomes the forest's base !! (level 0) and must outlive the controller. !! - levels : fills the level heights for a new surface mesh (see SELF_ExtrudedLevels). !! - bottomBC, topBC : the boundary condition ids of the column bottoms and tops, as !! passed to ExtrudedMesh. !! - refineThreshold > coarsenThreshold : sigma = log10 modal-energy-ratio cut-offs of the !! 3-D refinement indicator; ivar its driving variable (or SELF_AMR_ALLVARS). !! - maxLevel >= 0 caps the refinement depth; nHalo >= 0 the refine halo, in columns. !! - model2d (optional) : runs on surf, and is regridded and transferred with model3d. !! - relativeEnergyFloor, energyWeights, significantEnergyFloor : the indicator's !! amplitude gate, as for AMRController3D. implicit none class(ExtrudedAMRController),intent(out) :: this class(DGModel3D_t),intent(in) :: model3d type(Mesh2D),intent(in),target :: surf procedure(SELF_ExtrudedLevels) :: levels integer,intent(in) :: bottomBC integer,intent(in) :: topBC real(prec),intent(in) :: refineThreshold real(prec),intent(in) :: coarsenThreshold integer,intent(in) :: ivar integer,intent(in) :: maxLevel integer,intent(in) :: nHalo class(DGModel2D_t),intent(in),optional :: model2d real(prec),intent(in),optional :: relativeEnergyFloor real(prec),intent(in),optional :: energyWeights(:) real(prec),intent(in),optional :: significantEnergyFloor if(.not. associated(model3d%mesh)) then print*,__FILE__,':',__LINE__,' : Error : ExtrudedAMRController%Init requires an initialized 3-D model.' stop 1 endif if(model3d%mesh%nLayers < 1) then print*,__FILE__,':',__LINE__,' : Error : ExtrudedAMRController%Init requires a model on an extruded mesh.' stop 1 endif if(model3d%mesh%nElem /= model3d%mesh%nLayers*surf%nElem) then print*,__FILE__,':',__LINE__, & ' : Error : the 3-D model''s mesh was not extruded from the given surface mesh.' stop 1 endif if(maxLevel < 0 .or. nHalo < 0) then print*,__FILE__,':',__LINE__,' : Error : ExtrudedAMRController%Init requires maxLevel >= 0 and nHalo >= 0.' stop 1 endif if(present(model2d)) then if(.not. associated(model2d%mesh,surf)) then print*,__FILE__,':',__LINE__,' : Error : the 2-D model does not run on the given surface mesh.' stop 1 endif if(model2d%geometry%x%interp%N /= model3d%geometry%x%interp%N) then print*,__FILE__,':',__LINE__,' : Error : the 2-D and 3-D models must use the same degree.' stop 1 endif this%hasModel2D = .true. this%activeGeom2 => model2d%geometry endif this%levels => levels this%baseSurf => surf this%activeSurf => surf this%activeMesh => model3d%mesh this%activeGeom => model3d%geometry this%interp => model3d%geometry%x%interp this%ownsActive = .false. this%nz = model3d%mesh%nLayers this%bottomBC = bottomBC this%topBC = topBC this%ivar = ivar this%maxLevel = maxLevel this%nHalo = nHalo this%refineThreshold = refineThreshold this%coarsenThreshold = coarsenThreshold if(present(relativeEnergyFloor)) then this%relativeEnergyFloor = relativeEnergyFloor this%significantEnergyFloor = relativeEnergyFloor endif if(present(significantEnergyFloor)) this%significantEnergyFloor = significantEnergyFloor if(present(energyWeights)) then allocate(this%energyWeights(1:size(energyWeights))) this%energyWeights = energyWeights endif if(surf%decomp%nRanks > 1) then call InitForestFromDecomposedMesh(this%forest,surf) else call this%forest%Init(surf) endif call this%indicator%Init(this%interp,model3d%mesh%nElem,refineThreshold,coarsenThreshold) call this%ApplyIndicatorSettings() endsubroutine Init_ExtrudedAMRController subroutine ApplyIndicatorSettings_ExtrudedAMRController(this) !! Push the controller-owned amplitude-gate settings onto the indicator, after every !! indicator Init (which resets them). implicit none class(ExtrudedAMRController),intent(inout) :: this call this%indicator%SetRelativeEnergyFloor(this%relativeEnergyFloor, & significantEnergyFloor=this%significantEnergyFloor) if(allocated(this%energyWeights)) then call this%indicator%SetEnergyWeights(this%energyWeights) endif endsubroutine ApplyIndicatorSettings_ExtrudedAMRController subroutine Free_ExtrudedAMRController(this) !! Release the forest, the indicator and every mesh/geometry the controller emitted. The !! models' mesh and geometry pointers dangle afterwards if they were running on emitted !! ones, so free or stop using the models first. implicit none class(ExtrudedAMRController),intent(inout) :: this call ReleaseActive(this) call this%forest%Free() call this%indicator%Free() if(allocated(this%energyWeights)) deallocate(this%energyWeights) if(allocated(this%winFirst)) deallocate(this%winFirst,this%winLast) if(allocated(this%winFirst3)) deallocate(this%winFirst3,this%winLast3) this%baseSurf => null() this%activeSurf => null() this%activeMesh => null() this%activeGeom => null() this%activeGeom2 => null() this%interp => null() this%levels => null() endsubroutine Free_ExtrudedAMRController subroutine ReleaseActive(this) !! Free the active meshes and geometries if this controller emitted them. implicit none type(ExtrudedAMRController),intent(inout) :: this if(.not. this%ownsActive) return call this%activeGeom%Free() deallocate(this%activeGeom) call this%activeMesh%Free() deallocate(this%activeMesh) if(this%hasModel2D) then call this%activeGeom2%Free() deallocate(this%activeGeom2) endif call this%activeSurf%Free() deallocate(this%activeSurf) this%ownsActive = .false. endsubroutine ReleaseActive subroutine Adapt_ExtrudedAMRController(this,model3d,adapted,model2d) !! Perform one adaptation epoch (see the module documentation). On return adapted reports !! whether the mesh changed; when it did, both models are rebound to the new meshes with !! their solutions transferred, and the caller should re-evaluate the time step !! (RecommendedTimeStep) and any model state derived from the mesh. model2d must be given !! exactly when Init was given one. implicit none class(ExtrudedAMRController),intent(inout),target :: this class(DGModel3D_t),intent(inout) :: model3d logical,intent(out) :: adapted class(DGModel2D_t),intent(inout),optional :: model2d ! Local integer :: c,k,li,s,pass,node,nbr,ns,nf,nOld,nCol,nz,nR,changed integer :: cFirst,cLast,wF,wL,wF3,wL3,nRefine,nCoarsen integer,allocatable :: colFlag(:),flag(:),spread(:),oldLeaf(:),leafIdx(:) real(prec),allocatable :: zLevels(:,:,:,:) type(TransferPlan2D) :: plan type(Mesh2D),pointer :: newSurf type(Mesh3D),pointer :: newMesh type(SEMHex),pointer :: newGeom type(SEMQuad),pointer :: newGeom2 adapted = .false. nz = this%nz if(.not. associated(model3d%mesh,this%activeMesh)) then print*,__FILE__,':',__LINE__,' : Error : the 3-D model is not running on this controller''s active mesh.' stop 1 endif if(present(model2d) .neqv. this%hasModel2D) then print*,__FILE__,':',__LINE__,' : Error : Adapt must be given a 2-D model exactly when Init was.' stop 1 endif if(present(model2d)) then if(.not. associated(model2d%mesh,this%activeSurf)) then print*,__FILE__,':',__LINE__,' : Error : the 2-D model is not running on the active surface mesh.' stop 1 endif endif ! ---- 1. Element flags from the 3-D solution, reduced to column flags ---- if(model3d%mesh%decomp%mpiEnabled) then call this%indicator%Estimate(model3d%solution,this%ivar,comm=model3d%mesh%decomp%mpiComm) else call this%indicator%Estimate(model3d%solution,this%ivar) endif nCol = this%activeSurf%nElem allocate(colFlag(1:max(nCol,1))) do c = 1,nCol nRefine = 0 nCoarsen = 0 do k = 1,nz if(this%indicator%flag((c-1)*nz+k) == SELF_AMR_REFINE) nRefine = nRefine+1 if(this%indicator%flag((c-1)*nz+k) == SELF_AMR_COARSEN) nCoarsen = nCoarsen+1 enddo if(nRefine > 0) then colFlag(c) = SELF_AMR_REFINE elseif(nCoarsen == nz) then colFlag(c) = SELF_AMR_COARSEN else colFlag(c) = SELF_AMR_KEEP endif enddo nOld = this%forest%nLeaves allocate(flag(1:nOld)) if(this%activeSurf%decomp%nRanks > 1) then if(this%activeSurf%decomp%nElem /= nOld) then print*,__FILE__,':',__LINE__,' : Error : the active decomposition does not span the forest leaf list.' stop 1 endif call AllgatherPerElemInts(this%activeSurf%decomp,1,colFlag,flag) else flag(1:nOld) = colFlag(1:nOld) endif deallocate(colFlag) ! ---- 2. Cap and halo, as in AMRController2D ---- do li = 1,nOld if(flag(li) == SELF_AMR_REFINE .and. & this%forest%level(this%forest%leaf(li)) >= this%maxLevel) then flag(li) = SELF_AMR_KEEP endif enddo allocate(leafIdx(1:this%forest%nNodes)) leafIdx = 0 do li = 1,nOld leafIdx(this%forest%leaf(li)) = li enddo allocate(spread(1:nOld)) do pass = 1,this%nHalo spread(1:nOld) = flag(1:nOld) do li = 1,nOld if(flag(li) /= SELF_AMR_REFINE) cycle node = this%forest%leaf(li) do s = 1,4 call this%forest%FaceNeighbor(node,s,nbr,ns,nf) if(nbr == 0) cycle if(this%forest%child(1,nbr) /= 0) cycle if(this%forest%level(nbr) >= this%maxLevel) cycle spread(leafIdx(nbr)) = SELF_AMR_REFINE enddo enddo flag(1:nOld) = spread(1:nOld) enddo deallocate(spread,leafIdx) ! ---- 3. Mutate the forest; an unchanged leaf set ends the epoch ---- allocate(oldLeaf(1:nOld)) oldLeaf(1:nOld) = this%forest%leaf(1:nOld) call this%forest%AdaptFromFlags(flag) deallocate(flag) call this%forest%Balance2to1() changed = 1 if(this%forest%nLeaves == nOld) then changed = 0 do li = 1,nOld if(this%forest%leaf(li) /= oldLeaf(li)) changed = 1 enddo endif if(changed == 0) then deallocate(oldLeaf) return endif ! ---- 4. Transfer plan, new surface mesh, levels, extruded mesh and geometries ---- call BuildTransferPlan(this%forest,nOld,oldLeaf,plan) deallocate(oldLeaf) allocate(newSurf) call EmitMesh(this%forest,this%baseSurf,newSurf) allocate(zLevels(1:newSurf%nGeo+1,1:newSurf%nGeo+1,0:nz,1:newSurf%nElem)) call this%levels(newSurf,zLevels) allocate(newMesh) call newMesh%ExtrudedMesh(newSurf,zLevels,this%bottomBC,this%topBC) deallocate(zLevels) allocate(newGeom) call newGeom%Init(this%interp,newMesh%nElem) call newGeom%GenerateFromMesh(newMesh) newGeom2 => null() if(this%hasModel2D) then allocate(newGeom2) call newGeom2%Init(this%interp,newSurf%nElem) call newGeom2%GenerateFromMesh(newSurf) endif ! ---- 5. Migrate or stage, regrid, and transfer ---- cFirst = newSurf%decomp%offsetElem(newSurf%decomp%rankId+1)+1 cLast = newSurf%decomp%offsetElem(newSurf%decomp%rankId+2) nR = this%activeSurf%decomp%nRanks if(nR > 1) then if(.not. allocated(this%winFirst)) then allocate(this%winFirst(1:nR),this%winLast(1:nR),this%winFirst3(1:nR),this%winLast3(1:nR)) endif call PlanWindows(plan,nR,newSurf%decomp%offsetElem,this%winFirst,this%winLast) ! Columns are contiguous runs of nz elements, so a window of columns is a window of 3-D ! elements (an empty window wFirst > wLast stays empty). this%winFirst3 = (this%winFirst-1)*nz+1 this%winLast3 = this%winLast*nz wF = this%winFirst(newSurf%decomp%rankId+1) wL = this%winLast(newSurf%decomp%rankId+1) if(cLast < cFirst) then wF = 1 wL = 0 endif wF3 = (wF-1)*nz+1 wL3 = wL*nz call MigrateOldWindow_DGModel3D_t(model3d,this%winFirst3,this%winLast3,wF3,wL3, & this%nMigrateBytesRecv,this%nMigrateBytesSent, & this%nMigrateElemRemote) if(present(model2d)) then call model2d%MigrateOldWindow(this%winFirst,this%winLast,wF,wL,this%nMigrateBytesRecv, & this%nMigrateBytesSent,this%nMigrateElemRemote) endif else call StageSolutionForTransfer_DGModel3D_t(model3d) if(present(model2d)) call model2d%StageSolutionForTransfer() endif call model3d%Regrid(newMesh,newGeom) if(present(model2d)) call model2d%Regrid(newSurf,newGeom2) if(cLast >= cFirst) then call ApplyExtrudedTransferPlan_DGModel3D_t(model3d,plan,nz,this%interp,cFirst,cLast) if(present(model2d)) call model2d%ApplyTransferPlan(plan,this%interp,cFirst,cLast) endif call plan%Free() ! ---- 6. Retire the previous meshes/geometries and re-size the indicator ---- call ReleaseActive(this) this%activeSurf => newSurf this%activeMesh => newMesh this%activeGeom => newGeom if(this%hasModel2D) this%activeGeom2 => newGeom2 this%ownsActive = .true. call this%indicator%Free() call this%indicator%Init(this%interp,newMesh%nElem,this%refineThreshold,this%coarsenThreshold) call this%ApplyIndicatorSettings() adapted = .true. endsubroutine Adapt_ExtrudedAMRController function RecommendedTimeStep_ExtrudedAMRController(this,dtBase) result(dt) !! Level-based explicit time step bound: dtBase / 2**maxLevel, with dtBase stable on the !! base (level-0) mesh. Refinement halves the horizontal element size per level; the !! vertical spacing is unchanged, so this is conservative when the vertical CFL binds. implicit none class(ExtrudedAMRController),intent(in) :: this real(prec),intent(in) :: dtBase real(prec) :: dt dt = dtBase/real(2**this%maxLevel,prec) endfunction RecommendedTimeStep_ExtrudedAMRController endmodule SELF_AMRController_Extruded