main Program

program main
  !>--------------------------------------------------------------------------
  ! Program: Rebound
  ! Author: Mohamed Houssem Kasbaoui
  !
  ! Solver: CDIFS
  !
  ! Description:  A fully resolved sphere that rebounds on a planar IB
  !
  ! 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
  character(len=str64) :: configuration                                        !! Configuration
  logical              :: use_IB                                               !! Switch to activate Immersed Boundaries
  integer              :: CASE_ID                                              !! Case to run
  !> Configurations
  integer, parameter :: CASE_IB_COLLISION   = 1
  integer, parameter :: CASE_WALL_COLLISION = 2

  ! 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 = "wall" )
  select case (trim(adjustl(configuration)))
  case ('WALL')
    CASE_ID = CASE_WALL_COLLISION
    use_IB  = .false.
  case ('IB')
    CASE_ID = CASE_IB_COLLISION
    use_IB  = .true.
  case default
    call parallel%Stop("Unknown case: "//trim(adjustl(configuration)))
  end select

  ! Set the block info
  call SetUpCaseBlock()

  ! Set the initial fields
  call SetUpCaseFields()

  ! Set resolved particles
  call SetUpCaseResPart()

  ! Set the immersed boundary
  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)
      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(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.])

      ! 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 SetUpCaseResPart()
      !> 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)          :: xc(3)
      real(wp)          :: dl
      integer           :: n

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

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

      if (parallel%RankIsRoot()) then

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

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

        select type(particle => RP%p)
        type is (ResPart_obj)
          do n=1,1
            ! Particle globabl ID
            particle(n)%id   = int(n,kind=8)
            ! Diameter
            particle(n)%d    = diam
            ! Position
            particle(n)%p    = xc
            ! Velocity
            particle(n)%v    = 0.0_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 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 SetUpCaseResPart
    subroutine SetUpCaseIB
      !> Builds and writes immersed boundary.
      use immersed_boundaries
      implicit none
      type(marker_set) :: IB
      character(str64) :: filename
      character(len=3) :: normal
      real(wp)         :: xc(3)
      real(wp)         :: xc_default(3)
      real(wp)         :: width(3)
      real(wp)         :: dl

      ! Default position (if not supplied by user in input file)
      xc_default = [ 0.5_wp*(block%pmin(1)+block%pmax(1)), &
                     block%pmin(2)+4.0_wp*block%dx(2),     &
                     0.5_wp*(block%pmin(3)+block%pmax(3)) ]

      ! Get info from parser
      call parser%Get("IB IC file",   filename                    )
      call parser%Get("IB position",  xc,     default = xc_default)
      call parser%Get("IB normal",    normal, default = "+x2"     )

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

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

      ! Determine plane width to cover domain from end to end
      width = block%pmax-block%pmin
      select case (trim(adjustl(normal)))
      case ("+x1", "-x1")
        width(1) = 0.0_wp
      case ("+x2", "-x2")
        width(2) = 0.0_wp
      case ("+x3", "-x3")
        width(3) = 0.0_wp
      case default
        call parallel%Stop("Invalid plane normal")
      end select

      if (parallel%RankIsRoot()) then
        ! Add plane
        call IB%AddPlane(xc,normal,width,[0.0_wp,0.0_wp,0.0_wp],dl)
      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)
      real(wp)          :: 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

      select case (CASE_ID)
      case (CASE_WALL_COLLISION)
        call bcs%SetBC('yL', BC_WALL )
        call bcs%SetBC('yR', BC_WALL )

        call bcs%SetBC('yR', BC_SYMMETRY, 'ibVF')
        call bcs%SetBC('yL', BC_SYMMETRY, 'ibVF')

      case (CASE_IB_COLLISION)
        call bcs%SetBC('yL', BC_OUTFLOW)
        call bcs%SetBC('yR', BC_OUTFLOW)

        call bcs%SetBC('yL', BC_DIRICHLET,'ibVF')
        call bcs%SetBC('yR', BC_DIRICHLET,'ibVF')

        call bcs%GetBCPointer('yL','ibVF',BCVal); BCVal=1.0_wp
        call bcs%GetBCPointer('yR','ibVF',BCVal); BCVal=0.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