main Program

program main
  !>--------------------------------------------------------------------------
  ! Program: Vortex impingement
  ! Author: Mohamed Houssem Kasbaoui
  !
  ! Solver: CDIFS
  !
  ! Description: Vortex dipole impinging on bottom wall
  !
  ! 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

  ! 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()
      !> Setup and write the block file
      implicit none
      ! Work variables
      character(str64) :: filename
      real(wp)         :: L(3)
      integer          :: N(3)
      integer          :: ngc
      integer          :: Nb(3)
      real(wp)         :: r
      integer          :: i,j,k,dir

      ! 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("Stretching factor", r       )
      call parser%Get("Partition",         Nb      )

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

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

      ! Create a grid, locally refined at the 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

        ! Stretched y-axis
        do j=lo(2),hi(2)+1
          y(j)= (j-lo(2))*L(2)/real(hi(2)-lo(2)+1,wp)
          y(j)= L(2)*(1.0_wp - tanh(r*(L(2)-y(j))/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)
        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 SetUpCaseFields()
      !> Set the 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
      real(wp)             :: L(3)
      real(wp)             :: r
      real(wp)             :: a
      real(wp)             :: Vd
      real(wp)             :: Vr, Vt
      real(wp)             :: costheta
      real(wp)             :: sintheta
      real(wp)             :: center(3)
      real(wp)             :: rho
      real(wp)             :: mu
      real(wp)             :: Rey
      integer              :: i,j,k
      real(wp),parameter   :: twoPi=8.0_wp*atan(1.0_wp)
      ! Bessel firt root and constant C
      real(wp),parameter   :: a1 = 3.8317059702075123115_wp
      real(wp)             :: C

      ! Get info from parser
      call parser%Get("Fields IC file",  filename)
      call parser%Get("Domain size",     L       )
      call parser%Get('Dipole radius',   a       )
      call parser%Get('Dipole velocity', Vd      )
      call parser%Get('Dipole center',   center  )
      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)

        V(1) = 0.0_wp
        V(2) = 0.0_wp
        V(3) = 0.0_wp
        P    = 0.0_wp

        ! Compute the constant
        C = 2.0_wp/J0(a1)

        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)
              ! U - velocity
              r = sqrt((x(i)-center(1))**2+(ym(j)-center(2))**2)
              costheta= (ym(j)-center(2))/(r+epsilon(1.0_wp))
              sintheta= (x (i)-center(1))/(r+epsilon(1.0_wp))

              if ((r/a).le.1.0_wp) then
                r = a1*r/a
                Vr = Vd * ( C*J1(r)/(r+epsilon(1.0_wp)) - 1.0_WP ) * costheta
                Vt = Vd * ( 1.0_WP - C*(J0(r)-J1(r)/(r+epsilon(1.0_wp))) ) * sintheta
              else
                Vr = - Vd/(r/a)**2 * costheta
                Vt = - Vd/(r/a)**2 * sintheta
              end if
              V(1)%cell(i,j,k) = Vr*sintheta+Vt*costheta

              ! V - velocity
              r = sqrt((xm(i)-center(1))**2+(y(j)-center(2))**2)
              costheta= (y (j)-center(2))/(r+epsilon(1.0_wp))
              sintheta= (xm(i)-center(1))/(r+epsilon(1.0_wp))


              if ((r/a).le.1.0_wp) then
                r = a1*r/a
                Vr = Vd * ( C*J1(r)/(r+epsilon(1.0_wp)) - 1.0_WP ) * costheta
                Vt = Vd * ( 1.0_WP - C*(J0(r)-J1(r)/(r+epsilon(1.0_wp))) ) * sintheta
              else
                Vr = - Vd/(r/a)**2 * costheta
                Vt = - Vd/(r/a)**2 * sintheta
              end if
              V(2)%cell(i,j,k) = Vr*costheta-Vt*sintheta
            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 = rho*Vd*a/(mu+epsilon(1.0_wp))
        write(*,*) "Reynolds number = ", Rey
      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()
      !> Setup the boundary conditions
      use leapBC
      implicit none
      ! Work variables
      type(bc_set) :: bcs
      real(wp)     :: xhi(3),xlo(3)

      ! 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_OUTFLOW)
      call bcs%SetBC('yL', BC_WALL   )

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

      ! Clear data
      call bcs%Finalize()

      return
    end subroutine SetUpCaseBCS
    pure function J0(x) result(val)
      !> Bessel function J0(x)
      implicit none
      real(wp), intent(in) :: x
      real(wp)             :: val
      ! work variables
      real(wp) :: ax,xx,z,y
      real(wp), parameter :: r(6) = [ 57568490574.0_wp,-13362590354.0_wp, 651619640.7_wp,&
                                     -11214424.18_wp,   77392.33017_wp,  -184.9052456_wp]
      real(wp), parameter :: s(6) = [ 57568490411.0_wp, 1029532985.0_wp, 9494680.718_wp, &
                                      59272.64853_wp,   267.8532712_wp,  1.0_wp]
      real(wp), parameter :: p(5) = [ 1.0_wp,-0.1098628627E-2_wp, 0.2734510407E-4_wp,&
                                             -0.2073370639E-5_wp, 0.2093887211E-6_wp]
      real(wp), parameter :: q(5) = [-.1562499995E-1_wp,.1430488765E-3_wp,-.6911147651E-5_wp,&
                                      .7621095161E-6_wp,-.934945152E-7_wp]

      if (abs(x).lt.8.0_wp) then
         y   = x**2
         val = (r(1)+y*(r(2)+y*(r(3)+y*(r(4)+y*(r(5)+y*r(6)))))) / &
               (s(1)+y*(s(2)+y*(s(3)+y*(s(4)+y*(s(5)+y*s(6))))))
      else
         ax  = abs(x)
         z   = 8.0_wp/ax
         y   = z**2
         xx  = ax-.785398164_wp
         val = sqrt(.636619772_wp/ax)*(cos(xx)*(p(1)+y*(p(2)+y*(p(3)+y*(p(4)+y*p(5))))) - &
               z*sin(xx)*(q(1)+y*(q(2)+y*(q(3)+y*(q(4)+y*q(5))))))
      end if
      return
    end function J0
    pure function J1(x) result(val)
      !> Bessel function J1(x)
      implicit none
      real(wp), intent(in) :: x
      real(wp)             :: val
      ! work variables
      real(wp) :: ax,xx,z,y
      real(wp), parameter :: r(6) = [ 72362614232.0_wp,-7895059235.0_wp, 242396853.1_wp, &
                                     -2972611.439_wp,   15704.48260_wp, -30.16036606_wp]
      real(wp), parameter :: s(6) = [144725228442.0_wp, 2300535178.0_wp, 18583304.74_wp, &
                                     99447.43394_wp,    376.9991397_wp,  1.0_wp]
      real(wp), parameter :: p(5) = [1.0_wp, 0.183105E-2_wp,     -0.3516396496E-4_wp, &
                                             0.2457520174E-5_wp, -0.240337019E-6_wp]
      real(wp), parameter :: q(5) = [ .04687499995_wp,-.2002690873E-3_wp,.8449199096E-5_wp,&
                                     -.88228987E-6_wp,.105787412E-6_wp]

      if (abs(x).lt.8.0_wp) then
         y   = x**2
         val = x*(r(1)+y*(r(2)+y*(r(3)+y*(r(4)+y*(r(5)+y*r(6)))))) / &
               (s(1)+y*(s(2)+y*(s(3)+y*(s(4)+y*(s(5)+y*s(6))))))
      else
         ax  = abs(x)
         z   = 8.0_wp/ax
         y   = z**2
         xx  = ax-2.356194491_wp
         val = sqrt(.636619772_wp/ax)*(cos(xx)*(p(1)+y*(p(2)+y*(p(3)+y*(p(4)+y*p(5))))) - &
               z*sin(xx)*(q(1)+y*(q(2)+y*(q(3)+y*(q(4)+y*q(5))))))*sign(1.0_wp,x)
      end if
      return
    end function J1
end program main