main Program

program main
  !>--------------------------------------------------------------------------
  ! Program: Lid-driven cavity problem
  ! Author: Mohamed Houssem Kasbaoui
  !
  ! Solver: CDIFS
  !
  ! Description: Flow in cavity driven by a sliding lid
  !
  ! References:
  ! Ghia, Ghia, and Shin (1982), "High-Re solutions for incompressible flow
  ! using the Navier-Stokes equations and a multigrid method", Journal of
  ! Computational Physics, Vol. 48, pp. 387-411.
  ! --------------------------------------------------------------------------
  use leapKinds
  use leapParser
  use leapParallel
  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)         :: xlo(3)
      real(wp)         :: xhi(3)
      integer          :: ilo(3)
      integer          :: ihi(3)

      ! 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("Partition",   Nb      )

      ! Domain extents
      xlo=-0.5_wp*L ; xhi= 0.5_wp*L
      ilo=[1,1,1] ; ihi=N

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

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

      ! Create a uniform block
      call block%SetupUniformGrid(xlo,xhi,ilo,ihi)

      ! 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

      ! Get info from parser
      call parser%Get("Fields IC file",  filename)

      ! 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   )

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

      ! 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
      type(extent_obj)  :: extents
      integer           :: i,j,k
      real(wp)          :: xhi(3)
      real(wp)          :: xlo(3)
      real(wp)          :: L(3)
      real(wp)          :: Vlid
      real(wp), pointer :: V1bc(:,:,:)
      real(wp), pointer :: V2bc(:,:,:)
      real(wp), pointer :: V3bc(:,:,:)

      ! Get lid velocity and domain size
      call parser%Get("Lid velocity", Vlid, default = 1.0_wp  )
      call parser%Get("Domain size", L )

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

      ! Setup regions where boundary conditions apply
      ! - Left boundary
      xlo = [-0.5_wp*L(1),-0.5_wp*L(2),-0.5_wp*L(3)]
      xhi = [-0.5_wp*L(1), 0.5_wp*L(2), 0.5_wp*L(3)]
      call bcs%Add('Left', xlo, xhi, normal = '+x1')

      ! - Right boundary
      xlo = [ 0.5_wp*L(1),-0.5_wp*L(2),-0.5_wp*L(3)]
      xhi = [ 0.5_wp*L(1), 0.5_wp*L(2), 0.5_wp*L(3)]
      call bcs%Add('Right', xlo, xhi, normal = '-x1')

      ! - Bottom boundary
      xlo = [-0.5_wp*L(1),-0.5_wp*L(2),-0.5_wp*L(3)]
      xhi = [ 0.5_wp*L(1),-0.5_wp*L(2), 0.5_wp*L(3)]
      call bcs%Add('Bottom', xlo, xhi, normal = '+x2')

      ! - Top boundary
      xlo = [-0.5_wp*L(1), 0.5_wp*L(2),-0.5_wp*L(3)]
      xhi = [ 0.5_wp*L(1), 0.5_wp*L(2), 0.5_wp*L(3)]
      call bcs%Add('Top', xlo, xhi, normal = '-x2')

      call bcs%SetBC('Top',    BC_INFLOW)
      call bcs%SetBC('Bottom', BC_WALL  )
      call bcs%SetBC('Left',   BC_WALL  )
      call bcs%SetBC('Right',  BC_WALL  )

      ! Set the inflow values
      extents = bcs%GetExtents('Top')
      call bcs%GetBCPointer('Top','V1',V1bc)
      call bcs%GetBCPointer('Top','V2',V2bc)
      call bcs%GetBCPointer('Top','V3',V3bc)
      do k=extents%lo(3),extents%hi(3)
        do j=extents%lo(2),extents%hi(2)
          do i=extents%lo(1),extents%hi(1)
            V1bc(i,j,k) = Vlid
            V2bc(i,j,k) = 0.0_wp
            V3bc(i,j,k) = 0.0_wp
          end do
        end do
      end do

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

      ! Clear data
      V1bc => null()
      V2bc => null()
      V3bc => null()
      call bcs%Finalize()

      return
    end subroutine SetUpCaseBCS
end program main