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.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(ExtrudedAMRController), | intent(inout), | target | :: | this | ||
| class(DGModel3D_t), | intent(inout) | :: | model3d | |||
| logical, | intent(out) | :: | adapted | |||
| class(DGModel2D_t), | intent(inout), | optional | :: | model2d |
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