! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
!
! 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_ECDGModel2D_t
  !! Entropy-conserving DGSEM base class for 2-D conservation laws.
  !!
  !! Extends DGModel2D_t by replacing the standard volume flux divergence with
  !! the split-form (two-point) EC volume term following
  !!   Gassner, Winters, Kopriva (2016)
  !!   Winters, Kopriva, Gassner, Hindenlang (2021)
  !!
  !! Semidiscretization:
  !!   du/dt = source - (1/J) * [ 2 sum_n D_{n,i} F~^1(n,i,j)
  !!                             + 2 sum_n D_{n,j} F~^2(n,i,j)
  !!                             + M^{-1} B^T f_Riemann ]
  !!
  !! where F~^r is the Jacobian-weighted contravariant EC two-point flux in the
  !! r-th computational direction (computed by TwoPointFluxMethod via twopointflux2d),
  !! and f_Riemann is a standard Riemann solver flux at element faces (computed
  !! by BoundaryFlux via riemannflux2d).

  use SELF_DGModel2D_t
  use SELF_MappedTwoPointVector_2D

  implicit none

  type,extends(DGModel2D_t) :: ECDGModel2D_t

    type(MappedTwoPointVector2D) :: twoPointFlux

  contains

    procedure :: Init => Init_ECDGModel2D_t
    procedure :: Free => Free_ECDGModel2D_t

    procedure :: twopointflux2d => twopointflux2d_ECDGModel2D_t
    procedure :: TwoPointFluxMethod => TwoPointFluxMethod_ECDGModel2D_t
    procedure :: CalculateTendency => CalculateTendency_ECDGModel2D_t

  endtype ECDGModel2D_t

contains

  subroutine Init_ECDGModel2D_t(this,mesh,geometry)
    implicit none
    class(ECDGModel2D_t),intent(out) :: this
    type(Mesh2D),intent(in),target :: mesh
    type(SEMQuad),intent(in),target :: geometry

    ! Initialise all parent fields (solution, flux, source, ...)
    call Init_DGModel2D_t(this,mesh,geometry)

    ! Additional two-point flux field
    call this%twoPointFlux%Init(geometry%x%interp,this%nvar,mesh%nElem)
    call this%twoPointFlux%AssociateGeometry(geometry)

  endsubroutine Init_ECDGModel2D_t

  subroutine Free_ECDGModel2D_t(this)
    implicit none
    class(ECDGModel2D_t),intent(inout) :: this

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

  endsubroutine Free_ECDGModel2D_t

  pure function twopointflux2d_ECDGModel2D_t(this,sL,sR) result(flux)
    !! Entropy-conserving two-point flux function.
    !!
    !! Returns the physical-space two-point EC flux for the state pair (sL, sR).
    !! flux(ivar, d) is the d-th physical component (d=1: x, d=2: y) for variable ivar.
    !!
    !! Override this procedure in a concrete model to supply an EC flux.
    !! The stub returns zero (no flux), which is safe but trivial.
    implicit none
    class(ECDGModel2D_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:2)

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

  endfunction twopointflux2d_ECDGModel2D_t

  subroutine TwoPointFluxMethod_ECDGModel2D_t(this)
    !! Computes pre-projected SCALAR contravariant two-point fluxes for all
    !! node pairs and stores them in twoPointFlux%interior(n,i,j,iel,ivar,r).
    !!
    !! Following Trixi.jl for curved meshes, each computational direction r
    !! uses the correct partner node AND the correct averaged metric Ja^r:
    !!
    !!   r=1: pair (i,j)-(nn,j), Ja^1_avg = 0.5*(Ja^1(i,j) + Ja^1(nn,j))
    !!        interior(nn,i,j,...,1) = sum_d Ja^1_avg(d) * F_EC_d(sL, sR)
    !!
    !!   r=2: pair (i,j)-(i,nn), Ja^2_avg = 0.5*(Ja^2(i,j) + Ja^2(i,nn))
    !!        interior(nn,i,j,...,2) = sum_d Ja^2_avg(d) * F_EC_d(sL, sR)
    !!
    !! The result is a SCALAR per (nn,i,j,iel,ivar,r) — no cross-contamination
    !! between directions.
    implicit none
    class(ECDGModel2D_t),intent(inout) :: this
    ! Local
    integer :: nn,i,j,d,iEl,iVar
    real(prec) :: sL(1:this%nvar),sR(1:this%nvar)
    real(prec) :: Fphys(1:this%nvar,1:2)
    real(prec) :: Fc

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

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

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

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

    enddo

  endsubroutine TwoPointFluxMethod_ECDGModel2D_t

  subroutine CalculateTendency_ECDGModel2D_t(this)
    !! Computes du/dt = source - EC-DG flux divergence.
    !!
    !! Volume term : MappedDivergence of twoPointFlux (EC split-form sum, /J)
    !! Surface term: (1/J) * M^{-1} B^T f_Riemann  (weak-form Riemann flux)
    implicit none
    class(ECDGModel2D_t),intent(inout) :: this
    ! Local
    integer :: i,j,iEl,iVar
    real(prec) :: bM1,bM2,qwi,qwj,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() ! fills flux%boundaryNormal with Riemann fluxes

    ! EC volume divergence (includes J weighting via MappedDivergence)
    call this%TwoPointFluxMethod()
    call this%twoPointFlux%UpdateDevice()
    call this%twoPointFlux%MappedDivergence(this%fluxDivergence%interior)

    ! Add weak-form DG surface term: (1/J) * M^{-1} B^T f_Riemann
    ! Boundary side ordering: 1=South, 2=East, 3=North, 4=West
    do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
                  iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar)

      bM1 = this%solution%interp%bMatrix(i,2) ! right/north boundary value at node i
      bM2 = this%solution%interp%bMatrix(i,1) ! left/south boundary value at node i
      qwi = this%solution%interp%qWeights(i)
      qwj = this%solution%interp%qWeights(j)
      jac = this%geometry%J%interior(i,j,iEl,1)

      ! East (side 2) and West (side 4) contributions — vary in xi^1 (i-direction)
      this%fluxDivergence%interior(i,j,iEl,iVar) = &
        this%fluxDivergence%interior(i,j,iEl,iVar)+ &
        (bM1*this%flux%boundaryNormal(j,2,iEl,iVar)+ &
         bM2*this%flux%boundaryNormal(j,4,iEl,iVar))/(qwi*jac)+ &
        (this%solution%interp%bMatrix(j,2)*this%flux%boundaryNormal(i,3,iEl,iVar)+ &
         this%solution%interp%bMatrix(j,1)*this%flux%boundaryNormal(i,1,iEl,iVar))/(qwj*jac)

    enddo

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

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

    enddo

  endsubroutine CalculateTendency_ECDGModel2D_t

endmodule SELF_ECDGModel2D_t