! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
!
! 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_ECDGModel3D_t
  !! Entropy-conserving DGSEM base class for 3-D conservation laws.
  !! See SELF_ECDGModel2D_t for design documentation.

  use SELF_DGModel3D_t
  use SELF_MappedTwoPointVector_3D

  implicit none

  type,extends(DGModel3D_t) :: ECDGModel3D_t

    type(MappedTwoPointVector3D) :: twoPointFlux

  contains

    procedure :: Init => Init_ECDGModel3D_t
    procedure :: Free => Free_ECDGModel3D_t

    procedure :: twopointflux3d => twopointflux3d_ECDGModel3D_t
    procedure :: TwoPointFluxMethod => TwoPointFluxMethod_ECDGModel3D_t
    procedure :: CalculateTendency => CalculateTendency_ECDGModel3D_t

  endtype ECDGModel3D_t

contains

  subroutine Init_ECDGModel3D_t(this,mesh,geometry)
    implicit none
    class(ECDGModel3D_t),intent(out) :: this
    type(Mesh3D),intent(in),target :: mesh
    type(SEMHex),intent(in),target :: geometry

    call Init_DGModel3D_t(this,mesh,geometry)

    call this%twoPointFlux%Init(geometry%x%interp,this%nvar,mesh%nElem)
    call this%twoPointFlux%AssociateGeometry(geometry)

  endsubroutine Init_ECDGModel3D_t

  subroutine Free_ECDGModel3D_t(this)
    implicit none
    class(ECDGModel3D_t),intent(inout) :: this

    call this%twoPointFlux%DissociateGeometry()
    call this%twoPointFlux%Free()
    call Free_DGModel3D_t(this)

  endsubroutine Free_ECDGModel3D_t

  pure function twopointflux3d_ECDGModel3D_t(this,sL,sR) result(flux)
    !! Entropy-conserving two-point flux function (3-D).
    !! flux(ivar, d) is the d-th physical component (d=1:x, d=2:y, d=3:z).
    !! Override in concrete models. Stub returns zero.
    implicit none
    class(ECDGModel3D_t),intent(in) :: this
    real(prec),intent(in) :: sL(1:this%nvar)
    real(prec),intent(in) :: sR(1:this%nvar)
    real(prec) :: flux(1:this%nvar,1:3)

    flux = 0.0_prec
    if(.false.) flux(1,1) = flux(1,1)+sL(1)+sR(1) ! suppress unused-argument warning

  endfunction twopointflux3d_ECDGModel3D_t

  subroutine TwoPointFluxMethod_ECDGModel3D_t(this)
    !! Computes pre-projected SCALAR contravariant two-point fluxes for all
    !! node pairs, following Trixi.jl for curved meshes. Each direction r
    !! uses the correct partner AND the correct averaged metric Ja^r.
    implicit none
    class(ECDGModel3D_t),intent(inout) :: this
    ! Local
    integer :: nn,i,j,k,d,iEl,iVar
    real(prec) :: sL(1:this%nvar),sR(1:this%nvar)
    real(prec) :: Fphys(1:this%nvar,1:3)
    real(prec) :: Fc

    do concurrent(nn=1:this%solution%N+1,i=1:this%solution%N+1, &
                  j=1:this%solution%N+1,k=1:this%solution%N+1, &
                  iEl=1:this%mesh%nElem)

      sL = this%solution%interior(i,j,k,iEl,1:this%nvar)

      ! xi^1: pair (i,j,k)-(nn,j,k), project onto avg(Ja^1)
      sR = this%solution%interior(nn,j,k,iEl,1:this%nvar)
      Fphys = this%twopointflux3d(sL,sR)
      do iVar = 1,this%nvar
        Fc = 0.0_prec
        do d = 1,3
          Fc = Fc+0.5_prec*( &
               this%geometry%dsdx%interior(i,j,k,iEl,1,d,1)+ &
               this%geometry%dsdx%interior(nn,j,k,iEl,1,d,1))* &
               Fphys(iVar,d)
        enddo
        this%twoPointFlux%interior(nn,i,j,k,iEl,iVar,1) = Fc
      enddo

      ! xi^2: pair (i,j,k)-(i,nn,k), project onto avg(Ja^2)
      sR = this%solution%interior(i,nn,k,iEl,1:this%nvar)
      Fphys = this%twopointflux3d(sL,sR)
      do iVar = 1,this%nvar
        Fc = 0.0_prec
        do d = 1,3
          Fc = Fc+0.5_prec*( &
               this%geometry%dsdx%interior(i,j,k,iEl,1,d,2)+ &
               this%geometry%dsdx%interior(i,nn,k,iEl,1,d,2))* &
               Fphys(iVar,d)
        enddo
        this%twoPointFlux%interior(nn,i,j,k,iEl,iVar,2) = Fc
      enddo

      ! xi^3: pair (i,j,k)-(i,j,nn), project onto avg(Ja^3)
      sR = this%solution%interior(i,j,nn,iEl,1:this%nvar)
      Fphys = this%twopointflux3d(sL,sR)
      do iVar = 1,this%nvar
        Fc = 0.0_prec
        do d = 1,3
          Fc = Fc+0.5_prec*( &
               this%geometry%dsdx%interior(i,j,k,iEl,1,d,3)+ &
               this%geometry%dsdx%interior(i,j,nn,iEl,1,d,3))* &
               Fphys(iVar,d)
        enddo
        this%twoPointFlux%interior(nn,i,j,k,iEl,iVar,3) = Fc
      enddo

    enddo

  endsubroutine TwoPointFluxMethod_ECDGModel3D_t

  subroutine CalculateTendency_ECDGModel3D_t(this)
    implicit none
    class(ECDGModel3D_t),intent(inout) :: this
    ! Local
    integer :: i,j,k,iEl,iVar
    real(prec) :: bMi1,bMi2,bMj1,bMj2,bMk1,bMk2,qwi,qwj,qwk,jac

    call this%solution%BoundaryInterp()
    call this%solution%SideExchange(this%mesh)

    call this%PreTendency()
    call this%SetBoundaryCondition()

    if(this%gradient_enabled) then
      call this%CalculateSolutionGradient()
      call this%SetGradientBoundaryCondition()
      call this%solutionGradient%AverageSides()
    endif

    call this%SourceMethod()
    call this%BoundaryFlux()

    call this%TwoPointFluxMethod()
    call this%twoPointFlux%UpdateDevice()
    call this%twoPointFlux%MappedDivergence(this%fluxDivergence%interior)

    ! Add (1/J) * M^{-1} B^T f_Riemann
    ! 3D boundary side ordering: 1=Bottom, 2=South, 3=East, 4=North, 5=West, 6=Top
    do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
                  k=1:this%solution%N+1,iEl=1:this%mesh%nElem, &
                  iVar=1:this%solution%nVar)

      bMi1 = this%solution%interp%bMatrix(i,2)
      bMi2 = this%solution%interp%bMatrix(i,1)
      bMj1 = this%solution%interp%bMatrix(j,2)
      bMj2 = this%solution%interp%bMatrix(j,1)
      bMk1 = this%solution%interp%bMatrix(k,2)
      bMk2 = this%solution%interp%bMatrix(k,1)
      qwi = this%solution%interp%qWeights(i)
      qwj = this%solution%interp%qWeights(j)
      qwk = this%solution%interp%qWeights(k)
      jac = this%geometry%J%interior(i,j,k,iEl,1)

      this%fluxDivergence%interior(i,j,k,iEl,iVar) = &
        this%fluxDivergence%interior(i,j,k,iEl,iVar)+ &
        (bMk1*this%flux%boundaryNormal(i,j,6,iEl,iVar)+ & ! top
         bMk2*this%flux%boundaryNormal(i,j,1,iEl,iVar))/(qwk*jac)+ & ! bottom
        (bMi1*this%flux%boundaryNormal(j,k,3,iEl,iVar)+ & ! east
         bMi2*this%flux%boundaryNormal(j,k,5,iEl,iVar))/(qwi*jac)+ & ! west
        (bMj1*this%flux%boundaryNormal(i,k,4,iEl,iVar)+ & ! north
         bMj2*this%flux%boundaryNormal(i,k,2,iEl,iVar))/(qwj*jac) ! south

    enddo

    do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
                  k=1:this%solution%N+1,iEl=1:this%mesh%nElem, &
                  iVar=1:this%solution%nVar)

      this%dSdt%interior(i,j,k,iEl,iVar) = &
        this%source%interior(i,j,k,iEl,iVar)- &
        this%fluxDivergence%interior(i,j,k,iEl,iVar)

    enddo

  endsubroutine CalculateTendency_ECDGModel3D_t

endmodule SELF_ECDGModel3D_t
