main Program

program main
  !>--------------------------------------------------------------------------
  ! Program: Lamb Oseen
  ! Author: Mohamed Houssem Kasbaoui
  !
  ! Solver: CDIFS
  !
  ! Description: Simulation of 4 counter-rotating Lamb-Oseen vortices in a
  ! fully periodic box.
  !
  ! References:
  ! Shuai, Shuai, and M. Houssem Kasbaoui. 2022. “Accelerated Decay of a
  ! Lamb–Oseen Vortex Tube Laden with Inertial Particles in Eulerian–Lagrangian
  ! Simulations.” Journal of Fluid Mechanics 936.
  ! --------------------------------------------------------------------------
  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 initial fields
  call SetUpCaseFields()

  ! 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([.false.,.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 fields.
      implicit none
      ! Work variables
      type(Eulerian_set)   :: fields
      type(eulerian_obj_r) :: V(3)
      type(eulerian_obj_r) :: P
      character(str64)     :: filename
      real(wp)             :: L(3)
      real(wp)             :: radius
      real(wp)             :: circ
      real(wp)             :: rho
      real(wp)             :: mu
      real(wp)             :: Rey
      real(wp)             :: r
      real(wp)             :: sgn
      real(wp)             :: c(2)
      real(wp)             :: coef
      integer              :: i,j,k
      real(wp),parameter   :: twoPi=8.0_wp*atan(1.0_wp)

      ! Get info from parser
      call parser%Get("Fields IC file",  filename)
      call parser%Get("Domain size",     L       )
      call parser%Get("Vortex radius",   radius  )
      call parser%Get("Circulation",     circ    )
      call parser%Get("Fluid density",   rho     )
      call parser%Get("Fluid viscosity", mu      )

      ! 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   )


      associate (lo => block%lo, hi=> block%hi,          &
          x =>block%x , y =>block%y,  z => block%z, &
          xm=>block%xm, ym=>block%ym, zm=> block%zm)

        c = [0.5_wp*(block%pmax(1)+block%pmin(1)),0.5_wp*(block%pmax(2)+block%pmin(2))]

        sgn=(-1)
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)
              ! U component
              r = sqrt((x(i)-c(1))**2+(ym(j)-c(2))**2)
              coef = sgn*circ/twoPi/(r+epsilon(1.0_wp))*(1.0_WP-exp(-r**2/radius**2))
              V(1)%cell(i,j,k) = V(1)%cell(i,j,k) - coef*(ym(j)-c(2))/(r+epsilon(1.0_wp))

              ! V component
              r = sqrt((xm(i)-c(1))**2+(y(j)-c(2))**2)
              coef = sgn*circ/twoPi/(r+epsilon(1.0_wp))*(1.0_WP-exp(-r**2/radius**2))
              V(2)%cell(i,j,k) = V(2)%cell(i,j,k) + coef*(xm(i)-c(1))/(r+epsilon(1.0_wp))
            end do
          end do
        end do

        V(3) = 0.0_wp
        P    = 0.0_wp
      end associate

      ! Write some info to stdout
      if (parallel%RankIsRoot()) then
        ! Reynolds number
        Rey = (circ/twoPi/radius)*radius*rho/(mu+epsilon(1.0_wp))
        write(*,*) "Reynolds number = ", Rey
        write(*,*) "Vortex time     = ", twoPi*radius**2/circ
        write(*,*) "radius/dx       = ", radius/(block%dx(1))
        write(*,*) "radius/dy       = ", radius/(block%dx(2))
        write(*,*) "radius/dz       = ", radius/(block%dx(3))
      end if

      ! 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 SetUpCaseBCS()
      !> Defines boundary conditions.
      use leapBC
      implicit none
      ! Work variables
      type(bc_set) :: bcs
      real(wp)     :: xhi(3)
      real(wp)     :: xlo(3)

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

      associate (pmin => block%pmin, pmax => block%pmax)
        ! - Left boundary
        xlo = [pmin(1),pmin(2),pmin(3)]
        xhi = [pmin(1),pmax(2),pmax(3)]
        call bcs%Add('xL',  xlo, xhi, normal = '+x1')

        ! - xR boundary
        xlo = [pmax(1),pmin(2),pmin(3)]
        xhi = [pmax(1),pmax(2),pmax(3)]
        call bcs%Add('xR',  xlo, xhi, normal = '-x1')

        ! - yL boundary
        xlo = [pmin(1),pmin(2),pmin(3)]
        xhi = [pmax(1),pmin(2),pmax(3)]
        call bcs%Add('yL', xlo, xhi, normal = '+x2')

        ! - yR boundary
        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, 'V1' )
      call bcs%SetBC('yR', BC_SYMMETRY, 'V2' )
      call bcs%SetBC('yR', BC_SYMMETRY, 'V3' )
      call bcs%SetBC('yR', BC_SYMMETRY, 'P'  )
      call bcs%SetBC('yR', BC_SYMMETRY, 'dP' )

      call bcs%SetBC('yL', BC_SYMMETRY, 'V1' )
      call bcs%SetBC('yL', BC_SYMMETRY, 'V2' )
      call bcs%SetBC('yL', BC_SYMMETRY, 'V3' )
      call bcs%SetBC('yL', BC_SYMMETRY, 'P'  )
      call bcs%SetBC('yL', BC_SYMMETRY, 'dP' )

      call bcs%SetBC('xL', BC_SYMMETRY, 'V1' )
      call bcs%SetBC('xL', BC_SYMMETRY, 'V2' )
      call bcs%SetBC('xL', BC_SYMMETRY, 'V3' )
      call bcs%SetBC('xL', BC_SYMMETRY, 'P'  )
      call bcs%SetBC('xL', BC_SYMMETRY, 'dP' )

      call bcs%SetBC('xR', BC_SYMMETRY, 'V1' )
      call bcs%SetBC('xR', BC_SYMMETRY, 'V2' )
      call bcs%SetBC('xR', BC_SYMMETRY, 'V3' )
      call bcs%SetBC('xR', BC_SYMMETRY, 'P'  )
      call bcs%SetBC('xR', BC_SYMMETRY, 'dP' )

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

      ! Clear data
      call bcs%Finalize()

      return
    end subroutine SetUpCaseBCS
end program main