main Program

program main
  !>--------------------------------------------------------------------------
  ! Program: Sphere Collisions
  ! Author: Mohamed Houssem Kasbaoui
  !
  ! Solver: GRANS
  !
  ! Description: Performs various collisions between a spherical resolved
  ! particle (RP) and immersed boundaries or other resolved particles.
  !
  ! 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
  character(len=str64) :: configuration                                        !! Configuration
  logical              :: use_RP                                               !! Switch to activate Resolved Particles
  logical              :: use_IB                                               !! Switch to activate Immersed Boundaries
  integer              :: CASE_ID                                              !! Case to run
  !> Configurations
  integer, parameter :: CASE_NORMAL_WALL  = 1
  integer, parameter :: CASE_NORMAL_PAIR  = 2
  integer, parameter :: CASE_NORMAL_IB    = 3
  integer, parameter :: CASE_CHANNEL_WALL = 4
  integer, parameter :: CASE_CHANNEL_IB   = 5

  ! Initialize parser
  call parser%Initialize()

  ! Parse input file
  call parser%ParseFile()

  ! Initialize parallel environment
  call parallel%Initialize()

  ! Chose configuration to run
  call parser%Get('Configuration',  configuration, default = "normal wall" )
  select case (trim(adjustl(configuration)))
  case ('normal wall')
    CASE_ID = CASE_NORMAL_WALL
    use_RP  = .true.
    use_IB  = .false.
  case ('normal pair')
    CASE_ID = CASE_NORMAL_PAIR
    use_RP  = .true.
    use_IB  = .false.
  case ('normal IB'  )
    CASE_ID = CASE_NORMAL_IB
    use_RP  = .true.
    use_IB  = .true.
  case ('channel wall')
    CASE_ID = CASE_CHANNEL_WALL
    use_RP  = .true.
    use_IB  = .false.
  case ('channel IB' )
    CASE_ID = CASE_CHANNEL_IB
    use_RP  = .true.
    use_IB  = .true.
  case default
    call parallel%Stop("Unknown case: "//trim(adjustl(configuration)))
  end select

  ! Set the block info
  call SetUpCaseBlock()

  ! Activate optional components
  if (use_RP) call SetUpCaseRP()
  if (use_IB) 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),xlo(3),xhi(3)
      integer          :: N(3), ngc, ilo(3),ihi(3), Nb(3)
      integer          :: dir
      integer          :: i,j,k
      real(wp)         :: r

      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(1), 0.0_wp,-0.5_wp*L(3)]
      xhi = [ 0.5_wp*L(1), 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([.true.,.false.,.true.])

      select case (CASE_ID)
      case (CASE_NORMAL_WALL, CASE_NORMAL_PAIR,CASE_NORMAL_IB)

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

      case (CASE_CHANNEL_WALL,CASE_CHANNEL_IB)

        ! Create a grid, locally refined at the wall
        ! ------------------------------------------
        call parser%Get("Stretching",  r  )

        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

          ! Stretched y-axis
          do j=lo(2),hi(2)+1
            y(j)= (j-lo(2))*L(2)/real(hi(2)-lo(2)+1,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

          ! 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

      end select

      ! 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
      type(ResPart_set) :: RP                                                  !! Resolved particles
      real(wp)          :: diam
      real(wp)          :: rhop(2)
      real(wp)          :: dl
      real(wp)          :: center(3,2)
      real(wp)          :: velp(3,2)
      character(str64)  :: filename
      integer           :: n

      select case (CASE_ID)
      case (CASE_NORMAL_WALL,CASE_NORMAL_IB,CASE_CHANNEL_WALL,CASE_CHANNEL_IB)
        call parser%Get('RP IC file',          filename    )
        call parser%Get("Particle diameter",   diam        )
        call parser%Get("Particle density",    rhop(1)     )
        call parser%Get("Particle center",     center(:,1) )
        call parser%Get("Particle velocity",   velp(:,1),  &
                           default = [0.0_wp,0.0_wp,0.0_wp])

        RP%count = 1
      case (CASE_NORMAL_PAIR)
        call parser%Get('RP IC file',          filename    )
        call parser%Get("Particle diameter",   diam        )
        call parser%Get("Particle 1 density",  rhop(1)     )
        call parser%Get("Particle 2 density",  rhop(2)     )
        call parser%Get("Particle 1 center",   center(:,1) )
        call parser%Get("Particle 2 center",   center(:,2) )
        call parser%Get("Particle 1 velocity", velp(:,1),  &
                           default = [0.0_wp,0.0_wp,0.0_wp])
        call parser%Get("Particle 2 velocity", velp(:,2),  &
                           default = [0.0_wp,0.0_wp,0.0_wp])

        RP%count = 2
      end select

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

      ! Seed particles on root
      if (parallel%RankIsRoot()) then
        ! Activate particle
        call RP%Resize(RP%count)

        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
            ! Position
            particle(n)%p    = center(:,n)
            ! Velocity
            particle(n)%v    = velp(:,n)
            ! Angular velocity
            particle(n)%w    = 0.0_wp
            ! Density
            particle(n)%rho  = rhop(n)
            ! 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 SetUpCaseIB()
      !> Builds and writes immersed boundary.
      use immersed_boundaries
      implicit none
      ! Work variables
      character(str64) :: filename
      type(marker_set) :: IB                                                   !! IB
      real(wp)         :: IBy(2)                                               !! IB plane vertical position
      real(wp)         :: center(3,2)                                          !! IB plane centroid
      character(len=3) :: normal(2)                                            !! IB plane normal direction
      real(wp)         :: width(3,2)                                           !! IB plane width
      real(wp)         :: vel(3,2)                                             !! IB plane velocity
      real(wp)         :: dl                                                   !! IB plane resolution (marker spacing)

      select case (CASE_ID)
      case (CASE_NORMAL_IB)
        call parser%Get('IB IC file',    filename   )
        call parser%Get("IB 1 position", IBy(1)    )
      case (CASE_CHANNEL_IB)
        call parser%Get('IB IC file',    filename   )
        call parser%Get("IB 1 position", IBy(1)    )
        call parser%Get("IB 2 position", IBy(2)    )
      end select

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

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

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

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

        ! Plane 2 (top)
        if (CASE_ID.eq.CASE_CHANNEL_IB) then
          center(:,2) = [ 0.5_wp*(block%pmin(1)+block%pmax(1)), &
                          IBy(2),                               &
                          0.5_wp*(block%pmin(3)+block%pmax(3))]
          normal(2)   = '-x2'
          width(:,2)  = [ block%pmax(1)-block%pmin(1),          &
                          0.0_wp,                               &
                          block%pmax(3)-block%pmin(3)           ]
          vel(:,2)    = 0.0_wp
          dl          = 0.5_wp*minval(block%dx)

          call IB%AddPlane(center(:,2),normal(2),width(:,2),vel(:,2),dl)
        end if
      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
    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('yL', BC_DIRICHLET,'ibVF')
      call bcs%SetBC('yR', BC_DIRICHLET,'ibVF')

      select case (CASE_ID)
      case (CASE_NORMAL_WALL, CASE_NORMAL_PAIR,CASE_CHANNEL_WALL)
        call bcs%GetBCPointer('yL','ibVF',BCVal); BCVal=0.0_wp
        call bcs%GetBCPointer('yR','ibVF',BCVal); BCVal=0.0_wp
      case (CASE_NORMAL_IB)
        call bcs%GetBCPointer('yL','ibVF',BCVal); BCVal=1.0_wp
        call bcs%GetBCPointer('yR','ibVF',BCVal); BCVal=0.0_wp
      case (CASE_CHANNEL_IB)
        call bcs%GetBCPointer('yL','ibVF',BCVal); BCVal=1.0_wp
        call bcs%GetBCPointer('yR','ibVF',BCVal); BCVal=1.0_wp
      end select

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

      ! Clear data
      call bcs%Finalize()

      return
    end subroutine SetUpCaseBCS
end program main