From f1be00116a3cd226e8f6ce0c63b32c9292cb7edf Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Thu, 18 Apr 2024 22:11:20 -0700 Subject: [PATCH] Proper implementation of multiple N_RADIAL capability for ZF test --- cgyro/src/cgyro_check.f90 | 2 +- cgyro/src/cgyro_equilibrium.F90 | 6 +- cgyro/src/cgyro_init_h.f90 | 11 ++-- cgyro/src/cgyro_init_kernel.F90 | 4 +- cgyro/src/cgyro_make_profiles.F90 | 14 +++-- cgyro/src/cgyro_write_initdata.f90 | 27 ++++----- f2py/pygacode/cgyro/data.py | 20 +++++-- f2py/pygacode/cgyro/data_plot.py | 75 +++++++++++++------------ f2py/pygacode/cgyro/data_plot_single.py | 10 ++++ 9 files changed, 91 insertions(+), 78 deletions(-) diff --git a/cgyro/src/cgyro_check.f90 b/cgyro/src/cgyro_check.f90 index b9003e345..dde7fad18 100644 --- a/cgyro/src/cgyro_check.f90 +++ b/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 diff --git a/cgyro/src/cgyro_equilibrium.F90 b/cgyro/src/cgyro_equilibrium.F90 index 5a47f42bd..453ea1681 100644 --- a/cgyro/src/cgyro_equilibrium.F90 +++ b/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 diff --git a/cgyro/src/cgyro_init_h.f90 b/cgyro/src/cgyro_init_h.f90 index 643a3e29d..7df95651b 100644 --- a/cgyro/src/cgyro_init_h.f90 +++ b/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 diff --git a/cgyro/src/cgyro_init_kernel.F90 b/cgyro/src/cgyro_init_kernel.F90 index 34272a52e..10d4e1033 100644 --- a/cgyro/src/cgyro_init_kernel.F90 +++ b/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) diff --git a/cgyro/src/cgyro_make_profiles.F90 b/cgyro/src/cgyro_make_profiles.F90 index 383c25ea9..843e75c36 100644 --- a/cgyro/src/cgyro_make_profiles.F90 +++ b/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 diff --git a/cgyro/src/cgyro_write_initdata.f90 b/cgyro/src/cgyro_write_initdata.f90 index 6c080bc5d..ded309459 100644 --- a/cgyro/src/cgyro_write_initdata.f90 +++ b/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 diff --git a/f2py/pygacode/cgyro/data.py b/f2py/pygacode/cgyro/data.py index 51e083de5..c55127e40 100644 --- a/f2py/pygacode/cgyro/data.py +++ b/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)) #----------------------------------------------------------------- diff --git a/f2py/pygacode/cgyro/data_plot.py b/f2py/pygacode/cgyro/data_plot.py index 248fd3475..2c649e58d 100644 --- a/f2py/pygacode/cgyro/data_plot.py +++ b/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 diff --git a/f2py/pygacode/cgyro/data_plot_single.py b/f2py/pygacode/cgyro/data_plot_single.py index 50442fb90..42406e7a4 100644 --- a/f2py/pygacode/cgyro/data_plot_single.py +++ b/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')