! Copyright (c) 2016-2021 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2019 Seyfettin Bilgi <seyfettin.bilgi@student.uni-siegen.de>
!
! 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.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “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 UNIVERSITY OF SIEGEN 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.
! ***************************************************************************** !
!> author: Kannan Masilamani
!! Module containing subroutines for building MUSUBI specific source 
!! variables
!!
module mus_source_var_module
  use, intrinsic :: iso_c_binding, only: c_ptr, c_f_pointer

  ! include treelm modules
  use env_module,               only: rk, long_k, labelLen, solSpecLen
  use tem_tools_module,         only: upper_to_lower
  use tem_aux_module,           only: tem_abort
  use tem_varSys_module,        only: tem_varSys_type,               &
    &                                 tem_varSys_append_derVar,      &
    &                                 tem_varSys_proc_point,         &
    &                                 tem_varSys_proc_element,       &
    &                                 tem_varSys_proc_setParams,     &
    &                                 tem_varSys_proc_getParams,     &
    &                                 tem_varSys_proc_setupIndices,  &
    &                                 tem_varSys_proc_getValOfIndex, &
    &                                 tem_varSys_dump
  use tem_varMap_module,        only: tem_possible_variable_type, &
    &                                 init, append, truncate,     &
    &                                 tem_variable_loadMapping
  use tem_stringKeyValuePair_module, only: init, truncate, &
    &                                      grw_stringKeyValuePairArray_type
  use tem_stencil_module,       only: tem_stencilHeader_type
  use tem_time_module,          only: tem_time_type
  use treelmesh_module,         only: treelmesh_type
  use tem_logging_module,       only: logUnit

  ! include musubi modules
  use mus_scheme_header_module, only: mus_scheme_header_type
  use mus_physics_module,       only: mus_convertFac_type
  use mus_derVarPos_module,     only: mus_derVarPos_type

  ! include aotus modules
  use aotus_module,     only: flu_State, aot_get_val, aoterr_Fatal,            &
    &                         aoterr_NonExistent, aoterr_WrongType
  use aot_table_module, only: aot_table_open, aot_table_close, aot_get_val

  implicit none
  private

  public :: mus_source_type
  public :: mus_source_op_type
  public :: mus_create_poss_srcVar
  public :: mus_load_source_var
  public :: mus_source_cleanup


  ! ************************************************************************** !
  !! \todo KM: 20210301 Move absorbLayer to seperate module.
  !! Make target_pressure and target_velocity as stFun.
  !> Contains additional information for absorblayer source
  type absorbLayer_config_type
    !> target pressure
    real(kind=rk) :: target_pressure
    !> target velocityX, velocityY and velocityZ
    real(kind=rk) :: target_velocity(3)
    !> Use time average for pressure. Default: false.
    logical :: isPressDyn = .false.
    !> Use time average for Velocity. Default: false.
    logical :: isVelDyn = .false.
    !> Number of iterations to record for time-averaging
    integer :: nRecord = 100
  end type absorbLayer_config_type
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Stores time average values of density and velocity for dynamic absorbLayer
  type absorbLayerDynamic_average_type
    !> density time average
    real(kind=rk), allocatable :: dens(:)
    !> velocity time average
    real(kind=rk), allocatable :: velX(:)
    real(kind=rk), allocatable :: velY(:)
    real(kind=rk), allocatable :: velZ(:)
    !> It is used to initialiye dynamic average density with initial condition
    logical :: isInitDens = .true.

    !> It is used to initialiye dynamic average velocity with initial condition
    logical :: isInitVel = .true.
  end type absorbLayerDynamic_average_type
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Contains information for absorblayer
  type mus_absorbLayer_type
    !> Information loaded from configuration file
    type(absorbLayer_config_type) :: config

    !> Smoothing factor for expoential moving average
    !! = 2 / (nRecord+1)
    real(kind=rk) :: smoothFac
  end type mus_absorbLayer_type
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Contains source elements position in state array and idx to access 
  !! data variable refered in config file.
  !! This type is defined for each level
  type mus_source_elems_type
    !> Number of source elements on this level.
    !! nFluids + nGhosts
    integer :: nElems

    !> Position of elements in state array to apply source terms.
    !! Position in state array is same as position in total list
    !! Size: nElems
    integer, allocatable :: posInTotal(:)

    !> Index to access point data type to retrieve values from variable 
    !! refered for source variable
    integer, allocatable :: idx(:)

    !> source field value obtained from ST_fun data variable.
    !! Filled only for elements where source is active i.e. elements in
    !! posInTotal. 
    !! size: nElems*nComponents
    !! \todo KM: might be not neccessary
!KM!    real(kind=rk), allocatable :: val(:)

    !> Contains time average values of density and velocity for dynamic 
    !! absorblayer.
    !! \todo KM: 02042021 Introduce method_data c_ptr and point to 
    !! dynAvg for absorbLayer and change intent(inout) to intent(in) in
    !! proc_addSrcToAuxField.
    type(absorbLayerDynamic_average_type) :: dynAvg
  end type mus_source_elems_type
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Description contains list of elements on which source is active and
  !! function pointer to update source
  type mus_source_op_type
    !> Position of this source term variable in the varSys
    integer :: srcTerm_varPos

    !> Position of data variable provided in config file in the varSys
    integer :: data_varPos

    !> Contains source elements position in state array/total list for 
    !! each level
    type(mus_source_elems_type), allocatable :: elemLvl(:)

    !> Function to update state with source term
    procedure(proc_apply_source), pointer :: applySrc => null()

    !> to use source field array
    !KM!logical :: useSrcFieldVal = .false.

    !> name of the source variable
    character(len=labelLen) :: varname

    !> Function pointer to append source field to auxilary variable
    procedure(proc_addSrcToAuxField), pointer :: addSrcToAuxField => null()

    !> Order of approximation for source like force, electric_field, 
    !! charge_density.
    !! Order = 1, uses force term in BE approximated by forward Euler method
    !! Order = 2, uses force term in BE approximated by Trapezoidal method.
    !! For order 2, macroscopic source is also added to auxField.
    !! For fluid, fluid_incompressible, multispecies_liquid: source is added
    !! to momentum and for poisson: source is added to potential.
    !! Default: order = 2.
    integer :: order


    !> Additional config information for absorbLayer
    type(mus_absorbLayer_type) :: absLayer
  end type mus_source_op_type
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> Description of musubi source type 
  type mus_source_type
    !> Contains source elements position in tree%treeID and
    !! function pointer to update source
    !! Size: varDict%nVals
    type(mus_source_op_type), allocatable :: method(:)

    !> Dictionary of source variable with 
    !! varDict%val()%key is the name of source variable and
    !! varDict%val()%value is the name of variable provided for the key
    type(grw_stringKeyValuePairArray_type) :: varDict
  end type mus_source_type
  ! *************************************************************************** !


  ! *************************************************************************** !
  abstract interface
    !> Abstract interface to update state with source terms
    subroutine proc_apply_source( fun, inState, outState, neigh, auxField, &
      & nPdfSize, iLevel, varSys, time, phyConvFac, derVarPos )
      import :: rk, mus_source_op_type, tem_varSys_type, tem_time_type, &
        &       mus_convertFac_type, mus_derVarPos_type

      !> Description of method to update source
      class(mus_source_op_type), intent(in) :: fun

      !> input  pdf vector
      !! \todo KM: instate is passed to compute auxField.
      !! Since auxField is precomputed from instate and passed to this routine.
      !! instate can be removed
      real(kind=rk), intent(in)  ::  inState(:)
  
      !> output pdf vector
      real(kind=rk), intent(inout) :: outState(:)
  
      !> connectivity Array corresponding to state vector
      integer,intent(in) :: neigh(:)

      !> auxField array
      real(kind=rk), intent(in) :: auxField(:)

      !> number of elements in state Array
      integer, intent(in) :: nPdfSize

      !> current level
      integer, intent(in) :: iLevel

      !> variable system
      type(tem_varSys_type), intent(in) :: varSys

      !> Point in time at which to evaluate the variable.
      type(tem_time_type), intent(in)  :: time

      !> Physics conversion factor for current level
      type(mus_convertFac_type), intent(in) :: phyConvFac

      !> position of derived quantities in varsys
      type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    end subroutine proc_apply_source

    !> Interface to add source to auxField vars if addSrcToAuxField is true 
    !! in source_op_type for all nSolve elements (nFluids+nGhostFromCoarser
    !! +nGhostFromFiner). Halo elements are exchanged
    subroutine proc_addSrcToAuxField(fun, auxField, iLevel, time, varSys, &
      & phyConvFac, derVarPos)
      import :: rk, tem_varSys_type, mus_derVarPos_type, mus_convertFac_type, &
        &       tem_time_type, mus_source_op_type

      !> Description of method to update source
      class(mus_source_op_type), intent(inout) :: fun
      !> output auxField array
      real(kind=rk), intent(inout)          :: auxField(:)
      !> current level
      integer, intent(in)                   :: iLevel
      !> current timing information
      type(tem_time_type), intent(in)       :: time
      !> variable system definition
      type(tem_varSys_type), intent(in)     :: varSys
      !> Physics conversion factor for current level
      type(mus_convertFac_type), intent(in) :: phyConvFac
      !> position of derived quantities in varsys
      type(mus_derVarPos_type), intent(in)  :: derVarPos(:)
    end subroutine proc_addSrcToAuxField 

  end interface
  ! *************************************************************************** !

contains


  ! *************************************************************************** !
  !> Routine initialize possible source variable depends on scheme kind
  subroutine mus_create_poss_srcVar(poss_srcVar, schemeHeader)
    ! --------------------------------------------------------------------------!
    !> possible source variables
    type(tem_possible_variable_type), intent(out) :: poss_srcVar

    !> Identifier of the scheme
    type(mus_scheme_header_type), intent(in) :: schemeHeader
    ! --------------------------------------------------------------------------!
    write(logUnit(10),*) 'Creating possible source terms '
    call init(me = poss_srcVar, length = 2 )

    select case(trim(schemeHeader%kind))
    case ('fluid', 'fluid_incompressible')
      ! body force
      call append(me          = poss_srcVar, &
        &         varName     = 'force',     &
        &         nComponents = 3            )
      ! absorb layer, STfun should return sponge_strength
      ! Target pressure and velocity are defined under absorb_layer_target
      call append(me          = poss_srcVar,    &
        &         varName     = 'absorb_layer', &
        &         nComponents = 1               )

      ! Absorn layer for inlet boundary with dynamic pressure and
      ! constant velocity.
      call append(me          = poss_srcVar,          &
        &         varName     = 'absorb_layer_inlet', &
        &         nComponents = 1                     )
      ! Absorn layer for outlet boundary with constant pressure and
      ! dynamic velocity.
      call append(me          = poss_srcVar,           &
        &         varName     = 'absorb_layer_outlet', &
        &         nComponents = 1                     )

    case ('passive_scalar')
      call append(me          = poss_srcVar,       &
        &         varName     = 'equal_injection', &
        &         nComponents = 1                  )

      call append(me          = poss_srcVar, &
        &         varName     = 'injection', &
        &         nComponents = 1            )

    case ('nernst_planck')
      call append(me          = poss_srcVar,      &
        &         varName     = 'electric_field', &
        &         nComponents = 3                 )

    case ('multispecies_liquid')
      call append(me          = poss_srcVar,      &
        &         varName     = 'electric_field', &
        &         nComponents = 3                 )

      call append(me          = poss_srcVar, &
        &         varName     = 'force',     &
        &         nComponents = 3            )

    case ('poisson')
      call append(me          = poss_srcVar,      &
        &         varName     = 'charge_density', &
        &         nComponents = 1                 )

    case default
      write(logUnit(1),*) 'No possible source term defined for scheme kind:' &
        &                 //trim(schemeHeader%kind)
    end select
    call truncate(poss_srcVar)

  end subroutine mus_create_poss_srcVar
  ! *************************************************************************** !


  ! ***************************************************************************!
  !> Routine load musubi source terms for given key.
  !! key is glob_source or source
  subroutine mus_load_source_var(me, possVars, conf, parent, key, varSys)
    ! --------------------------------------------------------------------------!
    !> Source variable type to initialize
    type(mus_source_type), intent(out) :: me
    !> possible source variables
    type(tem_possible_variable_type), intent(in) :: possVars
    !> flu state
    type( flu_State ) :: conf
    !> parent handle if scheme table is defined
    integer, intent(in), optional :: parent
    !> key to load source term
    character(len=*), intent(in) :: key
    !> Global variable system
    type(tem_varSys_type), intent(inout) :: varSys
    ! --------------------------------------------------------------------------!
    integer :: iSrc, iError, srchandle
    character(len=labelLen) :: varname_order
    ! --------------------------------------------------------------------------!
    ! initialize growing array stringKeyValuePair
    call init( me = me%varDict )

    ! load the global source variables
    call tem_variable_loadMapping( possVars = possVars,   &
      &                            conf     = conf,       &
      &                            parent   = parent,     &
      &                            key      = key,        &
      &                            varDict  = me%varDict, &
      &                            varSys   = varSys      )

    ! truncate varDict
    call truncate( me = me%varDict )

    ! If source is defined then reopen source table to load order of
    ! approximatation to use for defined source variable
    if (me%varDict%nVals > 0) then
      ! allocate source method
      allocate(me%method(me%varDict%nVals))

      call aot_table_open( L       = conf,      &
        &                  parent  = parent,    &
        &                  thandle = srchandle, &
        &                  key     = key        )

      do iSrc = 1, me%varDict%nVals
         varname_order = trim(me%varDict%val(iSrc)%key)//'_order'
         call aot_get_val( L       = conf,                  &
           &               thandle = srchandle,             &
           &               key     = varname_order,         &
           &               val     = me%method(iSrc)%order, &
           &               default = 2,                     &
           &               ErrCode = iError                 )

         if (btest(iError, aoterr_Fatal)) then
           write(logUnit(1),*)'FATAL Error occured, while retrieving source' &
             &                //'order'
           if (btest(iError, aoterr_WrongType)) then
             write(logUnit(1),*)'Variable has wrong type!'
             write(logUnit(1),*)'STOPPING'
             call tem_abort()
           endif
         end if
         write(logUnit(1),'(A,I0)') ' Order of approximation for '          &
           &                        //trim(me%varDict%val(iSrc)%key)//': ', &
           &                        me%method(iSrc)%order

         ! Load additional information for absorblayer
         select case (trim(me%varDict%val(iSrc)%key)) 
         case ('absorb_layer')
           call load_absorbLayer( me       = me%method(iSrc)%absLayer%config, &
             &                    conf     = conf,                            &
             &                    key      = 'absorb_layer_target',           & 
             &                    parent   = srchandle,                       &
             &                    loadPres = .true.,                          &
             &                    loadVel  = .true.                           )
         case ('absorb_layer_inlet')
           call load_absorbLayer( me       = me%method(iSrc)%absLayer%config, &
             &                    conf     = conf,                            &
             &                    key      = 'absorb_layer_inlet_target',     & 
             &                    parent   = srchandle,                       &
             &                    loadPres = .false.,                         &
             &                    loadVel  = .true.                           )
         case ('absorb_layer_outlet')
           call load_absorbLayer( me       = me%method(iSrc)%absLayer%config, &
             &                    conf     = conf,                            &
             &                    key      = 'absorb_layer_outlet_target',    & 
             &                    parent   = srchandle,                       &
             &                    loadPres = .true.,                          &
             &                    loadVel  = .false.                          )
         end select
      end do
    else
      allocate(me%method(0))
    end if  

  end subroutine mus_load_source_var
  ! ***************************************************************************!

  ! ************************************************************************** !
  !> This routine act as a destructor for source type.
  !! The arrays allocated in mus_init_sourceTerms are destroyed here
  subroutine mus_source_cleanup(me)
    ! --------------------------------------------------------------------------
    type(mus_source_type), intent(inout)   :: me
    ! --------------------------------------------------------------------------
    integer :: iSrc
    ! --------------------------------------------------------------------------
    ! KM: DO NOT DESTROY VARDICT IN SOURCE_TYPE AS IT CONTAINS CONFIG INFO
    do iSrc = 1, me%varDict%nVals 
      deallocate(me%method(iSrc)%elemLvl)
    end do  
  end subroutine mus_source_cleanup
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine load additional information for absorblayer 
  subroutine load_absorbLayer(me, conf, key, parent, loadPres, loadVel)
    ! -------------------------------------------------------------------------!
    !> Absorb layer 
    type(absorbLayer_config_type), intent(out) :: me
    !> flu state
    type( flu_State ) :: conf
    !> Table name to load target states
    character(len=*), intent(in) :: key
    !> parent source handle
    integer, intent(in) :: parent
    !> Load pressure if true else set to dynamic
    logical, intent(in) :: loadPres
    !> Load velocity if true else set to dynamic
    logical, intent(in) :: loadVel
    ! -------------------------------------------------------------------------!
    integer :: target_handle, iError, vError(3), errFatal(3)
    character(len=labelLen) :: tarStateAsStr
    ! -------------------------------------------------------------------------!
    errfatal = aotErr_Fatal
 
    call aot_table_open( L       = conf,                 &
      &                  parent  = parent,               &
      &                  thandle = target_handle,        &
      &                  key     = trim(key)             )

    ! Load target pressure and velocity only for static absorblayer
    write(logUnit(1),*) ' * Target state:'
    ! target state: pressure
    if (loadPres) then
      call aot_get_val( L       = conf,               &
        &               thandle = target_handle,      &
        &               key     = 'pressure',         &
        &               val     = me%target_pressure, &
        &               ErrCode = iError              )
      if (btest(iError, aoterr_Fatal)) then
        write(logUnit(1), *) 'Error loading target pressure as value'
        write(logUnit(5), *) 'Trying to load as string'
        call aot_get_val( L       = conf,          &
          &               thandle = target_handle, &
          &               key     = 'pressure',    &
          &               val     = tarStateAsStr, &
          &               ErrCode = iError         )
        if (btest(iError, aoterr_Fatal)) then
          write(logUnit(1),*) 'Error loading target pressure as string'
          call tem_abort()
        end if

        if (upper_to_lower(trim(tarStateAsStr)) == 'dynamic') then
          me%isPressDyn = .true.
          me%target_pressure = 0.0_rk
          write(logUnit(1),*) "    pressure = dynamic'"
        else
          call tem_abort('Target pressure is neither value or string "dynamic"')
        end if
      else
        write(logUnit(1),*) '    pressure =', me%target_pressure
      end if
    else
      me%isPressDyn = .true.
      me%target_pressure = 0.0_rk
      write(logUnit(1),*) "    pressure = 'dynamic'"
    end if

    ! target state: velocity
    if (loadVel) then
      call aot_get_val( L       = conf,               &
        &               thandle = target_handle,      &
        &               key     = 'velocity',         &
        &               val     = me%target_velocity, &
        &               ErrCode = vError              )
      if (any(btest( vError, errFatal )) ) then
        write(logUnit(1), *) 'Error loading target velocity as table'
        write(logUnit(5), *) 'Trying to load as string'
        call aot_get_val( L       = conf,          &
          &               thandle = target_handle, &
          &               key     = 'velocity',    &
          &               val     = tarStateAsStr, &
          &               ErrCode = iError         )
        if (btest(iError, aoterr_Fatal)) then
          write(logUnit(1),*) 'Error loading target velocity as string'
          call tem_abort()
        end if

        if (upper_to_lower(trim(tarStateAsStr)) == 'dynamic') then
          me%isVelDyn = .true.
          me%target_velocity = 0.0_rk
          write(logUnit(1),*) "    velocity = dynamic'"
        else
          call tem_abort('Target velocity is neither value or string "dynamic"')
        end if

      else
        write(logUnit(1),*) '    velocity =', me%target_velocity(1:3)
      end if
    else
      me%isVelDyn = .true.
      me%target_velocity = 0.0_rk
      write(logUnit(1),*) "    velocity = 'dynamic'"
    end if

    ! Load nRecord for time averaging
    if (me%isPressDyn .or. me%isVelDyn) then
      call aot_get_val( L       = conf,          &
        &               thandle = target_handle, &
        &               key     = 'nrecord',     &
        &               val     = me%nRecord,    &
        &               ErrCode = iError         )
      if (btest(iError, aoterr_Fatal)) then
        write(*,*) 'FATAL Error occured, when loading nrecord for absorb layer'
        call tem_abort()
      end if
      write(logUnit(1),*) '    nRecord =', me%nRecord
    end if

    call aot_table_close(conf, target_handle)
  
  end subroutine load_absorbLayer
  ! ************************************************************************** !

end module mus_source_var_module
! **************************************************************************** !
