! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
!
! 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_2D
!! Adaptive-mesh-refinement controller for 2-D DG models: the driver that closes the serial
!! AMR loop around a live, time-stepping model. One call to Adapt performs one adaptation
!! epoch between time steps:
!!
!!   1. estimate  - the Legendre modal-decay indicator flags each element refine/keep/coarsen
!!                  from the current solution (driving variable ivar, e.g. pressure);
!!   2. cap       - refine flags on leaves already at maxLevel are demoted to keep, bounding
!!                  the finest resolution (and the stable time step) a priori;
!!   3. halo      - refine flags spread to face neighbours for nHalo passes, so a feature
!!                  moving at speed c stays inside the refined band provided the adaptation
!!                  cadence satisfies k*dt*c <= nHalo * h_fine;
!!   4. mutate    - AdaptFromFlags + Balance2to1 update the quad-forest; if the leaf set is
!!                  unchanged the epoch is a no-op and the model is untouched;
!!   5. transfer  - BuildTransferPlan maps old leaves to new (copy / exact prolongation /
!!                  conservative restriction), EmitMesh produces the solver-ready
!!                  nonconforming mesh, and a new SEMQuad geometry is generated;
!!   6. regrid    - model%Regrid rebinds the model's storage and boundary conditions to the
!!                  new mesh, preserving its time state and parameters; the transferred
!!                  solution is applied and uploaded to the device.
!!
!! The controller owns the forest, the indicator, and the meshes/geometries it emits. The mesh
!! is rebuilt each epoch and the previous one freed after the model is rebound. Geometry instead
!! lives in two long-lived buffers that alternate (AMR Stage 6c), so it is resized rather than
!! reallocated and the previous epoch's geometry remains available while the new one is built.
!! The base mesh and geometry the model was initialized with belong to the caller and are
!! never freed here, but the base mesh must outlive the controller: it supplies the
!! boundary-condition metadata and the communicator for every emitted mesh. After
!! controller%Free the model's mesh/geometry pointers are dangling, so free or stop using the
!! model first.
!!
!! MPI (AMR Stage 5): the forest is rank-replicated. At Init the global base-mesh tables are
!! allgathered so every rank builds an identical forest; each epoch the rank-local indicator
!! flags are allgathered (one small collective) so every rank applies identical mutations and
!! computes identical transfer plans and global connectivity. EmitMesh re-decomposes the new
!! leaf list into contiguous (space-filling-curve) ranges, so repartitioning and load balance
!! are implicit in every epoch. Solution migration gathers the old rank-local solutions into a
!! global field and applies the plan to the new rank-local range only - simple and correct at
!! single-node scale; a point-to-point exchange is a drop-in replacement behind the same
!! interface. All collectives run between time steps at the adaptation cadence; nothing is
!! added to the time-stepping loop.
!!
!! Because refinement halves the element scale per level, an explicit-stability time step
!! chosen for the base mesh must shrink with the finest active level;
!! RecommendedTimeStep(dtBase) = dtBase / 2**MaxLevel gives the level-based bound to pass to
!! ForwardStep after each epoch.

  use iso_fortran_env,only:int64
  use SELF_Constants
  use SELF_Lagrange
  use SELF_Mesh_2D
  use SELF_Geometry_2D
  use SELF_DGModel2D_t
  use SELF_QuadTreeMesh_2D
  use SELF_AdaptiveMesh_2D
  use SELF_RefinementIndicator_2D
  use SELF_TransferPlan_2D
  use SELF_DomainDecomposition
  use mpi

  implicit none

  !! Debug switches for the Stage 6c incremental geometry path, resolved once from the
  !! environment on first use. Present so that a suspected geometry problem can be split
  !! between the reuse copy and the compacted generation without rebuilding:
  !!
  !!   SELF_AMR_GEOM_NO_REUSE=1  regenerate every element (the reuse predicate never fires)
  !!   SELF_AMR_GEOM_VERIFY=1    additionally generate the full geometry the old way and compare
  !!                             element by element, reporting the first quantity that differs
  logical,save :: geomDebugResolved = .false.
  logical,save :: geomNoReuse = .false.
  logical,save :: geomVerify = .false.
  logical,save :: geomFull = .false. !! SELF_AMR_GEOM_FULL=1: bypass the incremental path

  type :: AMRController2D
    type(QuadTreeMesh2D) :: forest
    type(RefinementIndicator2D) :: indicator
    type(Mesh2D),pointer :: baseMesh => null() !! caller-owned; metadata source for EmitMesh
    type(Mesh2D),pointer :: activeMesh => null() !! mesh the model currently runs on
    type(SEMQuad),pointer :: activeGeom => null() !! geometry the model currently runs on
    type(Lagrange),pointer :: interp => null() !! the model's solution interpolant
    logical :: ownsActive = .false. !! whether activeMesh was emitted by us (vs the caller's)
    !! Two long-lived geometry buffers, alternated each epoch (AMR Stage 6c). Geometry is now
    !! resized in place rather than allocated and freed per epoch, which is what lets Stage 6b's
    !! amortization apply to it; alternating means the PREVIOUS epoch's geometry is still intact
    !! while the new one is filled, which the incremental reuse path needs. geomSlot records
    !! which buffer activeGeom currently is, or 0 while it is still the caller's geometry.
    type(SEMQuad),pointer :: geomA => null()
    type(SEMQuad),pointer :: geomB => null()
    integer :: geomSlot = 0
    !! Scratch geometry holding just the elements an epoch actually changed, generated compacted
    !! and then scattered into place. Persistent and resized, so it allocates only when an epoch
    !! changes more elements than any epoch before it.
    type(SEMQuad),pointer :: genGeom => null()
    !! Cumulative geometry-reuse accounting (AMR Stage 6c), so a run can report what fraction of
    !! elements avoided regeneration rather than leaving the payoff to be estimated.
    integer(int64) :: nGeomReused = 0
    integer(int64) :: nGeomGenerated = 0
    integer :: ivar = SELF_AMR_ALLVARS !! driving variable for the indicator
    integer :: maxLevel = 1 !! refinement-level cap
    integer :: nHalo = 1 !! refine-flag halo-expansion passes
    real(prec) :: refineThreshold = 0.0_prec
    real(prec) :: coarsenThreshold = 0.0_prec

  contains
    procedure,public :: Init => Init_AMRController2D
    procedure,public :: Free => Free_AMRController2D
    procedure,public :: Adapt => Adapt_AMRController2D
    procedure,public :: RecommendedTimeStep => RecommendedTimeStep_AMRController2D
    procedure,private :: NextGeomBuffer => NextGeomBuffer_AMRController2D
    procedure,private :: BuildGeometry => BuildGeometry_AMRController2D

  endtype AMRController2D

contains

  subroutine Init_AMRController2D(this,model,refineThreshold,coarsenThreshold,ivar, &
                                  maxLevel,nHalo)
    !! Attach the controller to an initialized model. The model's current mesh becomes the
    !! forest's base mesh (level 0); its geometry interpolant drives the indicator and the
    !! solution transfer. Thresholds are the sigma = log10 modal-energy-ratio cut-offs of the
    !! refinement indicator (refineThreshold > coarsenThreshold; see
    !! SELF_RefinementIndicator_2D). ivar is the driving solution variable (or
    !! SELF_AMR_ALLVARS). maxLevel >= 0 caps the refinement depth; nHalo >= 0 sets the
    !! refine-flag halo width in elements.
    implicit none
    class(AMRController2D),intent(out) :: this
    class(DGModel2D_t),intent(in) :: model
    real(prec),intent(in) :: refineThreshold
    real(prec),intent(in) :: coarsenThreshold
    integer,intent(in) :: ivar
    integer,intent(in) :: maxLevel
    integer,intent(in) :: nHalo

    if(.not. associated(model%mesh)) then
      print*,__FILE__,':',__LINE__, &
        ' : Error : AMRController2D%Init requires an initialized model.'
      stop 1
    endif
    if(maxLevel < 0 .or. nHalo < 0) then
      print*,__FILE__,':',__LINE__, &
        ' : Error : AMRController2D%Init requires maxLevel >= 0 and nHalo >= 0.'
      stop 1
    endif

    this%baseMesh => model%mesh
    this%activeMesh => model%mesh
    this%activeGeom => model%geometry
    this%interp => model%geometry%x%interp
    this%ownsActive = .false.
    this%ivar = ivar
    this%maxLevel = maxLevel
    this%nHalo = nHalo
    this%refineThreshold = refineThreshold
    this%coarsenThreshold = coarsenThreshold

    ! Rank-replicated forest: on one rank, straight from the mesh; on several, from the
    ! allgathered global base tables (every rank builds the identical forest).
    if(model%mesh%decomp%nRanks > 1) then
      call InitForestFromDecomposedMesh(this%forest,model%mesh)
    else
      call this%forest%Init(model%mesh)
    endif
    call this%indicator%Init(this%interp,model%mesh%nElem,refineThreshold,coarsenThreshold)

  endsubroutine Init_AMRController2D

  subroutine InitForestFromDecomposedMesh(forest,mesh)
    !! Build a rank-replicated forest over a decomposed base mesh: allgather the global
    !! node coordinates, side table, and material ids (by the decomposition's contiguous
    !! element ownership) and initialize the forest from the global tables. Collective over
    !! the mesh communicator; runs once, at controller initialization.
    implicit none
    type(QuadTreeMesh2D),intent(out) :: forest
    type(Mesh2D),intent(in) :: mesh
    ! Local
    integer :: nG,nGeo,r,s
    real(prec),allocatable :: coordsG(:,:,:,:)
    integer,allocatable :: siG(:,:,:),matG(:)
    integer,allocatable :: nbr(:,:),nbrSide(:,:),flip(:,:),bc(:,:)

    nG = mesh%decomp%nElem ! global element count
    nGeo = mesh%nGeo

    allocate(coordsG(1:2,1:nGeo+1,1:nGeo+1,1:nG))
    allocate(siG(1:5,1:4,1:nG))
    allocate(matG(1:nG))
    call AllgatherPerElemReals(mesh%decomp,2*(nGeo+1)*(nGeo+1), &
                               mesh%nodeCoords,coordsG)
    call AllgatherPerElemInts(mesh%decomp,20,mesh%sideInfo,siG)
    call AllgatherPerElemInts(mesh%decomp,1,mesh%elemMaterial,matG)

    ! Decode the global side table (sideInfo(3) already carries global element ids).
    allocate(nbr(1:4,1:nG),nbrSide(1:4,1:nG),flip(1:4,1:nG),bc(1:4,1:nG))
    do r = 1,nG
      do s = 1,4
        nbr(s,r) = siG(3,s,r)
        nbrSide(s,r) = siG(4,s,r)/10
        flip(s,r) = mod(siG(4,s,r),10)
        bc(s,r) = siG(5,s,r)
      enddo
    enddo

    call forest%InitGlobal(nG,nGeo,mesh%quadrature,coordsG,nbr,nbrSide,flip,bc,matG)

    deallocate(coordsG,siG,matG,nbr,nbrSide,flip,bc)

  endsubroutine InitForestFromDecomposedMesh

  subroutine AllgatherPerElemInts(decomp,perElem,localArr,globalArr)
    !! Allgather an integer array with perElem entries per element from the decomposition's
    !! contiguous rank-local element ranges into the global element ordering.
    implicit none
    type(DomainDecomposition),intent(in) :: decomp
    integer,intent(in) :: perElem
    integer,intent(in) :: localArr(*)
    integer,intent(out) :: globalArr(*)
    ! Local
    integer :: r,ierror
    integer,allocatable :: counts(:),displs(:)

    allocate(counts(1:decomp%nRanks),displs(1:decomp%nRanks))
    do r = 1,decomp%nRanks
      counts(r) = perElem*(decomp%offsetElem(r+1)-decomp%offsetElem(r))
      displs(r) = perElem*decomp%offsetElem(r)
    enddo
    call mpi_allgatherv(localArr,counts(decomp%rankId+1),MPI_INTEGER, &
                        globalArr,counts,displs,MPI_INTEGER, &
                        decomp%mpiComm,ierror)
    deallocate(counts,displs)

  endsubroutine AllgatherPerElemInts

  subroutine AllgatherPerElemReals(decomp,perElem,localArr,globalArr)
    !! Allgather a real(prec) array with perElem entries per element from the decomposition's
    !! contiguous rank-local element ranges into the global element ordering.
    implicit none
    type(DomainDecomposition),intent(in) :: decomp
    integer,intent(in) :: perElem
    real(prec),intent(in) :: localArr(*)
    real(prec),intent(out) :: globalArr(*)
    ! Local
    integer :: r,ierror
    integer,allocatable :: counts(:),displs(:)

    allocate(counts(1:decomp%nRanks),displs(1:decomp%nRanks))
    do r = 1,decomp%nRanks
      counts(r) = perElem*(decomp%offsetElem(r+1)-decomp%offsetElem(r))
      displs(r) = perElem*decomp%offsetElem(r)
    enddo
    call mpi_allgatherv(localArr,counts(decomp%rankId+1),decomp%mpiPrec, &
                        globalArr,counts,displs,decomp%mpiPrec, &
                        decomp%mpiComm,ierror)
    deallocate(counts,displs)

  endsubroutine AllgatherPerElemReals

  subroutine BuildGeometry_AMRController2D(this,newMesh,plan,newGeom,nReused)
    !! Fill newGeom for the emitted mesh, reusing the previous epoch's geometry for every element
    !! that did not change and generating only the rest (AMR Stage 6c). nReused reports how many
    !! elements were copied rather than computed.
    !!
    !! Why the reuse is exact. The transfer plan marks a new leaf SELF_TRANSFER_COPY only when the
    !! walk depth is zero, i.e. when the new element and old element sourceElem(li) are the SAME
    !! forest node (see BuildTransferPlan). A leaf's mesh node coordinates come from LeafCoords,
    !! which is a pure function of the root element's coordinates, the leaf's level and its
    !! quadrant path, evaluated in a fixed order; root coordinates are never mutated and node ids,
    !! levels and quadrants are stable across forest mutations. So an unchanged leaf's coordinates
    !! are bit-identical between epochs, and because per-element geometry generation touches only
    !! that element's own coordinates, its whole geometry block is too.
    !!
    !! Multi-rank: sourceElem indexes the GLOBAL old element list while each rank holds only its
    !! own slice, so reuse additionally requires the source to be locally owned. That is a range
    !! test against the OLD decomposition, which is still valid here because Regrid has not run
    !! yet. On one rank every COPY element qualifies and the test is always true; on several ranks
    !! whatever migrated is simply regenerated. No communication is added either way.
    implicit none
    class(AMRController2D),intent(inout) :: this
    type(Mesh2D),intent(in) :: newMesh
    type(TransferPlan2D),intent(in) :: plan
    type(SEMQuad),intent(inout) :: newGeom
    integer,intent(out) :: nReused
    ! Local
    integer :: li,gi,src,nLocal,nGen,oldFirst,oldLast,eFirst,nGeo,k
    integer,allocatable :: srcIdx(:),dstIdx(:),genIdx(:)
    real(prec),allocatable :: genCoords(:,:,:,:)

    nLocal = newMesh%nElem
    nGeo = newMesh%nGeo
    eFirst = newMesh%decomp%offsetElem(newMesh%decomp%rankId+1)+1

    ! Rank-local range of the OLD element list, i.e. what activeGeom actually holds.
    oldFirst = this%activeMesh%decomp%offsetElem(this%activeMesh%decomp%rankId+1)+1
    oldLast = this%activeMesh%decomp%offsetElem(this%activeMesh%decomp%rankId+2)

    allocate(srcIdx(1:nLocal),dstIdx(1:nLocal),genIdx(1:nLocal))

    call ResolveGeomDebug()

    ! Diagnostic bypass: reproduce the pre-6c behaviour exactly (full regeneration on the target
    ! buffer) while keeping the persistent alternating buffers, to tell a defect in the
    ! incremental assembly apart from one in the buffer reuse itself.
    if(geomFull) then
      call newGeom%GenerateFromMesh(newMesh)
      nReused = 0
      deallocate(srcIdx,dstIdx,genIdx)
      return
    endif

    nReused = 0
    nGen = 0
    do li = 1,nLocal
      gi = eFirst+li-1 ! this element's index in the plan's global new-leaf arrays
      src = plan%sourceElem(gi)
      if(.not. geomNoReuse .and. plan%sourceKind(gi) == SELF_TRANSFER_COPY .and. &
         src >= oldFirst .and. src <= oldLast) then
        nReused = nReused+1
        srcIdx(nReused) = src-oldFirst+1 ! rank-local index into activeGeom
        dstIdx(nReused) = li
      else
        nGen = nGen+1
        genIdx(nGen) = li
      endif
    enddo

    ! Generate the changed elements, compacted, so the generation loops run over nGen elements
    ! instead of all of them. Their geometry is then scattered into place.
    if(nGen > 0) then
      allocate(genCoords(1:2,1:nGeo+1,1:nGeo+1,1:nGen))
      do k = 1,nGen
        genCoords(1:2,:,:,k) = newMesh%nodeCoords(1:2,:,:,genIdx(k))
      enddo

      if(.not. associated(this%genGeom)) allocate(this%genGeom)
      if(this%genGeom%nElem == 0) then
        call this%genGeom%Init(this%interp,nGen)
      else
        call this%genGeom%Resize(this%interp,nGen)
      endif
      call this%genGeom%GenerateFromNodeCoords(genCoords,nGeo,newMesh%quadrature,nGen)
      call this%genGeom%x%UpdateDevice()
      call this%genGeom%x%BoundaryInterp()
      call this%genGeom%x%UpdateHost()
      call this%genGeom%CalculateMetricTerms()

      do k = 1,nGen
        srcIdx(nReused+k) = k
        dstIdx(nReused+k) = genIdx(k)
      enddo
      call newGeom%CopyElements(this%genGeom,srcIdx(nReused+1:),dstIdx(nReused+1:),nGen)
      deallocate(genCoords)
    endif

    ! Carry the unchanged elements across from the previous epoch's geometry.
    if(nReused > 0) then
      call newGeom%CopyElements(this%activeGeom,srcIdx,dstIdx,nReused)
    endif

    call newGeom%UploadGeometry()

    if(geomVerify) call VerifyGeometry(this,newMesh,newGeom)

    deallocate(srcIdx,dstIdx,genIdx)

  endsubroutine BuildGeometry_AMRController2D

  subroutine ResolveGeomDebug()
    !! Read the Stage 6c geometry debug switches once.
    implicit none
    character(8) :: envstr
    integer :: envstat

    if(geomDebugResolved) return
    geomDebugResolved = .true.

    call get_environment_variable("SELF_AMR_GEOM_NO_REUSE",envstr,status=envstat)
    if(envstat == 0) geomNoReuse = (trim(envstr) == "1")
    call get_environment_variable("SELF_AMR_GEOM_VERIFY",envstr,status=envstat)
    if(envstat == 0) geomVerify = (trim(envstr) == "1")
    call get_environment_variable("SELF_AMR_GEOM_FULL",envstr,status=envstat)
    if(envstat == 0) geomFull = (trim(envstr) == "1")

    if(geomNoReuse) print*,"SELF_AMR_GEOM_NO_REUSE: geometry reuse disabled"
    if(geomVerify) print*,"SELF_AMR_GEOM_VERIFY: cross-checking incremental geometry"
    if(geomFull) print*,"SELF_AMR_GEOM_FULL: incremental geometry bypassed"

  endsubroutine ResolveGeomDebug

  subroutine VerifyGeometry(this,newMesh,newGeom)
    !! Generate the geometry the original way and report, per quantity, the largest discrepancy
    !! against the incrementally built one. Diagnostic only; gated by SELF_AMR_GEOM_VERIFY.
    !!
    !! Whole-array maxima rather than per-element reporting: naming the quantity that diverges is
    !! what localizes the defect, and it keeps this to F2008 array intrinsics.
    implicit none
    class(AMRController2D),intent(inout) :: this
    type(Mesh2D),intent(in) :: newMesh
    type(SEMQuad),intent(in) :: newGeom
    ! Local
    type(SEMQuad) :: ref

    call ref%Init(this%interp,newMesh%nElem)
    call ref%GenerateFromMesh(newMesh)

    print*,"GEOM_VERIFY nElem            =",newMesh%nElem
    print*,"GEOM_VERIFY max|d x%interior|     =", &
      maxval(abs(newGeom%x%interior-ref%x%interior))
    print*,"GEOM_VERIFY max|d x%boundary|     =", &
      maxval(abs(newGeom%x%boundary-ref%x%boundary))
    print*,"GEOM_VERIFY max|d J%interior|     =", &
      maxval(abs(newGeom%J%interior-ref%J%interior))
    print*,"GEOM_VERIFY max|d J%boundary|     =", &
      maxval(abs(newGeom%J%boundary-ref%J%boundary))
    print*,"GEOM_VERIFY max|d dsdx%interior|  =", &
      maxval(abs(newGeom%dsdx%interior-ref%dsdx%interior))
    print*,"GEOM_VERIFY max|d dsdx%boundary|  =", &
      maxval(abs(newGeom%dsdx%boundary-ref%dsdx%boundary))
    print*,"GEOM_VERIFY max|d nHat%boundary|  =", &
      maxval(abs(newGeom%nHat%boundary-ref%nHat%boundary))
    print*,"GEOM_VERIFY max|d nScale%boundary|=", &
      maxval(abs(newGeom%nScale%boundary-ref%nScale%boundary))

    call ref%Free()

  endsubroutine VerifyGeometry

  subroutine NextGeomBuffer_AMRController2D(this,nElem,geom,slot)
    !! Select the geometry buffer to fill this epoch: whichever of the two is not currently
    !! active, so the previous epoch's geometry stays intact and readable (AMR Stage 6c).
    !!
    !! Each buffer is Init-ed the first time it is used and Resized on every subsequent epoch, so
    !! after the first two adaptations geometry performs no allocation at all - which is the point
    !! of the change. On the very first adaptation activeGeom is still the caller's geometry
    !! (geomSlot == 0); that object belongs to the caller and is never freed or reused here.
    implicit none
    class(AMRController2D),intent(inout) :: this
    integer,intent(in) :: nElem
    type(SEMQuad),pointer,intent(out) :: geom
    integer,intent(out) :: slot

    if(this%geomSlot == 1) then
      slot = 2
    else
      slot = 1
    endif

    if(slot == 1) then
      if(.not. associated(this%geomA)) allocate(this%geomA)
      geom => this%geomA
    else
      if(.not. associated(this%geomB)) allocate(this%geomB)
      geom => this%geomB
    endif

    if(geom%nElem == 0) then
      call geom%Init(this%interp,nElem)
    else
      call geom%Resize(this%interp,nElem)
    endif

  endsubroutine NextGeomBuffer_AMRController2D

  subroutine Free_AMRController2D(this)
    !! Release the forest, the indicator, and any controller-emitted mesh/geometry. The
    !! caller-owned base mesh/geometry are untouched. A model still pointing at a
    !! controller-emitted mesh must not be used after this call.
    implicit none
    class(AMRController2D),intent(inout) :: this

    call this%forest%Free()
    call this%indicator%Free()
    if(this%ownsActive) then
      call this%activeMesh%Free()
      deallocate(this%activeMesh)
    endif
    if(associated(this%geomA)) then
      call this%geomA%Free()
      deallocate(this%geomA)
    endif
    if(associated(this%geomB)) then
      call this%geomB%Free()
      deallocate(this%geomB)
    endif
    if(associated(this%genGeom)) then
      call this%genGeom%Free()
      deallocate(this%genGeom)
    endif
    this%baseMesh => null()
    this%activeMesh => null()
    this%activeGeom => null()
    this%geomA => null()
    this%geomB => null()
    this%genGeom => null()
    this%geomSlot = 0
    this%nGeomReused = 0
    this%nGeomGenerated = 0
    this%interp => null()
    this%ownsActive = .false.

  endsubroutine Free_AMRController2D

  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
    class(AMRController2D),intent(inout) :: 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,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.
    call this%indicator%Estimate(model%solution,this%ivar)
    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 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 (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: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

      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)

    adapted = .true.

  endsubroutine Adapt_AMRController2D

  function RecommendedTimeStep_AMRController2D(this,dtBase) result(dt)
    !! Level-based explicit-stability time step: refinement halves the element scale per
    !! level, so a time step dtBase that is stable on the base (level-0) mesh scales to
    !! dtBase / 2**MaxLevel on the current forest. Deterministic and exact for the quadtree
    !! (child elements are exact half-scale subdivisions); no geometry reduction is needed.
    implicit none
    class(AMRController2D),intent(in) :: this
    real(prec),intent(in) :: dtBase
    real(prec) :: dt

    dt = dtBase/real(2**this%forest%MaxLevel(),prec)

  endfunction RecommendedTimeStep_AMRController2D

endmodule SELF_AMRController_2D
