
!> definition of models which are simulated: Euler, Navier-Stokes, convection-diffusion
module modelNS_mod

 ! use main_data
  use model_mod
!  use f_mapping
!  use mesh_oper
!  use define_state
!  use blocks_integ
!  use model3DNS
!  use model2DNS
!  use modelTurb2e
!  use modelFE
!  use modelLaplace

  implicit none

  !> Euler or the Navier-Stokes equations
  type, EXTENDS( Model_t ), public :: NavierStokes_t
     !   real :: kappa, kappa1
     !   real :: Pr

   contains
     procedure :: init => initNavierStokes

  end type NavierStokes_t

  !> Eulerian model with mass diffusion
  type, EXTENDS( Model_t ), public :: Eulerian_t
     !   real :: kappa, kappa1
     !   real :: Pr

   contains
     procedure :: init => initEulerian

  end type Eulerian_t

  !> RANS ( Navier-Stokes equations) + k-\omega model of turbulence
  type, EXTENDS( Model_t ), public :: RANS_2e_t
     !   real :: kappa, kappa1
     !   real :: Pr


   contains
     procedure :: init => initRANS_2e

  end type RANS_2e_t



  !> pedestrian flow
  type, EXTENDS( Model_t ), public :: Pedestrian_t
     !real :: kappa, kappa1
     !real :: Pr
     !real :: p0

   contains
     procedure :: init => initPedestrian

  end type Pedestrian_t

   !> pedestrian flow
  type, EXTENDS( Model_t ), public :: ShallowWater_t
     !real :: kappa, kappa1
     !real :: Pr
     !real :: p0

   contains
     procedure :: init => initShallowWater

  end type ShallowWater_t

  public :: initNavierStokes
  public :: initEulerian
  public :: initRANS_2e
  public :: initPedestrian
  public :: initShallowWater


 contains

 !> initialization of the NS case
 !> isca, t2 not used
  subroutine initNavierStokes( this, Re, isca, t2)
    class (NavierStokes_t), intent(inout) :: this
    real, intent(in) ::Re
    integer, intent(in), optional :: isca
    real, intent(in), optional :: t2

    this%ndim = 4

    this%Re = Re
    this%kappa = 1.4
    this%kappa1 = this%kappa - 1.
    this%Pr = 0.72   ! Prandtl number (almost all gases)
    this%convective = .true.

    ! gravity accelaration
    this%g981 = 9.81

    ! heat capacity, constant volume
    this%cV = 1.    ! dimensionless quantity

    ! heat capacity, constant pressure
    this%cP =  this%cV * this%kappa 

    
    this%varying_time_term = .false.     ! .true. = a function in front of the time derivative term
    this%precomputed_arrays = .false.    ! .true. capacity and conductivity precomputed

    if ( Re == 0.) then
        print*,' # Compressible Euler equations'
        this%Re1 = 0.

        this%discretizationOfSlipBC_Euler = 3 !2 ! 1 !3 != standard version in Adgfem
        !this%discretizationOfSlipBC_Euler = 3 !  BY Filip

        if (this%discretizationOfSlipBC_Euler == 1) then
            print*,' # H_BC(w) = p(w) (0,n1,n2,0)'
        else if (this%discretizationOfSlipBC_Euler == 2) then
             print*, "H_BC(w) = p( u_\Gamma(w)) (0,n1,n2,0) , where u_Gamma = UpdateMirror_W"
        else if (this%discretizationOfSlipBC_Euler == 3) then
             print*, "H_BC(w) =  H_VS(w,\Mir(w),n), where \Mir = Mirror_W"
        else
            stop "Problem in setting the dicretization of the slip BC for Euler eq!"
        end if

     elseif  ( Re > 0.) then
        this%Re1 = 1./this%Re
        print*,' # Compressible Navier-Stokes equations, Re=',this%Re

        this%discretizationOfSlipBC_Euler = 3 !2 ! 1 !3 != standard version in Adgfem
        print*, "H_BC(w) =  H_VS(w,\Mir(w),n), where \Mir = Mirror_W"

     else
        print*,'# Reynolds number is negative',this%Re,' STABILIZATION !!!'
        !stop
     endif

  end subroutine initNavierStokes

  !>  Eulerian model with mass diffusion, initialization
  !> isca, t2 not used
  subroutine initEulerian( this, Re, isca, t2)
    class (Eulerian_t), intent(inout) :: this
    real, intent(in) ::Re
    integer, intent(in), optional :: isca
    real, intent(in), optional :: t2

    ! TO BE MODIFIED
    this%ndim = 4

    this%Re = Re
    this%kappa = 1.4
    this%kappa1 = this%kappa - 1.
    this%Pr = 0.72   ! Prandtl number (almost all gases)

    if(t2 <= 0.) stop "non-positive alpha for Eulerian meodel"
    this%alpha = t2  ! parameter alpha of the diffusivity

    ! gravity accelaration
    this%g981 = 9.81

    ! heat capacity, constant volume
    this%cV = 1.    ! dimensionless quantity

    ! heat capacity, constant pressure
    this%cP =  this%cV * this%kappa 


    this%convective = .true.

    this%varying_time_term = .false.     ! .true. = a function in front of the time derivative term
    this%precomputed_arrays = .false.    ! .true. capacity and conductivity precomputed

    if ( Re == 0.) then
        print*,' # Eulerian model with mass diffusuion'
        this%Re1 = 0.

        this%discretizationOfSlipBC_Euler = 3 !2 ! 1 !3 != standard version in Adgfem

        if (this%discretizationOfSlipBC_Euler == 1) then
            print*,' # H_BC(w) = p(w) (0,n1,n2,0)'
        else if (this%discretizationOfSlipBC_Euler == 2) then
             print*, "H_BC(w) = p( u_\Gamma(w)) (0,n1,n2,0) , where u_Gamma = UpdateMirror_W"
        else if (this%discretizationOfSlipBC_Euler == 3) then
             print*, "H_BC(w) =  H_VS(w,\Mir(w),n), where \Mir = Mirror_W"
        else
            stop "Problem in setting the dicretization of the slip BC for Euler eq!"
        end if

     elseif  ( Re > 0.) then
        this%Re1 = 1./this%Re
        print*,' # Eulerian model with mass diffusuion, Re=',this%Re

        this%discretizationOfSlipBC_Euler = 1 !2 ! 1 !3 != standard version in Adgfem
        print*,' # H_BC(w) = p(w) (0,n1,n2,0)'
        !print*, "H_BC(w) =  H_VS(w,\Mir(w),n), where \Mir = Mirror_W"

     else
        print*,'# Reynolds number is negative',this%Re,' STABILIZATION !!!'
        !stop
     endif

   end subroutine initEulerian

 !> initialization of the RANS case with k-\omega model of turbulence
 !> isca, t2 not used
  subroutine initRANS_2e( this, Re, isca, t2)
    class (RANS_2e_t), intent(inout) :: this
    real, intent(in) ::Re
    integer, intent(in), optional :: isca
    real, intent(in), optional :: t2

    this%ndim = 4

    this%Re = Re
    this%kappa = 1.4
    this%kappa1 = this%kappa - 1.
    this%Pr = 0.72   ! Prandtl number (almost all gases)
    this%Pr_t = 0.9  ! turbulence Prandtl number (almost all gases)
    this%sigmak = 1./2
    this%sigmaw = 1./2
    this%cp = 1.4
    this%betak = 0
    this%betaw = 0
    this%alphaw = 0

    this%convective = .true.

    this%varying_time_term = .false.     ! .true. = a function in front of the time derivative term
    this%precomputed_arrays = .false.    ! .true. capacity and conductivity precomputed


    this%discretizationOfSlipBC_Euler = 3 !2 ! 1 !3 != standard version in Adgfem

    if  ( Re > 0.) then
       this%Re1 = 1./this%Re
       print*,' # Reynold averagede Navier-Stokes equations, Re=',this%Re

     else
        print*,'# Reynolds number for RANS is non-positive',this%Re
        stop
     endif

   end subroutine initRANS_2e


 !> initialization of the pedestrian flow problem
  subroutine initPedestrian( this, Re, isca, t2)
    class (Pedestrian_t), intent(inout) :: this
    real, intent(in) ::Re
    integer, intent(in), optional :: isca
    real, intent(in), optional :: t2

    this%ndim = 3
    this%convective = .true.

    this%varying_time_term = .false.     ! .true. = a function in front of the time derivative term
    this%precomputed_arrays = .false.    ! .true. capacity and conductivity precomputed

    this%Re = 0.  ! Re is p0
    this%Re1 = 0.
    !!this%Pr = 0.005  ! used as the minimal allowed density
    !this%Pr = 1E-3  ! used as the minimal allowed density  WORKS for P_0 approximation
    this%Pr = 5E-3  ! used as the minimal allowed density
    !this%Pr = 1E-2  ! used as the minimal allowed density


    this%kappa = t2
    this%kappa1 = this%kappa - 1.

    ! chacateristic density
    this%rho_char = 1.
    !this%rho_char = 0.1

    !this%p0 = Re
    this%p0 = Re * this%rho_char**this%kappa

    !this%v_max = 2       ! given in 'eikonal.f90', module PedestrianSettings
    !this%rho_max = 9
    !this%alpha_ped = 7.5
    this%tau = 0.61

    this%discretizationOfSlipBC_Euler = 3 !2 ! 1 !3 != standard version in Adgfem


    write(*,'(a34, f8.2,a8,es12.4)')' # Pedestrian flow (SWE): kappa =', this%kappa,',   p0 =', this%p0
    !  elseif  ( Re > 0.) then
    !     this%Re1 = 1./this%Re
    !     print*,' # Compressible Navier-Stokes equations, Re=',this%Re

    !  else
    !     print*,'# Reynolds number is negative',this%Re,' STABILIZATION !!!'
    !     !stop
    !  endif

  end subroutine initPedestrian



 !> initialization of the shallow water problem
 !> isca, t2 not used
  subroutine initShallowWater( this, Re, isca, t2)
    class (ShallowWater_t), intent(inout) :: this
    real, intent(in) ::Re
    integer, intent(in), optional :: isca
    real, intent(in), optional :: t2

    this%ndim = 3
    this%convective = .true.

    this%varying_time_term = .false.     ! .true. = a function in front of the time derivative term
    this%precomputed_arrays = .false.    ! .true. capacity and conductivity precomputed

    this%Re = 0.  ! Re is p0
    this%Re1 = 0.
    !!this%Pr = 0.005  ! used as the minimal allowed density
    !this%Pr = 1E-3  ! used as the minimal allowed density  WORKS for P_0 approximation
    this%Pr = 5E-3  ! used as the minimal allowed density
    !this%Pr = 1E-2  ! used as the minimal allowed density

    ! topology of the bottom
    this%icase = isca

    ! gravity accelaration
    this%g981 = 9.81

    this%kappa = t2
    this%kappa1 = this%kappa - 1.

    ! chacateristic density
    this%rho_char = 1.
    !this%rho_char = 0.1

    !this%p0 = Re
    this%p0 = Re * this%rho_char**this%kappa

    !this%v_max = 2       ! given in 'eikonal.f90', module PedestrianSettings
    !this%rho_max = 9
    !this%alpha_ped = 7.5
    this%tau = 0.61

    write(*,'(a34, f8.2,a8,es12.4)')' # Pedestrian flow (SWE): kappa =', this%kappa,',   p0 =', this%p0
    !  elseif  ( Re > 0.) then
    !     this%Re1 = 1./this%Re
    !     print*,' # Compressible Navier-Stokes equations, Re=',this%Re

    !  else
    !     print*,'# Reynolds number is negative',this%Re,' STABILIZATION !!!'
    !     !stop
    !  endif

  end subroutine initShallowWater

end module modelNS_mod
