! Copyright (c) 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2019 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2021 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@dlr.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.
!> This module contains functions for MRT relaxation paramater for different
!! stencil layouts.
!! NOTE: The order of relaxation entries is consistent with moments
!! transformation matrix used in compute kernel.
!! author: Kannan Masilamani
! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2011 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2011 Konstantin Kleinheinz <k.kleinheinz@grs-sim.de>
! Copyright (c) 2011-2012 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2012 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2013-2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@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.
module mus_mrtRelaxation_module
  use env_module,       only: rk
  use tem_aux_module,   only: tem_abort
  use tem_float_module, only: operator(.fne.)

  use mus_moments_type_module,  only: mus_moment_type
  use mus_scheme_header_module,      only: mus_scheme_header_type

  implicit none
  private

  public :: mus_assign_mrt_ptr
  public :: mus_proc_mrt

  abstract interface
    !> function pointers to obtain relaxation matrix for MRT collision operator
    pure function mus_proc_mrt(omegaKine, omegaBulk, QQ) result (s_mrt)
      import :: rk

      !> omega related to kinematic viscosity
      real(kind=rk), intent(in) :: omegaKine

      !> omega related to bulk viscosity for compressible model
      real(kind=rk), intent(in) :: omegaBulk

      !> number of directions
      integer, intent(in) :: QQ

      !> output: MRT diagonal relaxation parameter
      real(kind=rk) :: s_mrt(QQ)
    end function mus_proc_mrt

  end interface

contains

  ! ************************************************************************** !
  !> This function returns mrt function pointer according to scheme definition.
  !! In Jonas Toelke paper (2006) about MRT, the following notions are used:\n
  !!  s(a) = s(2)
  !!  s(b) = s(3)
  !!  s(c) = s(5) = s(7) = s(9)
  !!  s(d) = s(11) = s(13
  !!  s(e) = s(17) = s(18) = s(19)
  !!  s(w) = s(10) = s(12) = s(14) = s(15) = s(16)
  !! It is suggested that, for D3Q19,
  !!  s(a) = s(b) = s(c) = s(d) = s(e) = max( s(w), -1.0 )
  !! Notice that the collision matrix S used in this papar corresponds to
  !! -omega in BGK model, because it express the LB equation is slightly
  !! different way.
  !!
  subroutine mus_assign_mrt_ptr(mrtPtr, schemeHeader)
    ! --------------------------------------------------------------------------
    !> mrt function pointer
    procedure(mus_proc_mrt), pointer :: mrtPtr
    !> Scheme header information
    type(mus_scheme_header_type), intent(in) :: schemeHeader
    ! --------------------------------------------------------------------------
    mrtPtr => null()
    select case (trim(schemeHeader%relaxation))
    case ( 'mrt', 'mrt_generic' )
      select case (trim(schemeHeader%layout))
      case ('d2q9')
        if (trim(schemeHeader%kind) == 'fluid_incompressible') then
          mrtPtr => mrt_d2q9_incomp
        else
          mrtPtr => mrt_d2q9
        end if
      case ('d3q15')
        mrtPtr => mrt_d3q15
      case ('d3q19')
        mrtPtr => mrt_d3q19
      case ('d3q27')
        mrtPtr => mrt_d3q27
      case default
        call tem_abort('Error: Unknown layout for mrt relaxation')
      end select
    case default
      ! set all relaxation paramter to same omega value as default
      mrtPtr => mrt_bgk
    end select
  end subroutine mus_assign_mrt_ptr
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This function returns mrt relaxation diagonal matrix for d2q9 layout
  !! Parameters are taken from:
  !! Lallemand, P., & Luo, L. (2000).
  !! Theory of the lattice boltzmann method: dispersion, dissipation,
  !! isotropy, galilean invariance, and stability. Physical Review. E,
  !! Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary
  !! Topics, 61(6 Pt A), 6546–62.
  !!
  !! Another paper for d2q9 with different set of parameters
  !! Chen, S., Peng, C., Teng, Y., Wang, L. P., & Zhang, K. (2016).
  !! Improving lattice Boltzmann simulation of moving particles in a
  !! viscosity flow using local grid refinement. Computers and Fluids,
  !! 136, 228–246.
  pure function mrt_d2q9_incomp(omegaKine, omegaBulk, QQ) result(s_mrt)
    ! --------------------------------------------------------------------------
    !> omega related to kinematic viscosity
    real(kind=rk), intent(in) :: omegaKine
    !> omega related to bulk viscosity
    real(kind=rk), intent(in) :: omegaBulk
    !> number of directions
    integer, intent(in) :: QQ
    !> output mrt diagonal matrix
    real(kind=rk) :: s_mrt(QQ)
    ! --------------------------------------------------------------------------
    s_mrt = 0.0_rk
    ! set relaxation by bulk viscostiy
    ! set to 1.63 for both compressible and incompressible model
    ! KM: Tried omegaBulk for compressible model and it is unstable
    s_mrt(2)   = 1.63_rk
    s_mrt(3)   = 1.14_rk
    s_mrt(5)   = 1.92_rk
    s_mrt(7)   = 1.92_rk
    s_mrt(8:9) = omegaKine

  end function mrt_d2q9_incomp
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This function returns mrt relaxation diagonal matrix for d2q9 layout
  !! Parameters are taken from:
  !! Lallemand, P., & Luo, L. (2000).
  !! Theory of the lattice boltzmann method: dispersion, dissipation,
  !! isotropy, galilean invariance, and stability. Physical Review. E,
  !! Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary
  !! Topics, 61(6 Pt A), 6546–62.
  !!
  !! Another paper for d2q9 with different set of parameters
  !! Chen, S., Peng, C., Teng, Y., Wang, L. P., & Zhang, K. (2016).
  !! Improving lattice Boltzmann simulation of moving particles in a
  !! viscosity flow using local grid refinement. Computers and Fluids,
  !! 136, 228–246.
  ! modified by GGSpinelli. All V^n moments are relaxed with omegaBulk
  pure function mrt_d2q9(omegaKine, omegaBulk, QQ) result(s_mrt)
    ! --------------------------------------------------------------------------
    !> omega related to kinematic viscosity
    real(kind=rk), intent(in) :: omegaKine
    !> omega related to bulk viscosity
    real(kind=rk), intent(in) :: omegaBulk
    !> number of directions
    integer, intent(in) :: QQ
    !> output mrt diagonal matrix
    real(kind=rk) :: s_mrt(QQ)
    ! --------------------------------------------------------------------------
    s_mrt = 0.0_rk
    ! set relaxation by bulk viscostiy
    ! set to 1.63 for both compressible and incompressible model
    ! KM: Tried omegaBulk for compressible model and it is unstable
    ! GGS: omega2 = omegaBulk & omega3 = 1.14 => not stable!
    ! GGS: omega2 = omega3 = omegaBulk => works!
    s_mrt(2:3) = omegaBulk
    s_mrt(5)   = 1.92_rk
    s_mrt(7)   = 1.92_rk
    s_mrt(8:9) = omegaKine

  end function mrt_d2q9
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This function returns mrt relaxation diagonal matrix for d3q15 layout
  !! Parameters are taken from:
  !! D’Humières, D., Ginzburg, I., Krafczyk, M., Lallemand, P., & Luo, L.-S.
  !! (2002). Multiple-relaxation-time lattice Boltzmann models in three
  !! dimensions. Philosophical Transactions. Series A, Mathematical,
  !! Physical, and Engineering Sciences, 360(1792), 437–51.
  pure function mrt_d3q15(omegaKine, omegaBulk, QQ) result(s_mrt)
    ! --------------------------------------------------------------------------
    !> omega related to kinematic viscosity
    real(kind=rk), intent(in) :: omegaKine
    !> omega related to bulk viscosity
    real(kind=rk), intent(in) :: omegaBulk
    !> number of directions
    integer, intent(in) :: QQ
    !> output mrt diagonal matrix
    real(kind=rk) :: s_mrt(QQ)
    ! --------------------------------------------------------------------------
    s_mrt = 0.0_rk
    s_mrt(2)     = omegaBulk
    s_mrt(3)     = 1.20_rk
    s_mrt(5)     = 1.60_rk
    s_mrt(7)     = 1.60_rk
    s_mrt(9)     = 1.60_rk
    s_mrt(10:14) = omegaKine
    s_mrt(15)    = 1.20_rk

  end function mrt_d3q15
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This function returns mrt relaxation diagonal matrix for d3q19 layout
  !! Parameters are taken from:
  !! D’Humières, D., Ginzburg, I., Krafczyk, M., Lallemand, P., & Luo, L.-S.
  !! (2002). Multiple-relaxation-time lattice Boltzmann models in three
  !! dimensions. Philosophical Transactions. Series A, Mathematical,
  !! Physical, and Engineering Sciences, 360(1792), 437–51.
  pure function mrt_d3q19(omegaKine, omegaBulk, QQ) result(s_mrt)
    ! --------------------------------------------------------------------------
    !> omega related to kinematic viscosity
    real(kind=rk), intent(in) :: omegaKine
    !> omega related to bulk viscosity
    real(kind=rk), intent(in) :: omegaBulk
    !> number of directions
    integer, intent(in) :: QQ
    !> output mrt diagonal matrix
    real(kind=rk) :: s_mrt(QQ)
    ! --------------------------------------------------------------------------
    s_mrt = 0.0_rk
    s_mrt(2)     = omegaBulk
    s_mrt(3)     = 1.40_rk
    s_mrt(5)     = 1.20_rk
    s_mrt(7)     = 1.20_rk
    s_mrt(9)     = 1.20_rk
    s_mrt(10)    = omegaKine
    s_mrt(11)    = 1.40_rk
    s_mrt(12)    = omegaKine
    s_mrt(13)    = 1.40_rk
    s_mrt(14:16) = omegaKine
    s_mrt(17:19) = 1.98_rk

  end function mrt_d3q19
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This function returns mrt relaxation diagonal matrix for d3q27 layout
  pure function mrt_d3q27(omegaKine, omegaBulk, QQ) result(s_mrt)
    ! --------------------------------------------------------------------------
    !> omega related to kinematic viscosity
    real(kind=rk), intent(in) :: omegaKine
    !> omega related to bulk viscosity
    real(kind=rk), intent(in) :: omegaBulk
    !> number of directions
    integer, intent(in) :: QQ
    !> output mrt diagonal matrix
    real(kind=rk) :: s_mrt(QQ)
    ! --------------------------------------------------------------------------
    s_mrt = 0.0_rk
    s_mrt(5:9)   = omegaKine
    s_mrt(10)    = omegaBulk
    s_mrt(11:13) = 1.50_rk
    s_mrt(14:17) = 1.74_rk
    s_mrt(18)    = 1.4_rk
    s_mrt(19:23) = 1.98_rk
    s_mrt(24:26) = 1.83_rk
    s_mrt(27)    = 1.61_rk

  end function mrt_d3q27
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> set all relaxation parameter to same omega, results in bgk collision
  pure function mrt_bgk(omegaKine, omegaBulk, QQ) result(s_mrt)
    ! --------------------------------------------------------------------------
    !> omega related to kinematic viscosity
    real(kind=rk), intent(in) :: omegaKine
    !> omega related to bulk viscosity
    real(kind=rk), intent(in) :: omegaBulk
    !> number of directions
    integer, intent(in) :: QQ
    !> output mrt diagonal matrix
    real(kind=rk) :: s_mrt(QQ)
    ! --------------------------------------------------------------------------
    s_mrt(:) = omegaKine
  end function mrt_bgk
  ! ************************************************************************** !

end module mus_mrtRelaxation_module
