main Program

program main
  !>--------------------------------------------------------------------------
  ! Program: Head-on Collision
  ! Author: Mohamed Houssem Kasbaoui
  !
  ! Solver: GRANS
  !
  ! Description: Head-on Collision. Performs head-on collisions between
  ! point particles (PP) and/or resolved particles (RP).
  !
  ! References:
  ! --------------------------------------------------------------------------
  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 point particles
  call SetUpCasePP()

  ! Set the resolved particles
  call SetUpCaseRP()

  ! 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([.true.,.true.,.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 SetUpCasePP()
      use particles_point
      !> Builds and writes point particles.
      implicit none
      ! Work variables
      type(particle_set) :: PP
      character(str64)   :: filename
      character(str64)   :: part_type
      real(wp)           :: L(3)
      real(wp)           :: dp
      real(wp)           :: rhop
      integer            :: n
      ! Get info from parser
      call parser%Get("PP IC file",  filename)
      call parser%Get("Particle diameter",   dp       )
      call parser%Get("Particle density",    rhop     )
      call parser%Get("Particle type",       part_type, default = 'default')

      ! Initialize point particle data data
      call PP%Initialize('PP',block,parallel)
      call PP%ChangePartType(part_type)


      if (parallel%RankIsRoot()) then

        L = block%pmax - block%pmin

        ! Activate Np particles
        call PP%Resize(2)

        ! Place particles in head-on collisions
        select type (particle => PP%p)
        class is (particle_obj)
          do n=1,PP%count_
            ! Particle global ID
            particle(n)%id = int(n, kind=8)
            ! Particle diameter
            particle(n)%d  = dp
            select case (n)
            case (1)
              ! Particle position
              particle(n)%p(1) = 0.5_wp*(block%pmin(1) + block%pmax(1))
              particle(n)%p(2) = 0.5_wp*(block%pmin(2) + block%pmax(2)) - particle(n)%d
              particle(n)%p(3) = 0.5_wp*(block%pmin(3) + block%pmax(3))
              ! Particle velocity
              particle(n)%v  = [0.0_wp, 1.0_wp, 0.0_wp]
            case (2)
              ! Particle position
              particle(n)%p(1) = 0.5_wp*(block%pmin(1) + block%pmax(1))
              particle(n)%p(2) = 0.5_wp*(block%pmin(2) + block%pmax(2)) + particle(n)%d
              particle(n)%p(3) = 0.5_wp*(block%pmin(3) + block%pmax(3))
              ! Particle velocity
              particle(n)%v  = [0.0_wp, -1.0_wp, 0.0_wp]
            end select
            ! Particle density
            particle(n)%rho= rhop
            ! Force and Torque on particle
            particle(n)%Fh = 0.0_wp
            particle(n)%Th = 0.0_wp
            particle(n)%Fc = 0.0_wp
            particle(n)%Tc = 0.0_wp
          end do
        end select

      end if

      ! Apply periodicity
      call PP%ApplyPeriodicity()

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

      ! Clear data
      call PP%Finalize()
      return
    end subroutine SetUpCasePP
    subroutine SetUpCaseRP()
      !> Builds and writes resolved particles.
      use particles_resolved
      implicit none
      ! Work variables
      type(ResPart_set) :: RP
      character(str64)  :: filename
      real(wp)          :: diam
      real(wp)          :: rhop
      real(wp)          :: dl
      integer           :: n

      call parser%Get('RP IC file',            filename  )
      call parser%Get("Particle diameter",     diam      )
      call parser%Get("Particle density",      rhop      )

      ! Initialze resolved particles
      call RP%Initialize('ResPart',block,parallel)

      if (parallel%RankIsRoot()) then

        ! Activate only 1 particle on this Rank
        call RP%Resize(2)

        ! IB resolution
        dl=0.5_wp*minval(block%dx)

        select type(particle => RP%p)
        type is (ResPart_obj)
          do n=1,RP%count_
            ! Particle globabl ID
            particle(n)%id   = int(n,kind=8)
            ! Diameter
            particle(n)%d    = diam
            select case (n)
            case (1)
              ! Particle position
              particle(n)%p(1) = 0.5_wp*(block%pmin(1) + block%pmax(1)) - particle(n)%d
              particle(n)%p(2) = 0.5_wp*(block%pmin(2) + block%pmax(2))
              particle(n)%p(3) = 0.5_wp*(block%pmin(3) + block%pmax(3))
              ! Particle velocity
              particle(n)%v  = [ 1.0_wp, 0.0_wp, 0.0_wp]
            case (2)
              ! Particle position
              particle(n)%p(1) = 0.5_wp*(block%pmin(1) + block%pmax(1)) + particle(n)%d
              particle(n)%p(2) = 0.5_wp*(block%pmin(2) + block%pmax(2))
              particle(n)%p(3) = 0.5_wp*(block%pmin(3) + block%pmax(3))
              ! Particle velocity
              particle(n)%v  = [-1.0_wp, 0.0_wp, 0.0_wp]
            end select
            ! Angular velocity
            particle(n)%w    = 0.0_wp
            ! Density
            particle(n)%rho  = rhop
            ! Zero force and torque
            particle(n)%Fh   = 0.0_wp
            particle(n)%Th   = 0.0_wp
            particle(n)%Fc   = 0.0_wp
            particle(n)%Tc   = 0.0_wp

            ! Add surface markers
            call RP%ib%AddSphere(particle(n)%p,particle(n)%d/2.0_wp,particle(n)%v,dl,particle(n)%id)
          end do
        end select

      end if

      ! Treatment for periodicity
      call RP%ApplyPeriodicity
      call RP%ib%ApplyPeriodicity

      ! Send to the right rank
      call RP%Communicate
      call RP%ib%Communicate

      ! Localize centers and markers on grid
      call RP%Localize
      call RP%ib%Localize

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

      ! Finalize
      call RP%Finalize
    end subroutine SetUpCaseRP
    subroutine SetUpCaseBCS()
      !> Defines boundary conditions.
      use leapBC
      implicit none
      ! Work variables
      type(bc_set) :: bcs

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

      ! Fully-periodic, nothing to do

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

      ! Clear data
      call bcs%Finalize()

      return
    end subroutine SetUpCaseBCS
end program main