! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
!
! 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_SolutionTransfer_3D
!! Solution transfer between an element and its eight h-refinement children (AMR Stage 3).
!!
!! When a 3-D hexahedral element is refined into eight children (SELF child ordering = CGNS
!! corner order, child c covering the octant with half-interval indices
!! (octAxc(c),octAyc(c),octAzc(c)) of the parent reference cube; see
!! SELF_RefinementPrimitives_3D), the prognostic solution must move with the mesh:
!!
!!   - **Prolongation** (parent -> 8 children): sample the parent's degree-N nodal polynomial at
!!     the children's nodes. Each child occupies a reference half-interval in each direction, so
!!     this is exactly a triple tensor product of the 2:1 mortar restriction operator
!!     Lagrange%mortarR (built and identity-checked for the nonconforming interface machinery).
!!     It is an exact interpolation: the child fields reproduce the parent polynomial with no
!!     loss.
!!
!!   - **Restriction** (8 children -> parent): the L2 projection of the (generally discontinuous
!!     across children) fine solution back onto the parent's degree-N polynomial space, a triple
!!     tensor product of the mortar projection operator Lagrange%mortarP (the M-inner-product
!!     adjoint of mortarR). Each 1-D mortarP carries the 1/2 per-direction sub-interval
!!     Jacobian, so the triple product carries the 1/8 sub-cell Jacobian appropriate for
!!     projecting *solution* fields: a conservative L2 projection.
!!
!! These inherit the discrete identities established for the 1-D mortar operators
!! (test/mortarprojection_identity.f90):
!!
!!   * Consistency / reversibility:  RestrictFromChildren(ProlongToChildren(u)) = u exactly, from
!!     sum_k P_k R_k = I applied in each tensor direction. Refining then immediately coarsening
!!     does not perturb the solution.
!!   * Conservation:  sum_ijk w_i w_j w_k u_parent(i,j,k) =
!!     (1/8) sum_c sum_ijk w_i w_j w_k u_child_c(i,j,k), i.e. the reference-cell integral of the
!!     restricted parent equals the sum of the children's reference-cell integrals (each child
!!     is an eighth of the parent). Weighted by the geometry Jacobian this is conservation of
!!     the cell-integrated prognostic quantity.
!!
!! The routines are element-local and portable (host `do concurrent` over children/variables);
!! the AMR driver maps forest parent/child relationships onto the element index ranges it passes
!! in. Transfer runs between time steps, so it is not a per-step hot path.
!!
!! These routines remain the reference implementation and are what the CPU backend and the
!! multi-rank migration path execute, and they are what the transfer tests pin. On a single-rank
!! GPU build the same mathematics run in TransferSolution_3D_gpu
!! (src/gpu/SELF_SolutionTransfer.cpp), so the contractions below are performed on the device
!! rather than on one CPU core and the transferred field never crosses the host link; see
!! ApplyTransferPlan_DGModel3D.

  use SELF_Constants
  use SELF_Lagrange

  implicit none

  ! Child octant -> (x-half, y-half, z-half) with 0 = lower ([-1,0]), 1 = upper ([0,1]).
  ! Half index h maps to mortar sub-interval k = h+1 (sub-interval 1 = [-1,0], 2 = [0,1]).
  ! Identical to octAxc/octAyc/octAzc in SELF_RefinementPrimitives_3D.
  integer,parameter :: transferAxc(1:8) = [0,1,1,0,0,1,1,0]
  integer,parameter :: transferAyc(1:8) = [0,0,1,1,0,0,1,1]
  integer,parameter :: transferAzc(1:8) = [0,0,0,0,1,1,1,1]

contains

  subroutine ProlongToChildren(interp,nVar,uParent,uChildren)
    !! Prolong (interpolate) a parent element's nodal solution onto its eight children.
    !! uChildren(:,:,:,:,c) is the field on child c (SELF child ordering = CGNS corner order).
    !! Exact for the parent's degree-N polynomial representation.
    implicit none
    type(Lagrange),intent(in) :: interp
    integer,intent(in) :: nVar
    real(prec),intent(in) :: uParent(1:interp%N+1,1:interp%N+1,1:interp%N+1,1:nVar)
    real(prec),intent(out) :: uChildren(1:interp%N+1,1:interp%N+1,1:interp%N+1,1:nVar,1:8)
    ! Local
    integer :: c,v,Np

    Np = interp%N+1

    do concurrent(c=1:8,v=1:nVar)
      block
        integer :: i,j,k,ii,jj,kk,kx,ky,kz
        real(prec) :: acc
        real(prec) :: tmp1(1:Np,1:Np,1:Np) ! tmp1(childX, parentY, parentZ) after the x pass
        real(prec) :: tmp2(1:Np,1:Np,1:Np) ! tmp2(childX, childY, parentZ) after the y pass

        kx = transferAxc(c)+1
        ky = transferAyc(c)+1
        kz = transferAzc(c)+1

        ! x-direction: tmp1(i,jj,kk) = sum_ii mortarR(ii,i,kx) * uParent(ii,jj,kk)
        do k = 1,Np
          do j = 1,Np
            do i = 1,Np
              acc = 0.0_prec
              do ii = 1,Np
                acc = acc+interp%mortarR(ii,i,kx)*uParent(ii,j,k,v)
              enddo
              tmp1(i,j,k) = acc
            enddo
          enddo
        enddo
        ! y-direction: tmp2(i,j,kk) = sum_jj mortarR(jj,j,ky) * tmp1(i,jj,kk)
        do k = 1,Np
          do j = 1,Np
            do i = 1,Np
              acc = 0.0_prec
              do jj = 1,Np
                acc = acc+interp%mortarR(jj,j,ky)*tmp1(i,jj,k)
              enddo
              tmp2(i,j,k) = acc
            enddo
          enddo
        enddo
        ! z-direction: uChild(i,j,k) = sum_kk mortarR(kk,k,kz) * tmp2(i,j,kk)
        do k = 1,Np
          do j = 1,Np
            do i = 1,Np
              acc = 0.0_prec
              do kk = 1,Np
                acc = acc+interp%mortarR(kk,k,kz)*tmp2(i,j,kk)
              enddo
              uChildren(i,j,k,v,c) = acc
            enddo
          enddo
        enddo
      endblock
    enddo

  endsubroutine ProlongToChildren

  subroutine RestrictFromChildren(interp,nVar,uChildren,uParent)
    !! Restrict (L2-project) the solution on eight children back onto their parent element.
    !! Conservative, and the exact left inverse of ProlongToChildren.
    implicit none
    type(Lagrange),intent(in) :: interp
    integer,intent(in) :: nVar
    real(prec),intent(in) :: uChildren(1:interp%N+1,1:interp%N+1,1:interp%N+1,1:nVar,1:8)
    real(prec),intent(out) :: uParent(1:interp%N+1,1:interp%N+1,1:interp%N+1,1:nVar)
    ! Local
    integer :: v,Np

    Np = interp%N+1

    do concurrent(v=1:nVar)
      block
        integer :: c,i,j,k,ii,jj,kk,kx,ky,kz
        real(prec) :: acc
        real(prec) :: tmp1(1:Np,1:Np,1:Np) ! tmp1(parentX, childY, childZ) after the x projection
        real(prec) :: tmp2(1:Np,1:Np,1:Np) ! tmp2(parentX, parentY, childZ) after the y projection
        real(prec) :: up(1:Np,1:Np,1:Np)

        up(1:Np,1:Np,1:Np) = 0.0_prec
        do c = 1,8
          kx = transferAxc(c)+1
          ky = transferAyc(c)+1
          kz = transferAzc(c)+1
          ! x-direction: tmp1(i,jj,kk) = sum_ii mortarP(ii,i,kx) * uChild(ii,jj,kk,c)
          do kk = 1,Np
            do jj = 1,Np
              do i = 1,Np
                acc = 0.0_prec
                do ii = 1,Np
                  acc = acc+interp%mortarP(ii,i,kx)*uChildren(ii,jj,kk,v,c)
                enddo
                tmp1(i,jj,kk) = acc
              enddo
            enddo
          enddo
          ! y-direction: tmp2(i,j,kk) = sum_jj mortarP(jj,j,ky) * tmp1(i,jj,kk)
          do kk = 1,Np
            do j = 1,Np
              do i = 1,Np
                acc = 0.0_prec
                do jj = 1,Np
                  acc = acc+interp%mortarP(jj,j,ky)*tmp1(i,jj,kk)
                enddo
                tmp2(i,j,kk) = acc
              enddo
            enddo
          enddo
          ! z-direction and accumulate over the eight children
          do k = 1,Np
            do j = 1,Np
              do i = 1,Np
                acc = 0.0_prec
                do kk = 1,Np
                  acc = acc+interp%mortarP(kk,k,kz)*tmp2(i,j,kk)
                enddo
                up(i,j,k) = up(i,j,k)+acc
              enddo
            enddo
          enddo
        enddo
        uParent(1:Np,1:Np,1:Np,v) = up(1:Np,1:Np,1:Np)
      endblock
    enddo

  endsubroutine RestrictFromChildren

endmodule SELF_SolutionTransfer_3D
