ginit_random_sine Subroutine

private subroutine ginit_random_sine()

Construct the initial condition from a sum of sine waves scaled with random amplitude up to phiinit (ky > 0) or zf_init (ky == 0). Use mode numbers up to minimum of max_mode or last resolved wave (i.e. ntheta/8).

A sine wave is selected in order to ensure that the initial distribution satisfies the parallel boundary conditions (i.e. that either gnew or h = 0).

Arguments

None

Contents

Source Code


Source Code

  subroutine ginit_random_sine
    use theta_grid, only: ntgrid, ntheta, theta
    use kt_grids, only: naky, ntheta0, aky, reality
    use le_grids, only: forbid
    use dist_fn_arrays, only: g, gnew
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx
    use constant_random, only: init_constant_random, finish_constant_random
    use mp, only: proc0, broadcast
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi
    real, dimension(-ntgrid:ntgrid) :: sin_theta_mode
    real :: a, b
    real, parameter :: sqrt_two = sqrt(2.0)
    integer :: iglo, ik, it, il, imode, the_max_mode

    call zero_array(phi)

    ! Make sure the max mode is at least 1, but no bigger than the smaller
    ! of the input max_mode and ntheta/8
    the_max_mode = max(1, min(max_mode, ntheta/8))

    ! We choose to only initialise on proc0 and then broadcast to make
    ! sure that all procs have the same values. If done on a processor
    ! by processor basis, then it's possible that different processors
    ! would produce different values. This could lead to a non-Maxwellian
    ! perturbation (i.e. something not constant in our velocity space).
    if (proc0) then
       if (constant_random_flag) call init_constant_random
       ! Here we go from 2*min_mode (min_mode == 1) to 2*the_max_mode
       ! but then divide imode by two in order to allow us to capture
       ! n + 1/2 mode numbers and hence provide an even parity
       ! contribution to the initial condition.
       do imode = 2, the_max_mode*2
          sin_theta_mode = sin(imode*(theta-theta(-ntgrid))/2.0)
          do ik = 1, naky
             do it = 1, ntheta0
                call set_a_and_b_uniform_random_weighting(a, b)
                phi(:, it, ik) = phi(:, it, ik) + cmplx(a*sin_theta_mode, b*sin_theta_mode)
             end do
          end do
       end do
       if (constant_random_flag) call finish_constant_random

       !Apply reality condition (i.e. -kx mode is conjugate of +kx mode)
       if (reality) then
          do it = 1, ntheta0/2
             phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1))
          enddo
       end if

    end if

    call broadcast(phi)

    ! Now set the distribution function from "phi".
    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo, iglo)
       it = it_idx(g_lo, iglo)
       il = il_idx(g_lo, iglo)
       if (aky(ik) == 0.0) then
          g(:, 1, iglo) = zf_init * phiinit * phi(:, it, ik)

          ! Zero the ky=kx=0 mode
          if (it == 1) g(:, 1, iglo) = 0.0
       else
          g(:, 1, iglo) = phiinit * phi(:, it, ik)
       end if

       where (forbid(:, il)) g(:, 1, iglo) = 0.0
       g(:, 2, iglo) = g(:, 1, iglo)
    end do
    call copy(g, gnew)
  contains
    !> Real and imaginary weights are uniformly distributed between
    !> -1/sqrt(2) and 1/sqrt(2) using either the main random number
    !> generator or the constant random one.
    subroutine set_a_and_b_uniform_random_weighting(a, b)
      use ran, only: ranf
      use constant_random, only: constant_ranf=>ranf
      implicit none
      real, intent(out) :: a, b
      if (constant_random_flag) then
         a = (constant_ranf() - 0.5)
         b = (constant_ranf() - 0.5)
      else
         a = (ranf() - 0.5)
         b = (ranf() - 0.5)
      end if
      a = a * sqrt_two
      b = b * sqrt_two
    end subroutine set_a_and_b_uniform_random_weighting

  end subroutine ginit_random_sine