main Program

program main
  !>--------------------------------------------------------------------------
  ! Program: Homogeneous Isotropic Turbulence
  ! Author: Mohamed Houssem Kasbaoui
  !
  ! Solver: CDIFS
  !
  ! Description: Generates synthetic homogeneous isotropic turbulence using
  ! the following energy spectrum:
  !     E(k) = A k^4 exp(-2(k/kp)^2)
  ! where A is an amplitude determined from the box size and target
  ! Taylor-microscale Reynolds number, and kp is the peak wavenumber in
  ! the box (associated by default to mode 4).
  !
  ! NOTE: Due to FFTW constraints, this initialization must be performed
  ! in serial.
  !
  ! References:
  ! R. S. Rogallo, Numerical experiments in homogeneous turbulence,
  ! NASA STI/Recon Technical Report N 81, 31508 (1981).
  !
  ! Rosales, C. & Meneveau, C. Linear forcing in numerical simulations of
  ! isotropic turbulence: Physical space implementations and convergence
  ! properties. Physics of Fluids 17, 095106 (2005).
  ! --------------------------------------------------------------------------
  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
  integer            :: Ng(3)                                                  !! Global number of grid points in each direction


  ! Initialize parser
  call parser%Initialize()

  ! Parse input file
  call parser%ParseFile()

  ! Initialize parallel environment
  call parallel%Initialize()

  if (parallel%nproc.ne.1) &
    call parallel%Stop("HIT initialization must be done in serial")

  ! 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          :: 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", Ng      )
      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=Ng

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

      ! Setup the domain periodicity
      call block%SetPeriodicity([.true.,.true.,.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 fields.
      use iso_fortran_env, only : stdout => output_unit
      implicit none
      ! Work variables
      type(Eulerian_set)   :: fields
      type(eulerian_obj_r) :: V(3)
      type(eulerian_obj_r) :: P
      character(str64)     :: filename
      real(wp)             :: rho
      real(wp)             :: mu
      real(wp)             :: nu
      real(wp)             :: dissipation_rate
      real(wp)             :: Re_lambda
      integer              :: peak_mode
      real(wp)             :: ell
      real(wp)             :: L(3)
      integer              :: Nk(3)
      complex(wp),          &
               allocatable :: Vk(:,:,:,:)
      real(wp)             :: lambda
      real(wp)             :: Vrms
      real(wp)             :: eta_k
      real(wp)             :: dV
      integer              :: i,j,k
      real(wp), parameter  :: C=0.19_wp

      ! Get info from parser
      call parser%Get("Fields IC file",   filename               )
      call parser%Get("Fluid density",    rho                    )
      call parser%Get("Fluid viscosity",  mu                     )
      call parser%Get("Reynolds number",  Re_lambda              )
      call parser%Get("Peak mode",        peak_mode, default = 4 )

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

      ! Kinematic viscosity
      nu = mu/rho

      ! Domain size
      L = block%pmax - block%pmin

      ! Integral length scale based on Rosales & Meneveau, Physics of Fluids, 2005
      ell = C*minval(L)

      ! Dissipation rate
      dissipation_rate = (nu/15.0_wp)**3*Re_lambda**6/ell**4

      ! Velocity fluctuation rms
      Vrms = (dissipation_rate*ell)**(1.0_wp/3.0_wp)

      call BuildFourierModes(Vrms,peak_mode,Nk,Vk)
      call InvertFourierModes(Nk,Vk,V)

      deallocate(Vk)

      ! Rescale
      Vrms = 0.0_wp
      do k=block%lo(3),block%hi(3)
        do j=block%lo(2),block%hi(2)
          do i=block%lo(1),block%hi(1)
            dV = (block%xm(i+1)-block%xm(i)) &
               * (block%ym(j+1)-block%ym(j)) &
               * (block%zm(k+1)-block%zm(k))
            Vrms = Vrms +(V(1)%cell(i,j,k)**2 &
                        + V(2)%cell(i,j,k)**2 &
                        + V(3)%cell(i,j,k)**2)*dV
          end do
        end do
      end do
      Vrms = sqrt((1.0_wp/3.0_wp)/(L(1)*L(2)*L(3))*Vrms)

      V(1)%cell = V(1)%cell* sqrt(nu*Re_lambda/ sqrt(15.0_wp*nu/dissipation_rate))/Vrms
      V(2)%cell = V(2)%cell* sqrt(nu*Re_lambda/ sqrt(15.0_wp*nu/dissipation_rate))/Vrms
      V(3)%cell = V(3)%cell* sqrt(nu*Re_lambda/ sqrt(15.0_wp*nu/dissipation_rate))/Vrms

      ! Compute urms
      Vrms = 0.0_wp
      do k=block%lo(3),block%hi(3)
        do j=block%lo(2),block%hi(2)
          do i=block%lo(1),block%hi(1)
            dV = (block%xm(i+1)-block%xm(i)) &
               * (block%ym(j+1)-block%ym(j)) &
               * (block%zm(k+1)-block%zm(k))
            Vrms = Vrms +(V(1)%cell(i,j,k)**2 &
                        + V(2)%cell(i,j,k)**2 &
                        + V(3)%cell(i,j,k)**2)*dV
          end do
        end do
      end do
      Vrms = sqrt((1.0_wp/3.0_wp)/(L(1)*L(2)*L(3))*Vrms)

      lambda    = sqrt(15.0_wp*nu/dissipation_rate)*Vrms
      Re_lambda = lambda*Vrms/nu
      eta_k     = (nu**0.75_wp)/(dissipation_rate**0.25_wp)

      write(stdout,*) "Taylor Microscale Reynolds : ", Re_lambda
      write(stdout,*) "Velocity fluctuation rms   : ", Vrms
      write(stdout,*) "Taylor microscale (lambda) : ", lambda
      write(stdout,*) "Kolomogorov length (eta_k) : ", eta_k
      write(stdout,*) "Dissipation rate           : ", dissipation_rate
      write(stdout,*) "Eddy turounover time       : ", Vrms**2/dissipation_rate
      write(stdout,*) "HIT linear forcing         : ", dissipation_rate/(3.0_wp*Vrms**2)
      write(stdout,*) "Resolution (kmax*eta_k)    : ", 4.0_wp*atan(1.0_wp)*minval(Ng/L)*eta_k


      ! 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

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

      ! Fully-periodic, nothing to do

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

      ! Clear data
      call bcs%Finalize()

      return
    end subroutine SetUpCaseBCS
    function GetEnergySpectrum(kmag,kp,A) result (res)
      !> Computes E(k)
      implicit none
      real(wp), intent(in) :: kmag                                             !! Wavevector magnitude
      real(wp), intent(in) :: kp                                               !! Peak
      real(wp), intent(in) :: A                                                !! Spectrum amplitude
      real(wp)             :: res                                              !! Result

      res = A*(kmag**4)*exp(-2.0_wp*(kmag/kp)**2)

      return
    end function GetEnergySpectrum
    subroutine BuildFourierModes(Vrms,peak_mode,Nk,Vk)
      !> Builds the Fourier modes, scaled with the provided
      ! energy spectrum.
      ! Note that Since the FFT of a real field satisfies the Hermitian
      ! symmetry, one needs to describe only half of the wavenumber space.
      ! Direction 1 is truncated in half to conform with FFTW.
      implicit none
      real(wp),     intent(in)  :: Vrms                                        !! Wavevector magnitude
      integer,      intent(in)  :: peak_mode                                   !! Peak
      integer,      intent(out) :: Nk(3)                                       !! Size of fourier modes
      complex(wp), allocatable,  &
                    intent(out) :: Vk(:,:,:,:)                                 !! Fourier modes
      ! Work variables
      real(wp)             :: A
      real(wp)             :: kp
      real(wp)             :: Ek
      real(wp)             :: kmag
      real(wp)             :: phi
      real(wp)             :: kappa(3)
      real(wp)             :: rand(3)
      real(wp)             :: dk(3)
      integer              :: mi,mj,mk
      real(wp),  parameter :: Pi = 4.0_wp*atan(1.0_wp)
      complex(wp),          &
                 parameter :: ii = (0.0_wp,1.0_wp)

      ! Infinitesimal wavenumbers
      dk = 2.0_wp*Pi/(block%pmax-block%pmin)

      ! Peak wavenumber
      kp = real(peak_mode,wp)*maxval(dk)

      ! Spectrum prefactor
      A = 16.0_wp/sqrt(0.5_wp*Pi)*Vrms**2/kp**5

      ! Number of fourier modes in each direction
      Nk = [Ng(1)/2+1, Ng(2), Ng(3)]

      allocate(Vk(3,Nk(1),Nk(2),Nk(3)),mold=(0.0_wp,0.0_wp))

      ! Build Fourier modes following Rogallo's method
      do mk=1,Nk(3)

        ! Wavenumber in 3-dir
        kappa(3) = real(mk-1,wp)*dk(3)
        if (mk.gt.Ng(3)/2+1) kappa(3) = -real(Ng(3)-(mk-1),wp)*dk(3)

        do mj=1,Nk(2)

          ! Wavenumber in 2-dir
          kappa(2) = real(mj-1,wp)*dk(2)
          if (mj.gt.Ng(2)/2+1) kappa(2) = -real(Ng(2)-(mj-1),wp)*dk(2)

          do mi=1,Nk(1)

            ! Wavenumber in 1-dir
            kappa(1) = real(mi-1,wp)*dk(1)

            ! Wavevector magnitude
            kmag=norm2(kappa)

            if (kmag.gt.1.0e-10_wp) then
              ! Energy at this shell
              Ek = GetEnergySpectrum(kmag,kp,A)

              ! Generate a random phase
              call random_number(phi)
              phi    = 2.0_wp*Pi*phi

              ! Generate a complex vector with a random amplitude and the previous random phase
              call random_number(rand)
              Vk(:,mi,mj,mk) = (rand - 0.5_wp) * exp(ii * phi)

              ! Make it orthogonal to the wavevector to satisfy continuity
              Vk(:,mi,mj,mk) = Vk(:,mi,mj,mk) - (1.0_wp/kmag**2)*dot_product(kappa,Vk(:,mi,mj,mk))*kappa

              ! Rescale to match target energy spectrum
              Vk(:,mi,mj,mk) = Vk(:,mi,mj,mk) * sqrt(Ek / (4.0_wp * pi * kmag**2))

            else
              Vk(:,mi,mj,mk) = (0.0_wp,0.0_wp)
            end if

          end do
        end do
      end do

      return
    end subroutine BuildFourierModes
    subroutine InvertFourierModes(Nk,Vk,V)
      ! Computes velocities in physical space from fourier modes.
      use, intrinsic :: iso_c_binding
      implicit none
      integer,      intent(in)    :: Nk(3)                                     !! Size of fourier modes
      complex(wp),  intent(in)    :: Vk(:,:,:,:)                               !! Fourier modes
      type(eulerian_obj_r),        &
                    intent(inout) :: V(3)
      ! Work variables
      type(C_PTR)             :: plan       = C_NULL_PTR
      type(C_PTR)             :: p_in       = C_NULL_PTR
      type(C_PTR)             :: p_out      = C_NULL_PTR
      complex(C_DOUBLE_COMPLEX),  &
                      pointer :: in(:,:,:)  => null()
      real(C_DOUBLE), pointer :: out(:,:,:) => null()
      integer                 :: mi,mj,mk
      integer                 :: i,j,k,dir

      include 'fftw3.f03'

      ! Allocate memory using FFTW to ensure SIMD allignment
      ! Note that the output must include +2 padding cells in the
      ! third dimension. This is a requirement from FFTW.
      p_in  = fftw_alloc_complex(int(Nk(1)*Nk(2)*Nk(3),         C_SIZE_T))
      p_out = fftw_alloc_real(   int(2*(Ng(1)/2+1)*Ng(2)*Ng(3), C_SIZE_T))

      ! Convert C pointers to Fortran arrays
      call c_f_pointer(p_in,  in,  [Nk(1),Nk(2),Nk(3)])
      call c_f_pointer(p_out, out, [Ng(1),Ng(2),Ng(3)])

      ! Create plan - dimensions must be inverted per FFTW.
      plan = fftw_plan_dft_c2r_3d(Ng(3), Ng(2), Ng(1), in, out, FFTW_MEASURE)


      do dir=1,3

        do mk=1,Nk(3)
          do mj=1,Nk(2)
            do mi=1,Nk(1)
               in(mi,mj,mk) = cmplx(Vk(dir,mi,mj,mk),kind=C_DOUBLE_COMPLEX)
            end do
          end do
        end do

        ! Execute plan (build inverse)
        call fftw_execute_dft_c2r(plan, in, out)

        ! Normalize (FFTW does not automatically normalize)
        do k=1,Ng(3)
          do j=1,Ng(2)
            do i=1,Ng(1)
              V(dir)%cell(i,j,k) = real(out(i,j,k),wp)/real(Ng(1)*Ng(2)*Ng(3),wp)
            end do
          end do
        end do

      end do

      ! Clean up
      call dfftw_destroy_plan(plan)
      call fftw_free(p_out)
      call fftw_free(p_in)

      return
    end subroutine InvertFourierModes
end program main