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