Perform one adaptation epoch on the model (see the module documentation). On return, adapted reports whether the mesh changed; when it did, the model is already rebound to the new mesh with the solution transferred (conservatively), and the caller should re-evaluate its time step (RecommendedTimeStep) before the next ForwardStep. When the leaf set is unchanged the model is untouched.
Where the transferred solution lives on return: on a GPU build the transfer is performed on the device (Stage 6a) and the result is left there, so solution%interior (the host mirror) is STALE afterwards. This matches the rest of the time loop, where the device is authoritative and a caller that wants host data calls solution%UpdateHost() first - as Write_DGModel2D_t does before writing a snapshot. Before the device transfer existed the mirror happened to be fresh here because the transfer ran on the host; do not rely on that. On CPU builds host and device are the same storage and the question does not arise.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(AMRController2D), | intent(inout), | target | :: | this | ||
| class(DGModel2D_t), | intent(inout) | :: | model | |||
| logical, | intent(out) | :: | adapted |
subroutine Adapt_AMRController2D(this,model,adapted)
!! Perform one adaptation epoch on the model (see the module documentation). On return,
!! adapted reports whether the mesh changed; when it did, the model is already rebound to
!! the new mesh with the solution transferred (conservatively), and the caller should
!! re-evaluate its time step (RecommendedTimeStep) before the next ForwardStep. When the
!! leaf set is unchanged the model is untouched.
!!
!! Where the transferred solution lives on return: on a GPU build the transfer is performed
!! on the device (Stage 6a) and the result is left there, so solution%interior (the host
!! mirror) is STALE afterwards. This matches the rest of the time loop, where the device is
!! authoritative and a caller that wants host data calls solution%UpdateHost() first - as
!! Write_DGModel2D_t does before writing a snapshot. Before the device transfer existed the
!! mirror happened to be fresh here because the transfer ran on the host; do not rely on
!! that. On CPU builds host and device are the same storage and the question does not arise.
implicit none
! target: the migration window buffer is a component of this, and the windowed apply is fed
! through a pointer remapped onto it (see step 6), which requires the target attribute here.
class(AMRController2D),intent(inout),target :: this
class(DGModel2D_t),intent(inout) :: model
logical,intent(out) :: adapted
! Local
integer :: li,s,pass,node,nbr,ns,nf,nOld,Np,changed,eFirst,eLast,iv
integer :: nR,nWin,wFirst,wLast
real(prec),pointer :: uWin(:,:,:,:)
integer,allocatable :: flag(:),spread(:)
integer,allocatable :: oldLeaf(:)
integer,allocatable :: leafIdx(:)
type(TransferPlan2D) :: plan
type(Mesh2D),pointer :: newMesh
type(SEMQuad),pointer :: newGeom
integer :: newSlot
integer :: nReused
real(prec),allocatable :: uOld(:,:,:,:)
adapted = .false.
if(.not. associated(model%mesh,this%activeMesh)) then
print*,__FILE__,':',__LINE__, &
' : Error : the model is not running on this controller''s active mesh.'
stop 1
endif
! ---- 1. Indicator flags from the current solution ----
! The indicator is rank-local; the (replicated) forest needs the global per-leaf flags, so
! on nRanks > 1 they are allgathered by the active decomposition's element ranges. From
! here on every rank applies identical mutations to its identical forest copy.
! The indicator's amplitude gate normalizes each element's energy by the peak element energy,
! so on several ranks that peak must be global: otherwise the gate - and hence the flags and
! the adapted mesh - would depend on the decomposition. That is one extra small collective per
! epoch, alongside the flag allgather below, and none inside the time-stepping loop.
if(model%mesh%decomp%mpiEnabled) then
call this%indicator%Estimate(model%solution,this%ivar, &
comm=model%mesh%decomp%mpiComm)
else
call this%indicator%Estimate(model%solution,this%ivar)
endif
nOld = this%forest%nLeaves
allocate(flag(1:nOld))
if(model%mesh%decomp%nRanks > 1) then
if(model%mesh%decomp%nElem /= nOld) then
print*,__FILE__,':',__LINE__, &
' : Error : the active decomposition does not span the forest leaf list.'
stop 1
endif
call AllgatherPerElemInts(model%mesh%decomp,1,this%indicator%flag,flag)
else
flag(1:nOld) = this%indicator%flag(1:nOld)
endif
! ---- 2. Cap refinement at maxLevel ----
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
! ---- 3. Halo expansion: spread refine flags to face neighbours ----
! A leaf neighbour of a refine-flagged leaf is also flagged (up to the level cap) so the
! refined band extends nHalo elements beyond where the indicator fires; internal (finer)
! neighbours are already refined and need nothing.
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 ! physical boundary
if(this%forest%child(1,nbr) /= 0) cycle ! finer neighbour, already refined
if(this%forest%level(nbr) >= this%maxLevel) cycle ! at the cap
spread(leafIdx(nbr)) = SELF_AMR_REFINE
enddo
enddo
flag(1:nOld) = spread(1:nOld)
enddo
deallocate(spread,leafIdx)
! ---- 4. Mutate the forest; detect a no-op 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
! ---- 5. Transfer plan, emitted mesh, and geometry ----
call BuildTransferPlan(this%forest,nOld,oldLeaf,plan)
deallocate(oldLeaf)
allocate(newMesh)
call EmitMesh(this%forest,this%baseMesh,newMesh)
! Take the geometry buffer that is NOT currently active, so the previous epoch's geometry
! stays readable (Stage 6c). Each buffer is Init-ed once and resized thereafter.
call this%NextGeomBuffer(newMesh%nElem,newGeom,newSlot)
call this%BuildGeometry(newMesh,plan,newGeom,nReused)
this%nGeomReused = this%nGeomReused+int(nReused,int64)
this%nGeomGenerated = this%nGeomGenerated+int(newMesh%nElem-nReused,int64)
! ---- 6. Regrid the model and transfer (migrate) the solution ----
! The solution is staged before Regrid (which releases the storage it lives in) and
! transferred onto the new mesh afterwards. Both steps are type-bound and backend-specific:
! the portable implementation stages on the host and runs ApplyTransferPlanRange, while the
! GPU backend stages device-to-device and applies the plan in a kernel, so an adapting run
! on one GPU moves no solution data across the host link at all (Stage 6a).
!
! On several ranks the rank that will own a new element and the rank that owned its old
! source need not be the same, so old-field data has to move. Both partitions are contiguous
! ranges of the same leaf order and the plan is rank-replicated, so each rank computes the
! contiguous WINDOW of old elements its own new range references (PlanWindows) and receives
! exactly that window point-to-point (ExchangeOldWindow) - Stage-5 v2. Per-rank traffic and
! memory are then set by what actually moves rather than by the size of the global field.
! SELF_AMR_MIGRATE_GATHER=1 selects the v1 allgather instead; the two are bit-identical, and
! SELF_AMR_MIGRATE_VERIFY=1 asserts that in-process.
!
! The migration is device-resident on a GPU build (#172): the window is assembled in device
! memory and the plan applied to it by the kernel, so the multi-rank path moves no solution
! data across the host link either, and the per-element interpolation never runs on a CPU
! core. That rests on GPU-aware MPI, which the per-step halo exchange has always required.
! SELF_AMR_TRANSFER_HOST=1 forces the portable windowed path back on.
call ResolveGeomDebug()
eFirst = -1 ! set below, once newMesh's decomposition is known
if(model%mesh%decomp%nRanks > 1 .and. .not. migrateGather) then
Np = this%interp%N+1
nR = model%mesh%decomp%nRanks
if(.not. allocated(this%winFirst)) allocate(this%winFirst(1:nR),this%winLast(1:nR))
call PlanWindows(plan,nR,newMesh%decomp%offsetElem,this%winFirst,this%winLast)
eFirst = newMesh%decomp%offsetElem(newMesh%decomp%rankId+1)+1
eLast = newMesh%decomp%offsetElem(newMesh%decomp%rankId+2)
wFirst = this%winFirst(newMesh%decomp%rankId+1)
wLast = this%winLast(newMesh%decomp%rankId+1)
if(eLast < eFirst) then ! this rank owns no new elements: empty window
wFirst = 1
wLast = 0
endif
! Migrate the window, then apply the plan from it. Both steps are type-bound and
! backend-specific, and that is the whole point: the portable implementation migrates into
! host memory and runs the windowed apply on the host, while the GPU backend assembles the
! window in DEVICE memory - its own run device-to-device, the peers' runs received straight
! into it - and applies the plan as a kernel. So on a GPU build the per-element
! tensor-product interpolation runs on the device on any number of ranks, and no solution
! data crosses the host link. SELF_AMR_TRANSFER_HOST=1 forces the portable path back on,
! which is how the two are timed against each other in one binary.
!
! The migration runs BEFORE Regrid: EmitMesh has already decomposed the new mesh, so both
! partitions are known, the sends can read the still-live pre-regrid solution, and no
! point-to-point traffic is left in flight across Regrid - which works on the same
! communicator with its own tag conventions.
if(transferHost) then
call MigrateOldWindow_DGModel2D_t(model,this%winFirst,this%winLast, &
wFirst,wLast,this%nMigrateBytesRecv, &
this%nMigrateBytesSent,this%nMigrateElemRemote)
else
call model%MigrateOldWindow(this%winFirst,this%winLast,wFirst,wLast, &
this%nMigrateBytesRecv,this%nMigrateBytesSent, &
this%nMigrateElemRemote)
endif
! Either diagnostic needs the window on the host, and needs it BEFORE the apply, which
! consumes the migration marker. One download serves both.
if(migrateVerify .or. transferVerify) then
nWin = Np*Np*max(wLast-wFirst+1,0)*model%nvar
if(.not. allocated(this%xferWin)) allocate(this%xferWin(1:max(nWin,1)))
if(nWin > size(this%xferWin)) then
deallocate(this%xferWin)
allocate(this%xferWin(1:nWin))
endif
uWin(1:Np,1:Np,wFirst:wLast,1:model%nvar) => this%xferWin(1:nWin)
if(transferHost) then
call DownloadOldWindow_DGModel2D_t(model,wFirst,wLast,uWin)
else
call model%DownloadOldWindow(wFirst,wLast,uWin)
endif
endif
if(migrateVerify) then
! Cross-check the migrated window against an allgathered reference, BIT FOR BIT. The
! window is model state now and may live in device memory, hence the download above; the
! comparison itself is unchanged and stays exact, because migration is pure data
! movement. Reads the pre-regrid old field, so it must run before Regrid.
call model%solution%UpdateHost()
call VerifyMigration(model%mesh%decomp,Np,model%nvar,nOld,model%solution%nElem, &
model%solution%interior,wFirst,wLast,uWin)
endif
call model%Regrid(newMesh,newGeom)
! A rank that owns no new elements has nothing to fill; it still took part in the
! migration above, because peers may need old elements it owns.
if(eLast >= eFirst) then
if(transferHost) then
call ApplyTransferPlan_DGModel2D_t(model,plan,this%interp,eFirst,eLast)
else
call model%ApplyTransferPlan(plan,this%interp,eFirst,eLast)
endif
endif
if(transferVerify .and. eLast >= eFirst) then
call VerifyWindowedApply(model,plan,this%interp,eFirst,eLast,wFirst,wLast,uWin)
endif
if(migrateVerify .or. transferVerify) uWin => null()
elseif(model%mesh%decomp%nRanks > 1) then
Np = this%interp%N+1
nR = model%mesh%decomp%nRanks
allocate(uOld(1:Np,1:Np,1:nOld,1:model%nvar))
call model%solution%UpdateHost()
do iv = 1,model%nvar
call AllgatherPerElemReals(model%mesh%decomp,Np*Np, &
model%solution%interior(:,:,:,iv),uOld(:,:,:,iv))
enddo
! v1 volume, counted the same way as v2 so the two are directly comparable: every rank
! receives every element it does not own, and sends its own elements to every other rank.
this%nMigrateElemRemote = this%nMigrateElemRemote+ &
int(nOld-model%solution%nElem,int64)
this%nMigrateBytesRecv = this%nMigrateBytesRecv+int(nOld-model%solution%nElem,int64)* &
Np*Np*model%nvar*(storage_size(1.0_prec)/8)
this%nMigrateBytesSent = this%nMigrateBytesSent+int(model%solution%nElem,int64)*(nR-1)* &
Np*Np*model%nvar*(storage_size(1.0_prec)/8)
call model%Regrid(newMesh,newGeom)
eFirst = newMesh%decomp%offsetElem(newMesh%decomp%rankId+1)+1
eLast = newMesh%decomp%offsetElem(newMesh%decomp%rankId+2)
call model%ApplyTransferPlan(plan,this%interp,eFirst,eLast,uOld)
deallocate(uOld)
else
call model%StageSolutionForTransfer()
call model%Regrid(newMesh,newGeom)
eFirst = newMesh%decomp%offsetElem(newMesh%decomp%rankId+1)+1
eLast = newMesh%decomp%offsetElem(newMesh%decomp%rankId+2)
call model%ApplyTransferPlan(plan,this%interp,eFirst,eLast)
endif
call plan%Free()
! ---- 7. Retire the previous mesh/geometry and re-size the indicator ----
! The geometry buffers are retained and reused; only the mesh is still rebuilt per epoch.
if(this%ownsActive) then
call this%activeMesh%Free()
deallocate(this%activeMesh)
endif
this%activeMesh => newMesh
this%activeGeom => newGeom
this%geomSlot = newSlot
this%ownsActive = .true.
call this%indicator%Free()
call this%indicator%Init(this%interp,newMesh%nElem, &
this%refineThreshold,this%coarsenThreshold)
! Only controller-owned indicator settings survive an epoch. A driver that called
! indicator%SetEnergyScale (or SetEnergyWeights) directly must re-apply it after any epoch
! that reported adapted = .true.
call this%ApplyIndicatorSettings()
adapted = .true.
endsubroutine Adapt_AMRController2D