Adapt_AMRController3D Subroutine

  subroutine Adapt_AMRController3D(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 (the 3-D analogue of the 2-D 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_DGModel3D_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. The multi-rank path still transfers on the host (see step 6),
    !! so both mirrors are fresh there.
    implicit none
    class(AMRController3D),intent(inout) :: this
    class(DGModel3D_t),intent(inout) :: model
    logical,intent(out) :: adapted
    ! Local
    integer :: li,s,pass,node,nbr,ns,nf,nOld,Np,changed,eFirst,eLast,iv
    integer,allocatable :: flag(:),spread(:)
    integer,allocatable :: oldLeaf(:)
    integer,allocatable :: leafIdx(:)
    type(TransferPlan3D) :: plan
    type(Mesh3D),pointer :: newMesh
    type(SEMHex),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,6
          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. 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. That moves the
    ! per-element tensor-product interpolation onto the device - which is where the measured
    ! saving comes from - and incidentally leaves no solution data crossing the host link.
    !
    ! On several ranks the old field must first be assembled globally, because each rank then
    ! fills exactly its new contiguous element range and elements that changed ranks are
    ! migrated by construction (the 2-D Stage-5 v1 migration). That allgather is a host
    ! operation, so the multi-rank path stays on the portable host transfer: a device transfer
    ! only pays off there once migration is point-to-point (Stage-5 v2).
    eFirst = -1 ! set below, once newMesh's decomposition is known
    if(model%mesh%decomp%nRanks > 1) then
      Np = this%interp%N+1
      allocate(uOld(1:Np,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*Np, &
                                   model%solution%interior(:,:,:,:,iv),uOld(:,:,:,:,iv))
      enddo

      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_AMRController3D