main Program

program main
  !>--------------------------------------------------------------------------
  ! Program: Granular Channel
  ! Author: Mohamed Houssem Kasbaoui
  !
  ! Solver: GRANS
  !
  ! Description: Simulates purely granular transport of fully resolved
  ! particles in an Immersed Boundary channel.
  !
  ! References:
  ! --------------------------------------------------------------------------
  use leapKinds
  use leapParser
  use leapParallel
  use leapBlock
  use leapEulerian
  implicit none
  type(parallel_obj) :: parallel                                               !! Utility that interfaces with MPI API
  type(parser_obj)   :: parser                                                 !! Parses input files
  type(block_obj)    :: block                                                  !! block object
  real(wp)           :: IBy(2)                                                 !! IB vertical positions

  ! Initialize parser
  call parser%Initialize()

  ! Parse input file
  call parser%ParseFile()

  ! Initialize parallel environment
  call parallel%Initialize()

  ! Set the block info
  call SetUpCaseBlock()

  ! Activate optional components
  call SetUpCaseRP()
  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), ngc, Nb(3)
      integer          :: dir
      integer          :: i,j,k
      integer          :: nwall
      integer          :: hif
      integer          :: lof
      integer          :: Nf

      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                    )
      call parser%Get("Number of wall cells", nwall,    default = 4 )

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

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

      ! Create a uniform grid, with extra cells
      ! in 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

        ! y-axis
        Nf  = N(2)-2*nwall  ! Total number of points within fluid
        lof = lo(2)+nwall   ! Minimum position within fluid
        hif = lof + Nf - 1  ! Maximum position within fluid

        ! Inner channel points
        do j=lof,hif+1
          y(j)= (j-lof)*L(2)/real(Nf,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
        ! Bottom wall points
        do j=lo(2),lof-1
          y(j) = y(lof)-(lof-j)*(y(lof+1)-y(lof))
        end do
        ! Top wall points
        do j=hif+2,hi(2)+1
          y(j) = y(hif+1)+(j-(hif+1))*(y(hif+1)-y(hif))
        end do

        ! Store IB positions
        IBy(1) = y(lof)
        IBy(2) = y(hif+1)

        ! 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 SetUpCaseRP()
      !> Builds and writes resolved particles.
      use particles_resolved
      implicit none
      ! Work variables
      character(str64)   :: filename
      type(ResPart_set)  :: RP                                                 !! Resolved particles
      real(wp)           :: L(3)
      real(wp)           :: diam
      real(wp)           :: rhop
      real(wp)           :: dl
      real(wp)           :: VFP
      real(wp)           :: Vp
      real(wp)           :: Lp
      real(wp)           :: ymin
      real(wp)           :: ymax
      real(wp)           :: Vel
      real(wp)           :: amp
      real(wp)           :: rand(3)
      integer            :: iimax,jjmax,kkmax
      integer            :: i,j,k
      integer            :: n
      real(wp),parameter :: Pi=4.0_wp*atan(1.0_wp)

      ! Get case info
      call parser%Get('RP IC file',            filename               )
      call parser%Get('Domain size',           L                      )
      call parser%Get("Particle diameter",     diam                   )
      call parser%Get("Particle density",      rhop                   )
      call parser%Get("Particle velocity",     Vel                    )
      call parser%Get('Particle fluctuations', amp  , default = 0.0_wp)
      call parser%Get("Avg VFP",               VFP                    )

      ymin  = -0.5_wp*L(2)
      ymax  =  0.5_wp*L(2)

      ! Compute average spacing for uniform distribution
      Vp = Pi/6.0_wp*diam**3
      Lp = (Vp/VFP)**(1.0_wp/3.0_wp)
      iimax = floor(L(1)/Lp)
      jjmax = floor((IBy(2)-IBy(1))/Lp)
      kkmax = floor(L(3)/Lp)

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

      ! Total count of particles to seed
      RP%count = iimax*jjmax*kkmax

      ! Activate particles on Root
      if (parallel%RankIsRoot()) then

        call RP%Resize(RP%count)

        dl=0.5_wp*minval(block%dx)

        select type(particle => RP%p)
        type is (ResPart_obj)
          do k=1,kkmax
            do j=1,jjmax
              do i=1,iimax
                n = i + (j-1)*iimax + (k-1)*jjmax*iimax
                ! Particle globabl ID
                particle(n)%id   = int(n,kind=8)
                ! Diameter
                particle(n)%d    = diam
                ! Position
                particle(n)%p(1) = block%pmin(1) + real(i-0.5_wp,wp)*L(1)/real(iimax,wp)
                particle(n)%p(2) = IBy(1)        + real(j-0.5_wp,wp)*(IBy(2)-IBy(1))/real(jjmax,wp)
                particle(n)%p(3) = block%pmin(3) + real(k-0.5_wp,wp)*L(3)/real(kkmax,wp)
                ! Velocity
                call random_number(rand)
                particle(n)%v(1) = Vel      + Vel*amp*(rand(1)-0.5_wp)
                particle(n)%v(2) = 0.0_wp   + Vel*amp*(rand(2)-0.5_wp)
                particle(n)%v(3) = 0.0_wp   + Vel*amp*(rand(3)-0.5_wp)
                ! 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 do
          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

      return
    end subroutine SetUpCaseRP
    subroutine SetUpCaseIB()
      !> Builds and writes immersed boundary.
      use immersed_boundaries
      implicit none
      ! Work variables
      character(str64)  :: filename
      type(marker_set)  :: IB                                                  !! IB
      real(wp)          :: center(3)                                           !! IB plane centroid
      character(len=3)  :: normal                                              !! IB plane normal direction
      real(wp)          :: width(3)                                            !! IB plane width
      real(wp)          :: vel(3)                                              !! IB plane velocity
      real(wp)          :: dl                                                  !! IB plane resolution (marker spacing)

      ! Get case info
      call parser%Get('IB IC file',    filename   )

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

      ! Root creates planes
      if (parallel%RankIsRoot()) then

        dl     = 0.5_wp*minval(block%dx)

        ! Plane 1 (bottom)
        center = [ 0.5_wp*(block%pmin(1)+block%pmax(1)), &
                   IBy(1),                               &
                   0.5_wp*(block%pmin(3)+block%pmax(3))]
        normal = '+x2'
        width  = [ block%pmax(1)-block%pmin(1),          &
                   0.0_wp,                               &
                   block%pmax(3)-block%pmin(3)           ]
        vel    = 0.0_wp

        call IB%AddPlane(center,normal,width,vel,dl)

        ! Plane 2 (top)
        center = [ 0.5_wp*(block%pmin(1)+block%pmax(1)), &
                   IBy(2),                               &
                   0.5_wp*(block%pmin(3)+block%pmax(3))]
        normal = '-x2'
        width  = [ block%pmax(1)-block%pmin(1),          &
                   0.0_wp,                               &
                   block%pmax(3)-block%pmin(3)           ]
        vel    = 0.0_wp

        call IB%AddPlane(center,normal,width,vel,dl)
      end if

      ! Treatment for periodicity
      call IB%ApplyPeriodicity

      ! Send to the right rank
      call IB%Communicate

      ! Localize centers and markers on grid
      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), pointer :: BCVal(:,:,:)

      ! 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_SYMMETRY,'ibVF')

      call bcs%SetBC('yL', BC_DIRICHLET,'ibVF')
      call bcs%GetBCPointer('yL','ibVF',BCVal); BCVal=1.0_wp

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

      ! Clear data
      call bcs%Finalize()

      return
    end subroutine SetUpCaseBCS
end program main