main Program

program main
  !>--------------------------------------------------------------------------
  ! Program: Flow Past Cylinder
  ! Author: Mohamed Houssem Kasbaoui
  !
  ! Solver: CDIFS
  !
  ! Description: Flow past a cylinder, represented as an IB.
  !
  ! References:
  ! Dave, H., Herrmann, M. & Kasbaoui, M. H. The volume-filtering immersed
  ! boundary method. Journal of Computational Physics 487, 112136 (2023).
  ! --------------------------------------------------------------------------
  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 the immersed boundary
  call SetUpCaseIB()

  ! 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.0_wp, -0.5_wp*L(2), -0.5_wp*L(3)]
      xhi=[L(1),    0.5_wp*L(2),  0.5_wp*L(3)]
      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 field.
      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 SetUpCaseIB
      !> Builds and writes immersed boundary.
      use immersed_boundaries
      implicit none
      type(marker_set)    :: IB
      character(str64)    :: filename
      real(wp)            :: diam
      real(wp)            :: xc(3)
      real(wp)            :: d_axis
      integer             :: n_axis
      real(wp)            :: d_theta
      integer             :: n_theta
      integer             :: i,j,m
      real(wp), parameter :: Pi=4.0*atan(1.0_wp)

      ! Get info from parser
      call parser%Get("IB IC file",        filename)
      call parser%Get("Cylinder diameter", diam    )
      call parser%Get("Cylinder center",   xc      )

      ! Initialize Immersed Boundaries
      call IB%Initialize('IB',block,parallel)

      ! Create a cylinder
      if (parallel%RankIsRoot()) then

        ! Count the number of markers
        d_axis = 0.5_wp*block%dx(3)
        n_axis = max(1,floor((block%pmax(3)-block%pmin(3))/d_axis))
        d_axis = (block%pmax(3)-block%pmin(3))/real(n_axis,wp)

        n_theta=ceiling(Pi*diam/d_axis)
        d_theta=2.0_wp*Pi/real(n_theta,wp)

        ! Create markers
        call IB%Resize(n_theta*n_axis)

        select type (markers=>IB%p)
        type is (marker_obj)
          do i=1,n_theta
            do j=1,n_axis
              ! ID
              m = n_axis*(i-1)+j
              markers(m)%id  = int(m,kind=8)
              ! Position
              markers(m)%p(1)= xc(1)+0.5_wp*diam*cos((i-1)*d_theta + d_theta/2.0_wp)
              markers(m)%p(2)= xc(2)+0.5_wp*diam*sin((i-1)*d_theta + d_theta/2.0_wp)
              markers(m)%p(3)= xc(3)+(j-1)*d_axis + d_axis/2.0_wp + block%pmin(3)
              ! Surface Area
              markers(m)%SA  = 0.5_wp*diam*d_theta*d_axis
              ! Normal
              markers(m)%n(1)= cos((i-1)*d_theta + d_theta/2.0_wp)
              markers(m)%n(2)= sin((i-1)*d_theta + d_theta/2.0_wp)
              markers(m)%n(3)= 0.0_wp
              ! velocity
              markers(m)%v   = 0.0_wp
            end do
          end do
        end select
      end if

      ! Treatment for periodicity
      call IB%ApplyPeriodicity()

      ! Send to the right rank
      call IB%Communicate()

      ! Localize marker points
      call IB%Localize()

      ! Write data to disk
      call IB%SetWriteFileName(filename)
      call IB%Write(0,0.0_WP)

      ! Finalize
      call IB%Finalize()
      return
    end subroutine SetUpCaseIB
    subroutine SetUpCaseBCS()
      !> Defines boundary conditions.
      use leapBC
      implicit none
      ! Work variables
      type(bc_set)      :: bcs
      real(wp)          :: xhi(3),xlo(3)
      real(wp)          :: Vin
      real(wp), pointer :: BCVal(:,:,:)

      ! Get info from parser
      call parser%Get("Inflow velocity", Vin,    0.0_wp )

      ! 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 = [pmin(1),pmax(2),pmax(3)]
        call bcs%Add('xL', xlo, xhi, normal = '+x1')

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

        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

      ! Flow boundary conditions
      call bcs%SetBC('yL', BC_OUTFLOW)
      call bcs%SetBC('yR', BC_OUTFLOW)
      call bcs%SetBC('xR', BC_OUTFLOW)
      call bcs%SetBC('xL', BC_INFLOW )

      call bcs%GetBCPointer('xL','V1',BCVal); BCVal=Vin
      call bcs%GetBCPointer('xL','V2',BCVal); BCVal=0.0_wp
      call bcs%GetBCPointer('xL','V3',BCVal); BCVal=0.0_wp


      ! Solid volume fraction boundary conditions
      call bcs%SetBC('xL', BC_DIRICHLET,'ibVF')
      call bcs%SetBC('xR', BC_DIRICHLET,'ibVF')
      call bcs%SetBC('yL', BC_DIRICHLET,'ibVF')
      call bcs%SetBC('yR', BC_DIRICHLET,'ibVF')

      call bcs%GetBCPointer('xL','ibVF',BCVal); BCVal=0.0_wp
      call bcs%GetBCPointer('xR','ibVF',BCVal); BCVal=0.0_wp
      call bcs%GetBCPointer('yL','ibVF',BCVal); BCVal=0.0_wp
      call bcs%GetBCPointer('yR','ibVF',BCVal); BCVal=0.0_wp

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

      ! Clear data
      call bcs%Finalize()

      return
    end subroutine SetUpCaseBCS
end program main