main Program

program main
  !>--------------------------------------------------------------------------
  ! Program: Homogeneous Cooling
  ! Author: Mohamed Houssem Kasbaoui
  !
  ! Solver: GRANS
  !
  ! Description: Simulates the cooling of an initially homogeneous system
  ! (i.e., randomly placed) of particles (either point particles (PP) and/or
  ! resolved particles (RP)) due to inelastic collisions.
  !
  ! References:
  ! Haff, P. K., 1983, “Grain Flow as a Fluid-Mechanical Phenomenon,”
  ! Journal of Fluid Mechanics, 134, pp. 401–430.
  ! --------------------------------------------------------------------------
  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 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)
      integer            :: Np
      real(wp)           :: dp
      real(wp)           :: VFP
      real(wp)           :: Vmag
      real(wp)           :: Vp
      real(wp)           :: Lp
      real(wp)           :: rhop
      real(wp)           :: rand(3)
      integer            :: i,j,k,n
      integer            :: iimax,jjmax,kkmax
      real(wp),parameter :: Pi=4.0_wp*atan(1.0_wp)

      ! Get info from parser
      call parser%Get("PP IC file",  filename)
      call parser%Get("Particle vol. frac.", VFP      )
      call parser%Get("Particle diameter",   dp       )
      call parser%Get("Particle density",    rhop     )
      call parser%Get("Particle velocity",   Vmag     )
      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
        ! Particle volume
        Vp    = Pi/6.0_wp*dp**3
        ! Particle spacing
        LP    = (VP/VFP)**(1.0_wp/3.0_wp)
        ! Particle lattice
        iimax = floor( L(1)/Lp)
        jjmax = floor( L(2)/Lp)
        kkmax = floor( L(3)/Lp)

        Np    = iimax*jjmax*kkmax

        ! Activate Np particles
        call PP%Resize(Np)

        ! Place particles in head-on collisions
        select type (particle => PP%p)
        class is (particle_obj)
          do k=1,kkmax
            do j=1,jjmax
              do i=1,iimax
                n = i + (j-1)*iimax + (k-1)*jjmax*iimax
                ! Particle global ID
                particle(n)%id = int(n, kind=8)
                ! Particle diameter
                particle(n)%d  = dp
                ! Particle position
                particle(n)%p(1) = block%pmin(1) + real(i-0.5_wp,wp)*L(1)/real(iimax,wp)
                particle(n)%p(2) = block%pmin(2) + real(j-0.5_wp,wp)*L(2)/real(jjmax,wp)
                particle(n)%p(3) = block%pmin(3) + real(k-0.5_wp,wp)*L(3)/real(kkmax,wp)
                ! Particle velocity
                call random_number(rand)
                particle(n)%v  = Vmag*(rand-0.5_wp)
                ! 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 do
          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 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