! See copyright notice in the COPYRIGHT file.

! ****************************************************************************** !
!> author: Jiaxing Qi
!! This module keeps all information about the nonNewtonian models
!!
!! For further information about the theory visit the
!! [non-newtonian theory page](../page/mus_nonNewtonianTheory.html).
!!
module mus_nonNewtonian_module

  ! include treelm modules
  use env_module,           only: rk, labellen
  use tem_tools_module,     only: tem_horizontalSpacer
  use tem_logging_module,   only: logUnit

  ! include aotus modules
  use aotus_module,     only: flu_State, aot_get_val
  use aot_table_module, only: aot_table_open, aot_table_close
  use aot_out_module,   only: aot_out_type, aot_out_val, aot_out_open_table,   &
    &                         aot_out_close_table

  implicit none

  private

  public :: mus_nNwtn_type
  public :: mus_nNwtn_load
  public :: mus_nNwtn_save2lua
  public :: mus_nNwtn_dump2OutUnit

  !> Identifier for Power-Law model
  integer, parameter, public :: mus_nNwtn_PL = 1

  !> Identifier for Casson model
  integer, parameter, public :: mus_nNwtn_CS = 2

  !> Identifier for Carreau-Yasuda model
  integer, parameter, public :: mus_nNwtn_CY = 3

  ! -----------------------------------------------------------------------------
  !> The nonNewtonian power law model parameter
  !!
  !! This date type gathers parameters of a power law (PL) model.
  !! It is encapsulated in mus_nNwtn_type.
  type mus_nNwtn_PL_type

    !> exponentiation parameter
    real(kind=rk) :: n = 0.5_rk

    !> Dynamic viscosity parameter when shear rate equals to 1
    real(kind=rk) :: k = 0.0035_rk

    !> parameter for computation
    real(kind=rk) :: n_1

  end type mus_nNwtn_PL_type
  ! -----------------------------------------------------------------------------

  ! -----------------------------------------------------------------------------
  !> The nonNewtonian power law model parameter
  !!
  !! This date type gathers parameters of the Carreau-Yasuda (CY) model
  !! It is encapsulated in mus_nNwtn_type.
  type mus_nNwtn_CY_type

    !> model parameter
    real(kind=rk) :: n = 0.2128_rk

    !> model parameter
    real(kind=rk) :: a = 0.64_rk

    !> model parameter
    real(kind=rk) :: lambda = 8.2_rk

    !> model parameter
    real(kind=rk) :: visc0 = 0.16_rk

    !> model parameter
    real(kind=rk) :: viscInf = 0.0035_rk

    !> calculated parameter for later usage, n_1_a = (n-1)/a
    real(kind=rk) :: n_1_a = (0.2128_rk - 1._rk) / 0.64_rk

  end type mus_nNwtn_CY_type
  ! -----------------------------------------------------------------------------

  ! -----------------------------------------------------------------------------
  !> The nonNewtonian power law model parameter
  !!
  !! This date type gathers parameters of the Casson model
  !! It is encapsulated in mus_nNwtn_type.
  type mus_nNwtn_CS_type

    !> model parameter
    real(kind=rk) :: k0 = 0.1937_rk

    !> model parameter
    real(kind=rk) :: k1 = 0.055_rk

  end type mus_nNwtn_CS_type
  ! -----------------------------------------------------------------------------

  ! -----------------------------------------------------------------------------
  !> The nonNewtonian fluid feature description type
  !!
  !! This date type gathers related parameters of a nonNewtonian fluid.
  type mus_nNwtn_type
    !> Indicator whether nonNewtonian feature is active
    !! maybe not useful. schemeHeader%kind can used to check if nNwtn is active
    logical :: active = .false.

    !> nonNewtonian fluid model label
    character(len=labellen) :: label

    !> nonNewtonian fluid model identifier
    integer :: model

    !> Power Law (PL) model type
    type( mus_nNwtn_PL_type ) :: PL

    !> Carreau-Yasuda (CY) model type
    type( mus_nNwtn_CY_type ) :: CY

    !> Casson model type
    type( mus_nNwtn_CS_type ) :: CS
  end type mus_nNwtn_type
  ! -----------------------------------------------------------------------------

contains

  ! **************************************************************************** !
  !> Read in the nonNewtonian table from Lua file and dump parameters to logUnit
  !! Specificly, this routine calls each model parameter loader.
  !!
  subroutine mus_nNwtn_load( me, conf, parent )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian type
    type( mus_nNwtn_type ), intent(out) :: me
    !> lua state
    type( flu_state ), intent(inout) :: conf
    !> parent handle
    integer, intent(in), optional  :: parent
    ! ---------------------------------------------------------------------------
    integer :: nonNwt_table
    integer :: iError
    CHARACTER(LEN=12) :: nonNwt_table_str
    ! ---------------------------------------------------------------------------

    nonNwt_table_str = "nonNewtonian"

    ! if nonNewtonian informations in scheme table parentHandle /= 0
    call aot_table_open( L       = conf,            &
      &                  parent  = parent,          &
      &                  thandle = nonNwt_table,    &
      &                  key     = nonNwt_table_str )

    ! when table exist, read in parameters from table
    if ( nonNwt_table /= 0 ) then
      ! Set nonNewtonian feature in on
      me%active = .true.

      ! load model label name
      call aot_get_val(L       = conf,        &
        &              tHandle = nonNwt_table,&
        &              key     = 'model',     &
        &              val     = me%label,    &
        &              default = 'power_law', &
        &              ErrCode = iError       )

      ! load model parameters by calling model loader
      ! set model identifier
      select case( trim(me%label) )
        case ( 'power_law', 'Power_Law', 'PL' )
          me%model = mus_nNwtn_PL
          call mus_nNwtn_PL_load( me%PL, conf, nonNwt_table )
        case ( 'carreau_yasuda', 'Carreau_Yasuda', 'CY' )
          me%model = mus_nNwtn_CY
          call mus_nNwtn_CY_load( me%CY, conf, nonNwt_table )
        case ( 'casson', 'Casson' )
          me%model = mus_nNwtn_CS
          call mus_nNwtn_CS_load( me%CS, conf, nonNwt_table )
        case default
          ! try to read in Power Law (PL) model parameters
          call mus_nNwtn_PL_load( me%PL, conf, nonNwt_table )
      end select

    end if

    call aot_table_close( L=conf, thandle = nonNwt_table )

  end subroutine mus_nNwtn_load
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> Read in the nonNewtonian Power Law (PL) model parameters from Lua file
  subroutine mus_nNwtn_PL_load( me, conf, tHandle )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian type
    type( mus_nNwtn_PL_type ), intent(out) :: me
    !> lua state
    type( flu_state ), intent(inout) :: conf
    !> nonNewtonian table handle
    integer, intent(inout) :: tHandle
    ! ---------------------------------------------------------------------------
    integer :: iError
    ! ---------------------------------------------------------------------------

    ! load n
    call aot_get_val( L       = conf,         &
      &               thandle = tHandle,      &
      &               key     = 'n',          &
      &               val     = me%n,         &
      &               default = 0.5_rk,       &
      &               ErrCode = iError        )

    ! load k
    call aot_get_val( L       = conf,         &
      &               thandle = tHandle,      &
      &               key     = 'k',          &
      &               val     = me%k,         &
      &               default = 0.0000035_rk, &
      &               ErrCode = iError        )

    me%n_1 = me%n - 1._rk

  end subroutine mus_nNwtn_PL_load
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> Read in the nonNewtonian Carreau-Yasuda (CY) model parameters from Lua file
  subroutine mus_nNwtn_CY_load( me, conf, tHandle )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian type
    type( mus_nNwtn_CY_type ), intent(out) :: me
    !> lua state
    type( flu_state ), intent(inout) :: conf
    !> nonNewtonian table handle
    integer, intent(inout) :: tHandle
    ! ---------------------------------------------------------------------------
    integer :: iError
    ! ---------------------------------------------------------------------------

    ! load visc0
    call aot_get_val(L       = conf,         &
      &              thandle = tHandle,      &
      &              key     = 'visc0',      &
      &              val     = me%visc0,     &
      &              default = 0.16_rk,      &
      &              ErrCode = iError        )

    ! load viscInf
    call aot_get_val( L       = conf,         &
      &               thandle = tHandle,      &
      &               key     = 'viscInf',    &
      &               val     = me%viscInf,   &
      &               default = 0.0035_rk,    &
      &               ErrCode = iError        )

    ! load lambda
    call aot_get_val( L       = conf,         &
      &               thandle = tHandle,      &
      &               key     = 'lambda',     &
      &               val     = me%lambda,    &
      &               default = 8.2_rk,       &
      &               ErrCode = iError        )

    ! load a
    call aot_get_val( L       = conf,         &
      &               thandle = tHandle,      &
      &               key     = 'a',          &
      &               val     = me%a,         &
      &               default = 0.64_rk,      &
      &               ErrCode = iError        )

    ! load n
    call aot_get_val( L       = conf,         &
      &               thandle = tHandle,      &
      &               key     = 'n',          &
      &               val     = me%n,         &
      &               default = 0.2128_rk,    &
      &               ErrCode = iError        )

    ! calculate intermediate parameter
    me%n_1_a = ( me%n - 1._rk ) / me%a

  end subroutine mus_nNwtn_CY_load
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> Read in the nonNewtonian Casson model parameters from Lua file
  subroutine mus_nNwtn_CS_load( me, conf, tHandle )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian type
    type( mus_nNwtn_CS_type ), intent(out) :: me
    !> lua state
    type( flu_state ), intent(inout) :: conf
    !> nonNewtonian table handle
    integer, intent(inout) :: tHandle
    ! ---------------------------------------------------------------------------
    integer :: iError
    ! ---------------------------------------------------------------------------

    ! load k0
    call aot_get_val(L       = conf,      &
      &              tHandle = tHandle,   &
      &              key     = 'k0',      &
      &              val     = me%k0,     &
      &              default = 0.1937_rk, &
      &              ErrCode = iError     )

    ! load k1
    call aot_get_val( L       = conf,     &
      &               tHandle = tHandle,  &
      &               key     = 'k1',     &
      &               val     = me%k1,    &
      &               default = 0.055_rk, &
      &               ErrCode = iError    )

  end subroutine mus_nNwtn_CS_load
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> write nonNewtonian fluid parameters into a lua file
  !!
  subroutine mus_nNwtn_save2lua( me, conf )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian parameters
    type( mus_nNwtn_type ), intent(in) :: me
    type( aot_out_type ) :: conf
    ! ---------------------------------------------------------------------------

    call aot_out_open_table( put_conf = conf, tname = 'nonNewtonian' )

    call aot_out_val( put_conf = conf,          &
      &               vname    = 'model',       &
      &               val      = trim(me%label) )

    select case( me%model )
      case ( mus_nNwtn_PL )
        call mus_nNwtn_PL_save( me%PL, conf )
      case ( mus_nNwtn_CY )
        call mus_nNwtn_CY_save( me%CY, conf )
      case ( mus_nNwtn_CS )
        call mus_nNwtn_CS_save( me%CS, conf )
    end select

    call aot_out_close_table( put_conf = conf )

  end subroutine mus_nNwtn_save2lua
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> write nonNewtonian Power Law (PL) parameters into a lua file
  !!
  subroutine mus_nNwtn_PL_save( me, conf )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian parameters
    type( mus_nNwtn_PL_type ), intent(in) :: me
    type( aot_out_type ) :: conf
    ! ---------------------------------------------------------------------------

    call aot_out_val( put_conf = conf,                                         &
      &               vname    = 'n',                                          &
      &               val      = me%n )
    call aot_out_val( put_conf = conf,                                         &
      &               vname    = 'k',                                          &
      &               val      = me%k )

  end subroutine mus_nNwtn_PL_save
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> write nonNewtonian (CY) parameters into a lua file
  !!
  subroutine mus_nNwtn_CY_save( me, conf )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian parameters
    type( mus_nNwtn_CY_type ), intent(in) :: me
    type( aot_out_type ) :: conf
    ! ---------------------------------------------------------------------------

    call aot_out_val( put_conf = conf,  &
      &               vname    = 'n',   &
      &               val      = me%n   )
    call aot_out_val( put_conf = conf,  &
      &               vname    = 'a',   &
      &               val      = me%a   )
    call aot_out_val( put_conf = conf,      &
      &               vname    = 'lambda',  &
      &               val      = me%lambda  )
    call aot_out_val( put_conf = conf,     &
      &               vname    = 'visc0',  &
      &               val      = me%visc0  )
    call aot_out_val( put_conf = conf,        &
      &               vname    = 'viscInf',   &
      &               val      = me%viscInf   )

  end subroutine mus_nNwtn_CY_save
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> write nonNewtonian Casson parameters into a lua file
  !!
  subroutine mus_nNwtn_CS_save( me, conf )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian parameters
    type( mus_nNwtn_CS_type ), intent(in) :: me
    type( aot_out_type ) :: conf
    ! ---------------------------------------------------------------------------

    call aot_out_val( put_conf = conf, &
      &               vname    = 'k0', &
      &               val      = me%k0 )
    call aot_out_val( put_conf = conf, &
      &               vname    = 'k1', &
      &               val      = me%k1 )

  end subroutine mus_nNwtn_CS_save
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> Dump nonNewtonian fluid parameters to outUnit
  !!
  subroutine mus_nNwtn_dump2outUnit( me, outUnit )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian parameters
    type( mus_nNwtn_type ), intent(in) :: me
    integer, intent(in) :: outUnit
    ! ---------------------------------------------------------------------------

    if ( me%active ) then
      write(outUnit,'(A)') 'nonNewtonian fluid parameters:'
      write(outUnit,'(A)') '  model label: ', trim(me%label)

      ! dump model parameters by calling model dumper
      select case( me%model )
        case ( mus_nNwtn_PL )
          call mus_nNwtn_PL_dump( me%PL, outUnit )
        case ( mus_nNwtn_CY )
          call mus_nNwtn_CY_dump( me%CY, outUnit )
        case ( mus_nNwtn_CS )
          call mus_nNwtn_CS_dump( me%CS, outUnit )
      end select

    else
      write(outUnit,'(A)') 'No nonNewtonian table defined.'
    end if

    call tem_horizontalSpacer( fUnit = outUnit )

  end subroutine mus_nNwtn_dump2outUnit
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> Dump nonNewtonian Power Law (PL) parameters to outUnit
  !!
  subroutine mus_nNwtn_PL_dump( me, outUnit )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian parameters
    type( mus_nNwtn_PL_type ), intent(in) :: me
    integer, intent(in) :: outUnit
    ! ---------------------------------------------------------------------------

    write(outUnit,"( '  n = ', F8.4)") me%n
    write(outUnit,"( '  k = ', F8.4)") me%k

  end subroutine mus_nNwtn_PL_dump
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> Dump nonNewtonian (CY) parameters to outUnit
  !!
  subroutine mus_nNwtn_CY_dump( me, outUnit )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian parameters
    type( mus_nNwtn_CY_type ), intent(in) :: me
    integer, intent(in) :: outUnit
    ! ---------------------------------------------------------------------------

    write(outUnit, "('  n       = ', F8.4)") me%n
    write(outUnit, "('  a       = ', F8.4)") me%a
    write(outUnit, "('  lambda  = ', F8.4)") me%lambda
    write(outUnit, "('  visc0   = ', F8.4)") me%visc0
    write(outUnit, "('  viscInf = ', F8.4)") me%viscInf

  end subroutine mus_nNwtn_CY_dump
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> Dump nonNewtonian (CY) parameters to outUnit
  !!
  subroutine mus_nNwtn_CS_dump( me, outUnit )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian parameters
    type( mus_nNwtn_CS_type ), intent(in) :: me
    integer, intent(in) :: outUnit
    ! ---------------------------------------------------------------------------

    write(outUnit,"('  k0 = ', F8.4)") me%k0
    write(outUnit,"('  k1 = ', F8.4)") me%k1

  end subroutine mus_nNwtn_CS_dump
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> nonNewtonian power-law model
  !!
  real(kind=rk) function viscPhy_PL( me, shearRate )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian parameters
    type( mus_nNwtn_type ), intent(in)  :: me
    real(kind=rk),          intent(in)  :: shearRate
    ! ---------------------------------------------------------------------------

    viscPhy_PL = ( shearRate ** me%PL%n_1 ) * me%PL%k

  end function viscPhy_PL
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> nonNewtonian Casson model
  !!
  real(kind=rk) function viscPhy_CS( me, shearRate )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian parameters
    type( mus_nNwtn_type ), intent(in)  :: me
    real(kind=rk),           intent(in)  :: shearRate
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: t

    t = ( me%CS%k0 + me%CS%k1 * sqrt(shearRate) )
    viscPhy_CS = t * t / shearRate

  end function viscPhy_CS
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> nonNewtonian Carreau-Yasuda model
  !!
  real(kind=rk) function viscPhy_CY( me, shearRate )
    ! ---------------------------------------------------------------------------
    !> nonNewtonian parameters
    type( mus_nNwtn_type ), intent(in)  :: me
    real(kind=rk),           intent(in)  :: shearRate
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: t

    t = ( 1._rk + (me%CY%lambda*shearRate) ** me%CY%a ) ** me%CY%n_1_a
    viscPhy_CY = me%CY%viscInf + ( me%CY%visc0 - me%CY%viscInf ) * t

  end function viscPhy_CY
  ! **************************************************************************** !

end module mus_nonNewtonian_module
! ****************************************************************************** !
