main Program

program main
  !>--------------------------------------------------------------------------
  ! Program: Channel flow
  ! Author: Mohamed Houssem Kasbaoui
  !
  ! Solver: CDIFS
  !
  ! Description: Laminar/Turbulent Channel Flow
  !
  ! References:
  ! 1. H. Dave, M.H. Kasbaoui, Mechanisms of drag reduction by semidilute inertial
  ! particles in turbulent channel flow, Phys. Rev. Fluids 8 (2023) 084305.
  ! 2. J. Kim, P. Moin, and R. Moser, Turbulence Statistics in Fully Developed
  ! Channel Flow at Low Reynolds Number, Journal of Fluid Mechanics 177, 133 (1987).
  ! 3. R. B. Dean, Reynolds Number Dependence of Skin Friction and Other Bulk Flow
  ! Variables in Two-Dimensional Rectangular Duct Flow, Journal of Fluids
  ! Engineering 100, 215 (1978).
  ! --------------------------------------------------------------------------
  use leapKinds
  use leapParser
  use leapParallel
  use leapIO
  use leapBlock
  use leapEulerian
  implicit none
  type(parallel_obj) :: parallel                                               !! Utility that handles parallel (MPI) functions
  type(parser_obj)   :: parser                                                 !! Utility that parses input files
  type(block_obj)    :: block                                                  !! Block object manages the Cartesian grid

  ! Initialize parser
  call parser%Initialize()

  ! Parse input file
  call parser%ParseFile()

  ! Initialize parallel environment
  call parallel%Initialize()

  ! Set the block info
  call SetUpCaseBlock()

  ! Set the initial fields
  call SetUpCaseFields()

  ! Set boundary conditions
  call SetUpCaseBCS()

  ! Free up data
  call block%Finalize()
  call parser%Finalize()
  call parallel%Finalize()
  contains
    subroutine SetUpCaseBlock()
      !> Builds and writes block file.
      implicit none
      ! Work variables
      character(str64) :: filename
      real(wp)         :: L(3)
      integer          :: N(3)
      integer          :: ngc
      integer          :: Nb(3)
      real(wp)         :: r
      integer          :: i,j,k,dir

      ! Get info from parser
      call parser%Get("Block file",        filename)
      call parser%Get("Domain size",       L       )
      call parser%Get("Grid points",       N       )
      call parser%Get("Ghost cells",       ngc     )
      call parser%Get("Stretching factor", r       )
      call parser%Get("Partition",         Nb      )

      ! Initialize the main block
      call block%Initialize(ngc,parallel)

      ! Setup the domain periodicity
      call block%SetPeriodicity([.true.,.false.,.true.])

      ! Create a grid, locally refined at the wall
      ! ------------------------------------------
      associate (lo =>block%lo, hi=>block%hi)
        ! Setup bounds
        lo = [1,1,1]
        hi = N

        ! Initialize axes
        do dir=1,3
          call block%axis(dir)%Initialize(lo(dir),hi(dir),ngc)
        end do
        call block%SetConveniencePointers()
      end associate

      associate (x=>block%x, y=>block%y, z=>block%z, &
        lo =>block%lo, hi=>block%hi)
        ! Uniform x-axis
        do i=lo(1),hi(1)+1
          x(i)= (i-lo(1))*L(1)/real(hi(1)-lo(1)+1,wp)
        end do

        ! Stretched y-axis
        do j=lo(2),hi(2)+1
          y(j)= (j-lo(2))*L(2)/real(hi(2)-lo(2)+1,wp) - 0.5_wp*L(2)
          y(j)= 0.5_wp*L(2)*tanh(r*y(j)/(0.5_wp*L(2)))/tanh(r)
        end do

        ! Uniform z-axis
        do k=lo(3),hi(3)+1
          z(k)= (k-lo(3))*L(3)/real(hi(3)-lo(3)+1,wp) - 0.5_wp*L(3)
        end do
      end associate

      ! Partition block for parallel initializations
      call block%Partition(Nb)

      ! Write block to disk
      call block%Write(filename)

      return
    end subroutine SetUpCaseBlock
    subroutine SetUpCaseFields()
      !> Builds and writes initial flow fields.
      implicit none
      ! Work variables
      type(Eulerian_set)   :: fields
      type(eulerian_obj_r) :: V(3)
      type(eulerian_obj_r) :: P
      character(str64)     :: filename
      real(wp)             :: L(3)
      real(wp)             :: Reb
      real(wp)             :: h
      real(wp)             :: Vbulk
      real(wp)             :: Vfric
      real(wp)             :: Vlaminar
      real(wp)             :: rand
      real(wp)             :: amp
      real(wp)             :: rho
      real(wp)             :: mu
      real(wp)             :: Cf
      real(wp)             :: tau_w
      real(wp)             :: pgrad
      real(wp)             :: ymin,ymax
      integer              :: i,j,k
      real(wp), parameter  :: twoPi=8.0_wp*atan(1.0_wp)
      real(wp), parameter  :: Re_transition = 650.0_wp

      ! Get info from parser
      call parser%Get("Fields IC file",        filename)
      call parser%Get("Domain size",           L       )
      call parser%Get('Bulk Reynolds',         Reb     )
      call parser%Get('Fluctuation amplitude', amp     )
      call parser%Get("Fluid density",         rho     )
      call parser%Get("Fluid viscosity",       mu      )

      ! Initialize fields container
      call fields%Initialize(block,parallel)

      ! Add fields to container (this will allocate data)
      call fields%Add('V1', 1, V(1))
      call fields%Add('V2', 2, V(2))
      call fields%Add('V3', 3, V(3))
      call fields%Add('P',  0, P   )


      ! Channel half height and mean velocity
      h     = 0.5_wp*L(2)
      Vbulk = Reb*(mu/rho)/h

      associate (lo => block%lo, hi=> block%hi,          &
          x =>block%x , y =>block%y,  z => block%z, &
          xm=>block%xm, ym=>block%ym, zm=> block%zm)

        V(1) = 0.0_wp
        V(2) = 0.0_wp
        V(3) = 0.0_wp
        P    = 0.0_wp

        ! Base laminar flow
        ymin = block%pmin(2)
        ymax = block%pmax(2)
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)
              ! Base laminar contribution
              Vlaminar = Vbulk*6.0_wp*(ym(j)-ymin)*(ymax-ym(j))/(ymax-ymin)**2
              ! Add perturbation to make transition faster
              call random_number(rand)
              V(1)%cell(i,j,k) = Vlaminar + Vbulk*amp*(rand-0.5_wp)*cos(8.0_wp*twoPi*zm(k)/L(3))
              call random_number(rand)
              V(3)%cell(i,j,k) = 0.0_wp  +  Vbulk*amp*(rand-0.5_wp)*cos(8.0_wp*twoPi*xm(i)/L(1))
            end do
          end do
        end do
      end associate

      ! Determine friction factor
      if (Reb.le.Re_transition) then
        Cf    = 6.0_wp/Reb
      else
        ! Dean's formula
        Cf    = 0.0614_wp/Reb**0.25_wp
      end if

      ! Determine other flow parameters
      tau_w = 0.5_wp*rho*Vbulk**2*Cf
      pgrad = -tau_w/h
      Vfric = sqrt(tau_w/rho)

      ! Write some info to stdout
      if (parallel%RankIsRoot()) then
        write(*,*) "Bulk Reynolds     = ", Reb
        write(*,*) "Friction Reynolds = ", rho*Vfric*h/mu
        write(*,*) "Bulk velocity     = ", Vbulk
        write(*,*) "Friction velocity = ", Vfric
        write(*,*) "Pressure gradient = ", pgrad
        write(*,*) "min dx(1)/del_v   = ", block%dx(1)*(rho*Vfric)/mu
        write(*,*) "min dx(2)/del_v   = ", block%dx(2)*(rho*Vfric)/mu
        write(*,*) "min dx(3)/del_v   = ", block%dx(3)*(rho*Vfric)/mu
      end if

      ! Write data to disk
      call fields%SetWriteFileName(filename)
      call fields%Write(0,0.0_wp)

      ! Clear data
      call fields%Finalize()

      return
    end subroutine SetUpCaseFields
    subroutine SetUpCaseBCS()
      !> Defines boundary conditions.
      use leapBC
      implicit none
      ! Work variables
      type(bc_set) :: bcs
      real(wp)     :: xhi(3),xlo(3)

      ! Initialize utility that handles boundary conditions
      call bcs%Initialize(block,parallel)

      ! Setup regions where boundary conditions apply
      associate (pmin => block%pmin, pmax => block%pmax)

        xlo = [pmin(1),pmin(2),pmin(3)]
        xhi = [pmax(1),pmin(2),pmax(3)]
        call bcs%Add('yL', xlo, xhi, normal = '+x2')

        xlo = [pmin(1),pmax(2),pmin(3)]
        xhi = [pmax(1),pmax(2),pmax(3)]
        call bcs%Add('yR', xlo, xhi, normal = '-x2')

      end associate

      call bcs%SetBC('yR', BC_WALL )
      call bcs%SetBC('yL', BC_WALL )

      ! Write boundary conditions
      call bcs%Write(0,0.0_wp)

      ! Clear data
      call bcs%Finalize()

      return
    end subroutine SetUpCaseBCS
end program main