SELF_TransferPlan_2D.f90 Source File

! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
!
! 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_TransferPlan_2D
!! Old-leaf -> new-leaf solution transfer plan for one adaptation epoch of the 2-D quad-forest:
!! the driver layer of AMR Stage 3 that connects the element-local transfer operators
!! (SELF_SolutionTransfer_2D) to the forest mutation (SELF_QuadTreeMesh_2D).
!!
!! One adaptation epoch mutates the forest between two leaf configurations:
!!
!!     nOld = forest%nLeaves ; oldLeaf = forest%leaf(1:nOld)   ! snapshot BEFORE mutating
!!     call forest%AdaptFromFlags(flag)                        ! at most one call, then
!!     call forest%Balance2to1()                               ! any number of further
!!                                                             ! refinements (RefineNode too)
!!     call BuildTransferPlan(forest,nOld,oldLeaf,plan)        ! AFTER the last mutation
!!
!! Solution element index = leaf-list position (the element ordering EmitMesh produces), so the
!! plan records, for every new leaf, where its data comes from in the old element ordering:
!!
!!   SELF_TRANSFER_COPY     : the leaf survived unchanged; copy old element sourceElem.
!!   SELF_TRANSFER_PROLONG  : the leaf descends from old leaf sourceElem; interpolate the old
!!                            element's degree-N polynomial down the quadtree path (one step per
!!                            level, so depth > 1 handles a fresh child that 2:1 balancing
!!                            refined again in the same epoch). Exact, no loss.
!!   SELF_TRANSFER_RESTRICT : the leaf is (an ancestor of) a coarsened family; L2-project the
!!                            four old children family(1:4) onto their parent, then prolong down
!!                            depth >= 0 further steps (depth > 0 occurs when a just-coarsened
!!                            parent is immediately re-refined by 2:1 balancing). Conservative.
!!
!! The reconstruction is possible after the fact because forest node ids are stable: refinement
!! only appends nodes and coarsening only detaches children, whose level/parent/quadrant entries
!! persist. Each new leaf therefore ascends its parent chain until it meets either an old leaf
!! or a node holding a complete four-child old-leaf family; anything else means the snapshot
!! does not describe the epoch that produced the forest, and the builder stops with an error.
!!
!! ApplyTransferPlan executes the plan on nodal data u(1:N+1,1:N+1,1:nElem,1:nVar) (the layout
!! of MappedScalar2D %interior; units are those of the transferred fields) and inherits the
!! Stage-3 operator identities: prolongation is exact polynomial interpolation and restriction
!! is the conservative L2 projection, so the Jacobian-weighted cell integral of every variable
!! is conserved and a refine-then-coarsen round trip is the identity to roundoff. Transfer runs
!! once per adaptation epoch, between time steps; it is not a per-step hot path.

  use SELF_Constants
  use SELF_Lagrange
  use SELF_QuadTreeMesh_2D
  use SELF_SolutionTransfer_2D

  implicit none

  integer,parameter :: SELF_TRANSFER_COPY = 0
  integer,parameter :: SELF_TRANSFER_PROLONG = 1
  integer,parameter :: SELF_TRANSFER_RESTRICT = 2

  type :: TransferPlan2D
    integer :: nOld = 0 !! old (pre-epoch) element count
    integer :: nNew = 0 !! new element count (= forest%nLeaves at build time)
    integer :: maxDepth = 0 !! largest prolongation depth in the plan
    integer,allocatable :: sourceKind(:) !! (nNew) SELF_TRANSFER_COPY / PROLONG / RESTRICT
    integer,allocatable :: sourceElem(:) !! (nNew) old element index (COPY/PROLONG; 0 otherwise)
    integer,allocatable :: family(:,:) !! (4,nNew) old element index per child quadrant (RESTRICT)
    integer,allocatable :: depth(:) !! (nNew) number of prolongation steps below the source
    integer,allocatable :: path(:,:) !! (>=maxDepth,nNew) quadrant taken at each step, top-down

  contains
    procedure,public :: Free => Free_TransferPlan2D

  endtype TransferPlan2D

contains

  subroutine Free_TransferPlan2D(this)
    implicit none
    class(TransferPlan2D),intent(inout) :: this

    if(allocated(this%sourceKind)) deallocate(this%sourceKind)
    if(allocated(this%sourceElem)) deallocate(this%sourceElem)
    if(allocated(this%family)) deallocate(this%family)
    if(allocated(this%depth)) deallocate(this%depth)
    if(allocated(this%path)) deallocate(this%path)
    this%nOld = 0
    this%nNew = 0
    this%maxDepth = 0

  endsubroutine Free_TransferPlan2D

  subroutine BuildTransferPlan(forest,nOld,oldLeaf,plan)
    !! Build the old->new transfer plan for the adaptation epoch that took the forest from the
    !! leaf configuration (nOld, oldLeaf) - a snapshot of (forest%nLeaves, forest%leaf) taken
    !! before mutating - to its current leaf configuration. See the module documentation for the
    !! allowed mutations within one epoch. The plan is valid until the forest is mutated again.
    implicit none
    type(QuadTreeMesh2D),intent(in) :: forest
    integer,intent(in) :: nOld
    integer,intent(in) :: oldLeaf(1:nOld)
    type(TransferPlan2D),intent(out) :: plan
    ! Local
    integer :: i,li,walk,d,mx,steps
    logical :: found
    integer,allocatable :: oldElemOfNode(:)
    integer,allocatable :: famElem(:,:)
    integer,allocatable :: rev(:)

    do i = 1,nOld
      if(oldLeaf(i) < 1 .or. oldLeaf(i) > forest%nNodes) then
        print*,__FILE__,':',__LINE__, &
          ' : Error : oldLeaf snapshot contains a node id outside the forest node pool.'
        stop 1
      endif
    enddo

    ! Invert the snapshot: old element index by node id, and each old family's element indices
    ! by (child quadrant, parent node id). Coarsening zeroes the parent's child pointers but
    ! never the children's parent/quadrant entries, so famElem recovers detached families.
    allocate(oldElemOfNode(1:forest%nNodes))
    oldElemOfNode = 0
    allocate(famElem(1:4,1:forest%nNodes))
    famElem = 0
    do i = 1,nOld
      oldElemOfNode(oldLeaf(i)) = i
      if(forest%parent(oldLeaf(i)) > 0) then
        famElem(forest%quadrant(oldLeaf(i)),forest%parent(oldLeaf(i))) = i
      endif
    enddo

    plan%nOld = nOld
    plan%nNew = forest%nLeaves
    mx = max(forest%MaxLevel(),1) ! prolongation depth is bounded by the deepest leaf level
    allocate(plan%sourceKind(1:plan%nNew))
    allocate(plan%sourceElem(1:plan%nNew))
    allocate(plan%family(1:4,1:plan%nNew))
    allocate(plan%depth(1:plan%nNew))
    allocate(plan%path(1:mx,1:plan%nNew))
    plan%sourceKind = SELF_TRANSFER_COPY
    plan%sourceElem = 0
    plan%family = 0
    plan%depth = 0
    plan%path = 0
    plan%maxDepth = 0

    allocate(rev(1:mx))

    do li = 1,plan%nNew

      ! Ascend from the new leaf until reaching its data source; collect the quadrant taken at
      ! each level (bottom-up in rev, reversed into plan%path top-down). The chain has at most
      ! level+1 nodes, ending at the root.
      walk = forest%leaf(li)
      d = 0
      found = .false.
      do steps = 0,forest%level(forest%leaf(li))
        if(oldElemOfNode(walk) > 0) then
          if(d == 0) then
            plan%sourceKind(li) = SELF_TRANSFER_COPY
          else
            plan%sourceKind(li) = SELF_TRANSFER_PROLONG
          endif
          plan%sourceElem(li) = oldElemOfNode(walk)
          found = .true.
          exit
        elseif(famElem(1,walk) > 0 .and. famElem(2,walk) > 0 .and. &
               famElem(3,walk) > 0 .and. famElem(4,walk) > 0) then
          plan%sourceKind(li) = SELF_TRANSFER_RESTRICT
          plan%family(1:4,li) = famElem(1:4,walk)
          found = .true.
          exit
        endif
        if(forest%parent(walk) == 0) exit
        d = d+1
        rev(d) = forest%quadrant(walk)
        walk = forest%parent(walk)
      enddo

      if(.not. found) then
        print*,__FILE__,':',__LINE__, &
          ' : Error : new leaf has no old-leaf ancestor or coarsened old family; the snapshot'// &
          ' does not describe one adaptation epoch of this forest.'
        stop 1
      endif

      plan%depth(li) = d
      do i = 1,d
        plan%path(i,li) = rev(d-i+1)
      enddo
      plan%maxDepth = max(plan%maxDepth,d)

    enddo

    deallocate(oldElemOfNode,famElem,rev)

  endsubroutine BuildTransferPlan

  subroutine ApplyTransferPlan(plan,interp,nVar,uOld,uNew)
    !! Execute a transfer plan on nodal element data: uNew(:,:,li,:) receives old element data
    !! copied, prolonged (exact interpolation), or restricted (conservative L2 projection)
    !! according to plan entry li. interp must be the solution interpolant the data lives on
    !! (its mortar operators drive the transfer). Runs once per adaptation epoch - not a
    !! per-time-step hot path - so clarity is preferred over fused loops here.
    implicit none
    type(TransferPlan2D),intent(in) :: plan
    type(Lagrange),intent(in) :: interp
    integer,intent(in) :: nVar
    real(prec),intent(in) :: uOld(1:interp%N+1,1:interp%N+1,1:plan%nOld,1:nVar)
    real(prec),intent(out) :: uNew(1:interp%N+1,1:interp%N+1,1:plan%nNew,1:nVar)

    call ApplyTransferPlanRange(plan,interp,nVar,uOld,1,plan%nNew,uNew)

  endsubroutine ApplyTransferPlan

  subroutine ApplyTransferPlanRange(plan,interp,nVar,uOld,eFirst,eLast,uNew)
    !! Execute the contiguous sub-range eFirst..eLast of a transfer plan: uNew(:,:,k,:)
    !! receives the data of new element eFirst+k-1. uOld is the full (global) old field; the
    !! output is only the requested slice. This is the decomposed-mesh (AMR Stage 5) entry
    !! point: each rank passes its own contiguous range of the new element ordering and a
    !! gathered global old solution, and fills exactly its rank-local storage.
    implicit none
    type(TransferPlan2D),intent(in) :: plan
    type(Lagrange),intent(in) :: interp
    integer,intent(in) :: nVar
    real(prec),intent(in) :: uOld(1:interp%N+1,1:interp%N+1,1:plan%nOld,1:nVar)
    integer,intent(in) :: eFirst
    integer,intent(in) :: eLast
    real(prec),intent(out) :: uNew(1:interp%N+1,1:interp%N+1,1:eLast-eFirst+1,1:nVar)
    ! Local
    integer :: li,lo,c,step,Np
    real(prec),allocatable :: buf(:,:,:),kids(:,:,:,:),fam(:,:,:,:)

    if(.not. allocated(plan%sourceKind)) then
      print*,__FILE__,':',__LINE__, &
        ' : Error : ApplyTransferPlan called with an unbuilt TransferPlan2D.'
      stop 1
    endif
    if(eFirst < 1 .or. eLast > plan%nNew .or. eLast < eFirst) then
      print*,__FILE__,':',__LINE__, &
        ' : Error : ApplyTransferPlanRange called with a range outside 1..nNew.'
      stop 1
    endif

    Np = interp%N+1
    allocate(buf(1:Np,1:Np,1:nVar))
    allocate(kids(1:Np,1:Np,1:nVar,1:4))
    allocate(fam(1:Np,1:Np,1:nVar,1:4))

    do li = eFirst,eLast
      lo = li-eFirst+1

      if(plan%sourceKind(li) == SELF_TRANSFER_COPY) then
        uNew(1:Np,1:Np,lo,1:nVar) = uOld(1:Np,1:Np,plan%sourceElem(li),1:nVar)
        cycle
      endif

      if(plan%sourceKind(li) == SELF_TRANSFER_RESTRICT) then
        do c = 1,4
          fam(1:Np,1:Np,1:nVar,c) = uOld(1:Np,1:Np,plan%family(c,li),1:nVar)
        enddo
        call RestrictFromChildren(interp,nVar,fam,buf)
      else ! SELF_TRANSFER_PROLONG
        buf(1:Np,1:Np,1:nVar) = uOld(1:Np,1:Np,plan%sourceElem(li),1:nVar)
      endif

      do step = 1,plan%depth(li)
        call ProlongToChildren(interp,nVar,buf,kids)
        buf(1:Np,1:Np,1:nVar) = kids(1:Np,1:Np,1:nVar,plan%path(step,li))
      enddo

      uNew(1:Np,1:Np,lo,1:nVar) = buf(1:Np,1:Np,1:nVar)

    enddo

    deallocate(buf,kids,fam)

  endsubroutine ApplyTransferPlanRange

endmodule SELF_TransferPlan_2D