Proper implementation of multiple N_RADIAL capability for ZF test

cgyro/src/cgyro_check.f90

@@ -36,7 +36,7 @@ subroutine cgyro_check
   endif
 
   if (n_species > 11) then
-     call cgyro_error('n_species <= 11.')
+     call cgyro_error('n_species cannot be greater than 11.')
      return
   endif
 

cgyro/src/cgyro_equilibrium.F90

@@ -235,13 +235,11 @@ subroutine cgyro_equilibrium
 
      do itor=nt1,nt2
        do ir=1,n_radial
-        k_perp(ic_c(ir,it),itor) = &
-               sqrt((2.0*pi*(px(ir)+px0)*geo_grad_r(it)/length &
-             + k_theta_base*itor*geo_gq(it)*geo_captheta(it))**2 &
-             + (k_theta_base*itor*geo_gq(it))**2)
         k_x(ic_c(ir,it),itor) = &
                2.0*pi*(px(ir)+px0)*geo_grad_r(it)/length &
              + k_theta_base*itor*geo_gq(it)*geo_captheta(it)
+        k_perp(ic_c(ir,it),itor) = &
+               sqrt(k_x(ic_c(ir,it),itor)**2+(k_theta_base*itor*geo_gq(it))**2)
        enddo
      enddo
 

cgyro/src/cgyro_init_h.f90

@@ -101,16 +101,15 @@ subroutine cgyro_init_h
               ir = ir_c(ic) 
               it = it_c(ic)
 
-              if (is == 1 .and. px(ir) /= 0) then
+              if (is == 1) then
                  arg = k_perp(ic,itor)*rho*vth(is)*mass(is)/(z(is)*bmag(it)) &
                       *vel2(ie)*sqrt(1.0-xi(ix)**2)
                  h_x(ic,iv_loc,itor) = 1e-6*bessel_j0(abs(arg))
 
-                 ! J0 here for the ions is equivalent to having
-                 ! the electrons deviate in density.
-                 ! Alternatively this is the result of instantaneous
-                 ! gyroaveraging after the deposition of particles in
-                 ! a certain k_radial mode.
+                 ! J0 here for the ions is equivalent to having the electrons
+                 ! deviate in density. Alternatively this is the result of
+                 ! instantaneous gyroaveraging after the deposition of particles
+                 ! in a certain k_radial mode.
 
               endif
            enddo

cgyro/src/cgyro_init_kernel.F90

@@ -47,8 +47,8 @@ subroutine cgyro_init_kernel
   ! 1. MPI setup
   call cgyro_mpi_grid
   if (error_status > 0) then
-          call cgyro_final_kernel
-          return
+     call cgyro_final_kernel
+     return
   endif
 
   call system_clock(aftermpi_time,kernel_count_rate,kernel_count_max)

cgyro/src/cgyro_make_profiles.F90

@@ -360,12 +360,14 @@ subroutine cgyro_make_profiles
      ! my_toroidal == 0
      ! k_theta == 0
 
-     call cgyro_info('Triggered zonal flow test')
-
-     if (n_radial /= 1) then
-        call cgyro_info('Zonal flow test with n_radial > 1')
+     if (n_radial == 1) then
+        call cgyro_info('ZF test with n_radial = 1 ; setting UP_RADIAL=0')
+     else
+        call cgyro_info('ZF test with n_radial > 1 ; setting UP_RADIAL=0')
      endif
 
+     up_radial = 0.0
+
   else if (n_toroidal == 1) then
 
      ! Single linear mode (assume n=1, compute rho)
@@ -420,7 +422,7 @@ subroutine cgyro_make_profiles
      omega_eb_base = k_theta_base*length*gamma_e/(2*pi)
      select case (shear_method)
      case (1)
-        call cgyro_info('ExB shear: Hammett discrete shift') 
+        call cgyro_info('ExB shear: Hammett discrete shift (do not use for production simulations)') 
      case (2)
         call cgyro_info('ExB shear: Wavenumber advection') 
         source_flag = 1
@@ -434,7 +436,7 @@ subroutine cgyro_make_profiles
   endif
 
   if (profile_shear_flag == 1) then
-     call cgyro_info('Profile shear: Continuous wavenumber advection') 
+     call cgyro_info('Profile shear: Wavenumber advection') 
      source_flag = 1
   else
      sdlnndr(1:n_species) = 0.0

cgyro/src/cgyro_write_initdata.f90

@@ -39,6 +39,13 @@ subroutine cgyro_write_initdata
           write(io,'(t3,i4,t12,i4,t21,i4,t29,i5,t36,i3)') nc_loc,nv_loc,nsplit,n_proc,n_omp
         endif
      endif
+     write(io,*)
+
+     if (kymax < -99.9) then
+        lfmt = '(a,i4,2x,2(f7.2,2x),2x,f6.2,5x,i4,2a)'
+     else
+        lfmt = '(a,i4,2x,2(f7.3,2x),2x,f6.2,5x,i4,2a)'
+     endif
 
      if (zf_test_mode == 0) then
 
@@ -48,16 +55,9 @@ subroutine cgyro_write_initdata
         else
            kymax = q/rmin*(n_toroidal-1)*rho
         endif
-
-        if (kymax < -99.9) then
-           lfmt = '(a,i4,2x,2(f7.2,2x),2x,f6.2,5x,i4,2a)'
-        else
-           lfmt = '(a,i4,2x,2(f7.3,2x),2x,f6.2,5x,i4,2a)'
-        endif
 
         if (nonlinear_flag == 0) then
 
-           write(io,*)
            write(io,*) '          n    Delta      Max     L/rho'
            write(io,lfmt) ' kx*rho:',&
                 n_radial,2*pi*rho/length,2*pi*rho*(n_radial/2-1)/length,length/rho
@@ -66,8 +66,6 @@ subroutine cgyro_write_initdata
 
         else
 
-
-           write(io,*)
            write(io,*) '          n    Delta      Max     L/rho    n_fft'
            call prime_factors(nx,msg)
            write(io,lfmt) ' kx*rho:',&
@@ -81,14 +79,9 @@ subroutine cgyro_write_initdata
 
         ! ZONAL FLOW TEST ONLY
 
-        if (n_radial==1) then
-           write(io,*)
-           write(io,'(t2,a,1pe10.3)') ' kx*rho: ',2*pi*rho/length
-        else
-           write(io,*) '          n          Delta            Max           L/rho'
-           write(io,'(a,i4,2x,2(g0.8,2x),2x,g0.8)') ' kx*rho:',&
-                n_radial,2*pi*rho/length,2*pi*rho*(n_radial/2-1)/length,length/rho
-        endif
+        write(io,*) '          n    Delta      Max     L/rho'
+        write(io,lfmt) ' kx*rho:',&
+             n_radial,2*pi*rho/length,2*pi*rho*n_radial/length,length/rho
 
      endif
 

f2py/pygacode/cgyro/data.py

@@ -369,8 +369,15 @@ def getgrid(self):
       l=11
 
       self.p = np.array(data[l:l+self.n_radial],dtype=int)
-      # kx*rho (ignore leftmost "special" element)
-      self.kx = 2*np.pi*self.p[1:]/self.length
+
+      if self.p[0] > 0:
+         # zonal flow test
+         self.kx = 2*np.pi*self.p[:]/self.length
+         self.zf_test = True
+      else:
+         # kx*rho (ignore leftmost "special" element)
+         self.kx = 2*np.pi*self.p[1:]/self.length
+         self.zf_test = False
 
       mark = l+self.n_radial
       self.theta = np.array(data[mark:mark+self.n_theta])
@@ -409,10 +416,11 @@ def getgrid(self):
 
       # Construct k_perp
       # NOTE: kx,ky and kperp are kx*rho, ky*rho, kperp*rho
-      nx = self.n_radial-1
-      ny = self.n_n
-      self.kperp = np.sqrt(np.outer(self.kx[:]**2,np.ones(ny))+
-                           np.outer(np.ones(nx),self.ky[:]**2))
+      if not self.zf_test:
+         nx = self.n_radial-1
+         ny = self.n_n
+         self.kperp = np.sqrt(np.outer(self.kx[:]**2,np.ones(ny))+
+                              np.outer(np.ones(nx),self.ky[:]**2))
 
       #-----------------------------------------------------------------
 

f2py/pygacode/cgyro/data_plot.py

@@ -353,6 +353,10 @@ def plot_phi(self,xin):
       moment  = xin['moment']
       spec    = xin['spec']
 
+      if self.n_n < 2:
+         print('ERROR: (plot_phi) This plot not suitable for a single toroidal mode.')
+         return [None,]*4
+
       t = self.getnorm(xin['norm'])  
 
       if xin['fig'] is None:
@@ -944,56 +948,55 @@ def plot_zf(self,xin):
       if xin['fig'] is None:
          fig = plt.figure(MYDIR,figsize=(xin['lx'],xin['ly']))
 
-      if self.n_n > 1:
-         print('ERROR: (plot_zf) This plot option valid for ZF test only.')
-         return
+      ax = fig.add_subplot(111)
+      ax.grid(which="both",ls=":")
+      ax.grid(which="major",ls=":")
+      ax.set_xlabel(TIME)
+      ax.set_ylabel(r'$\mathrm{Re}\left( \delta\phi/\delta\phi_0 \right)$')
 
-      t  = self.t
-      k0 = 2*np.pi/self.length
+      t = self.getnorm(xin['norm'])
+      self.getbigfield()
 
-      print('INFO: (plot_zf) Using index theta index n_theta/3+1')
+      print('INFO: (plot_zf) Using theta index n_theta/3+1')
+      nselect=0 
       if field == 0:
-         f = self.phib[0,self.n_theta//3,:]
+         f = self.kxky_phi[0,:,nselect,0,:]
       elif field == 1:
-         f = self.aparb[0,self.n_theta//3,:]
+         f = self.kxky_apar[0,:,nselect,0,:]
       else:
-         f = self.bparb[0,self.n_theta/3,:]
+         f = self.kxky_bpar[0,:,nselect,0,:]
 
-      # Initialization in CGYRO is with 1e-6*besselj0 # phic[0]
-      gfactor = 1e6*(1-np.i0(k0**2)*np.exp(-k0**2))/(np.i0(k0**2)*np.exp(-k0**2))
+      for i,k0 in enumerate(self.kx):
+         # Initialization in CGYRO is with 1e-6*besselj0 # phic[0]
+         gfactor = 1e6*(1-np.i0(k0**2)*np.exp(-k0**2))/(np.i0(k0**2)*np.exp(-k0**2))
 
-      y = f*gfactor
-
-      #----------------------------------------------------
-      # Average calculations
-      imin,imax=time_index(t,w)
-      ave  = time_average(y,t,imin,imax)
-      print('INFO: (plot_zf) Integral time-average = {:.4f}'.format(ave))
+         y = f[i,:]*gfactor
+         ax.plot(t,y,label=self.kxstr+'$={:.4f}$'.format(k0))
 
-      ave_vec = ave*np.ones(len(t))
-      #----------------------------------------------------
+         #----------------------------------------------------
+         # Average calculations
+         imin,imax=time_index(t,w)
+         ave = time_average(y,t,imin,imax)
+         print('INFO: (plot_zf) Time average = {:.4f} for kx = {:.4f}'.format(ave,k0))
+         ave_vec = ave*np.ones(len(t))
+         #----------------------------------------------------
 
-      ax = fig.add_subplot(111)
-      ax.grid(which="both",ls=":")
-      ax.grid(which="major",ls=":")
-      ax.set_xlabel(TIME)
-      ax.set_ylabel(r'$\mathrm{Re}\left( \delta\phi/\delta\phi_0 \right)$')
+      if len(self.p) == 1:
 
-      ax.plot(t,y,label=r'$k_x=%.3g$' % k0)
-
-      ax.plot(t[imin:],ave_vec[imin:],color='b',
-              label=r'$\mathrm{Average}$',linewidth=1)
+         ax.plot(t[imin:],ave_vec[imin:],color='b',
+                 label=r'$\mathrm{average}$',linewidth=1)
 
-      theory = 1.0/(1.0+1.6*self.q**2/np.sqrt(self.rmin/self.rmaj))
-      ax.plot([0,max(t)],[theory,theory],color='grey',
-              label=r'$\mathrm{RH \; theory}$',alpha=0.3,linewidth=4)
+         theory = 1.0/(1.0+1.6*self.q**2/np.sqrt(self.rmin/self.rmaj))
+         ax.plot([0,max(t)],[theory,theory],color='grey',
+                 label=r'$\mathrm{RH \; theory}$',alpha=0.3,linewidth=4)
 
-      theory2 = 1./(1.0+2.*self.q**2)
-      ax.plot([0,max(t)],[theory2,theory2],color='m',
-              label=r'$\mathrm{fluid \; theory}$',alpha=0.3,linewidth=4)
+         theory2 = 1./(1.0+2.*self.q**2)
+         ax.plot([0,max(t)],[theory2,theory2],color='m',
+                 label=r'$\mathrm{fluid \; theory}$',alpha=0.3,linewidth=4)
 
       ax.legend(loc=1,prop={'size':14})
-
+      ax.set_xlim(0,t[-1])
+
       fig.tight_layout(pad=0.3)
 
       return

f2py/pygacode/cgyro/data_plot_single.py

@@ -17,6 +17,12 @@
 
 plot_type = sys.argv[3]
 
+# Rule out most plots for ZF test
+if data_in.zf_test:
+   zf_test = ['zf','geo','error']
+   if plot_type not in zf_test:
+      plot_type = 'zf_not'
+
 xin = data_in.plot_dictinit()
 
 xin['fig']    = None
@@ -150,6 +156,10 @@
 
    head = data_in.plot_hball(xin)
 
+elif plot_type == 'zf_not':
+
+   print('ERROR: (data_plot_single) Not available with ZF test.')
+
 else:
 
    print('ERROR: (data_plot_single) Plot type not found')