ginit_noise Subroutine

private subroutine ginit_noise()

Initialise the distribution function with random noise. This is the recommended initialisation option. Each different mode is given a random amplitude between zero and one.

Arguments

None

Contents

Source Code


Source Code

  subroutine ginit_noise
    use species, only: spec, tracer_species
    use theta_grid, only: ntgrid 
    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, is_idx, proc_id
    use dist_fn, only: boundary_option_linked, boundary_option_switch
    use dist_fn, only: l_links, r_links
    use dist_fn, only: pass_right, init_pass_ends
    use redistribute, only: fill, delete_redist
    use ran, only: ranf
    use constant_random, only: init_constant_random
    use constant_random, only: finish_constant_random
    use constant_random, only: constant_ranf=>ranf
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, phit
    real :: a, b
    integer :: iglo, ig, ik, it, il, is, nn


    !CMR: need to document tracer phit parameter   ;-)
    call zero_array(phit)
    do it=2,ntheta0/2+1
       nn = it-1
! extra factor of 4 to reduce the high k wiggles for now
       phit (:, it, 1) = (-1)**nn*exp(-8.*(real(nn)/ntheta0)**2)
    end do

    
    ! keep old (it, ik) loop order to get old results exactly: 

    !Fill phi with random (complex) numbers between -0.5 and 0.5
    if (constant_random_flag) call init_constant_random
    do it = 1, ntheta0
       do ik = 1, naky
          do ig = -ntgrid, ntgrid
            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
             phi(ig,it,ik) = cmplx(a,b)
           end do
!CMR,28/1/13: 
! clean_init debrutalises influence of chop_side on phi with linked bc
          if (chop_side) then
             if (left) then
                if (clean_init.and.(boundary_option_switch .eq. boundary_option_linked)) then
                   !Does this tube have more links to the right than the left?
                   !If so then it's on the left so set to zero
                   if (l_links(it, ik) .lt. r_links(it, ik)) then
                      phi(:,it,ik) = 0.0 
                   !This tube is in the middle so only set the left half to 0
                   else if (l_links(it, ik) .eq. r_links(it, ik)) then
                      phi(:-1,it,ik) = 0.0
                   endif
                else
                   phi(:-1,it,ik) = 0.0
                endif
             else
                if (clean_init.and.(boundary_option_switch .eq. boundary_option_linked)) then
                   !Does this tube have more links to the left than the right?
                   !If so then it's on the right so set to zero
                   if (r_links(it, ik) .lt. l_links(it, ik)) then
                      phi(:,it,ik) = 0.0
                   !This tube is in the middle so only set the right half to 0
                   else if (l_links(it, ik) .eq. r_links(it, ik)) then
                      phi(1:,it,ik) = 0.0
                   endif
                else
                   phi(1:,it,ik) = 0.0
                endif
             endif
          end if
       end do
    end do
    if (constant_random_flag) call finish_constant_random
    
    !Wipe out all but one kx if requested
    if (ikx_init  > 0) call single_initial_kx(phi)
    
    !Sort out the zonal/self-periodic modes
    if (naky .ge. 1 .and. aky(1) == 0.0) then
       !Apply scaling factor
       phi(:,:,1) = phi(:,:,1)*zf_init
       
       !Set ky=kx=0.0 mode to zero in amplitude
       phi(:,1,1) = 0.0
       
       !CMR, 25/1/13: clean_init forces periodic zonal flows
       if (clean_init .and. boundary_option_switch.eq.boundary_option_linked) then
          phi(ntgrid,:,1)=phi(-ntgrid,:,1)
          phit(ntgrid,:,1)=phit(-ntgrid,:,1)
       end if
    end if

    !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))
          phit(:,it+(ntheta0+1)/2,1) = conjg(phit(:,(ntheta0+1)/2+1-it,1))
       enddo
    end if

    !Now set g using data in 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)
       is = is_idx(g_lo,iglo)

       !Handle tracer_species 
       if (spec(is)%type == tracer_species) then          
          g(:,1,iglo) =-phit(:,it,ik)*spec(is)%z*phiinit
       else
          g(:,1,iglo) = -phi(:,it,ik)*spec(is)%z*phiinit
       end if

       !Make sure there's no g in the forbidden regions
       where (forbid(:,il)) g(:,1,iglo) = 0.0
     
       !Set incoming and outgoing boundary conditions
       if ( clean_init .and. boundary_option_switch .eq. boundary_option_linked .and. ik .gt. 1) then
          !CMR: clean_init sets g=0 for incoming particles, 
          !                         and outgoing particles
          !                      as initial condition sets g(:,1,:) = g(:,2,:) )
          !If no links to the left then we're at the left boundary
          if ( l_links(it, ik) .eq. 0 ) g(-ntgrid,1,iglo)=0.0
          !If no links to the right then we're at the right boundary
          if ( r_links(it, ik) .eq. 0 ) g(ntgrid,1,iglo)=0.0
       endif
    end do

    
    !If clean_init is set and there are any links/connections then make sure repeated points are consistent
    if ( clean_init .and. boundary_option_switch .eq. boundary_option_linked .and. sum(r_links+l_links) .gt. 0 ) then
       !   
       !CMR, 23/1/13:  
       !  clean_init: ensure continuity in || direction with linked bcs
       !  set up redistribute type pass_right 
       !  this passes g(ntgrid,1,iglo) to g(-ntgrid,1,iglo_r) on right neighbour 

       !Initialise communications object !NOTE This should be moved to dist_fn in future
       call init_pass_ends(pass_right,'r',1,'c')

       !Pass right boundaries (ntgrid) of g to linked left boundaries (-ntgrid) of g
       call fill(pass_right,g,g)

       !Deallocate comm object as not used anywhere else !NOTE This will need removing if/when initialisation moved to dist_fn
       call delete_redist(pass_right)

       !Make leftwards and rightwards particles the same
       g(:,2,:)=g(:,1,:)

       !Set gnew to be the "fixed" data
       call copy(g, gnew)
    else 
! finally initialise isign=2
       g(:,2,:) = g(:,1,:)
       call copy(g, gnew)
    endif

  end subroutine ginit_noise