new diag

Author: Georgia Acton

POST_PROCESSING/post_processing/gzvs_movie.py

@@ -1,16 +1,19 @@
 import numpy as np
-from stella_plots import movie_2d
+from stella_plots import movie_2d, plot_2d
 from stella_data import gzvs, vpa, zed, ntime, g2nozonal_vs_zvpas, nzed, time
-from stella_data import gzvs_zonal
+from stella_data import gzvs_zonal, kx_stella
 from stella_data import outdir
+import matplotlib.pyplot as plt
 
 nvpa = vpa.size
 print('nvpa', nvpa)
 gmax = np.arange(ntime,dtype=float)
 gmin = np.arange(ntime,dtype=float)
 
+kx_val = 0.05
+ikx = next((i for i, value in enumerate(kx_stella) if value >=kx_val), None)
 #(t, species, vpa, tube, zed, kx, ri)
-g_zonal = np.array(gzvs_zonal[:, 0, :, 0, :, 1, 0] + 1j*gzvs_zonal[:, 0, :, 0, :, 1, 1], dtype=complex) 
+g_zonal = np.array(gzvs_zonal[:, 0, :, 0, :, ikx, 0] + 1j*gzvs_zonal[:, 0, :, 0, :, ikx, 1], dtype=complex) 
 
 for i in range(ntime):
     gmax[i] = np.abs(g_zonal[i,:,:]).max()
@@ -19,7 +22,16 @@
 
 ylabel = '$v_{\parallel}$'
 xlabel = '$z$'
-title = '$\int d\mu \int d^2 \mathbf{R} g^2$'
+title = '$\int d\mu \int d^2 \mathbf{R} g$'
 movie_file = outdir + 'gzvs.mp4'
 
-movie_2d(g_zonal,zed,vpa,gmin,gmax,ntime-1,movie_file,xlabel,ylabel,title,cmp='RdBu', abs_flag=True)
+#movie_2d(g_zonal,zed,vpa,gmin,gmax,ntime-1,movie_file,xlabel,ylabel,title,cmp='RdBu', abs_flag=False)
+
+time_avg = ntime//2
+g_zonal_avg = np.sum(g_zonal[time_avg:,:,:],axis=0)
+gmax = np.abs(g_zonal_avg).max()
+gmin = - gmax
+
+outfile = outdir + 'timeave.pdf'
+plot_2d(g_zonal_avg,zed,vpa,gmin,gmax,'zed','vpa',title,cmp='RdBu', outfile=outfile)
+

POST_PROCESSING/post_processing/kspectra_plots.py

@@ -48,49 +48,7 @@
 #f2out.write('#ky0: '+str(ky0)+' ky1: '+str(ky1)+'\n')
 f2out.close()
 
-#fig=plt.figure(figsize=(12,8))
-fig, axs = plt.subplots(ncols=2, layout='constrained', sharex='col', figsize=(12,8))
-
-# distinguish between zonal
-for ikx in range(nakx):
-    axs[0].semilogy(time,phi2[:,ikx,0])
-    print('phi2(t=0), kx=', np.sqrt(phi2[0,ikx,0])/kx[ikx]**2 )
-    axs[1].scatter(kx[ikx], np.sqrt(phi2[-1, ikx,0]/phi2[0,ikx,0]))
-    print('ratio =', np.sqrt(phi2[-1, ikx,0]/phi2[0,ikx,0]), 'kx=',kx[ikx])
-    
-axs[1].set_xscale('log')
-axs[0].set_xlabel('$t (v_{t}/a)$')
-axs[1].set_label('$k_x \rho_i$') 
-axs[0].set_xlim([time[0],time[ntime-1]])
-axs[1].set_xlim(0,2.1 )
-plt.title('$\Phi(k_x,k_y=0)$')
-
-file = outdir+file_prefix+'.stella_kspectra_zonal.pdf'
-pdf = PdfPages(file)
-for i in plt.get_fignums():
-    pdf.savefig(i)
-pdf.close()
-plt.close('all')
-
-
 fig=plt.figure(figsize=(12,8))
-
-for ikx in range(1, nakx):
-    plt.semilogy(time, entropy_vs_kx[:,ikx])
-
-plt.semilogy(time, np.sum(entropy_vs_kx[:,:], axis=1), '--') 
-plt.xlabel('$t (v_{t}/a)$')
-plt.ylabel('Free energy')
-file = outdir+file_prefix+'.stella_free_energy.pdf'
-pdf = PdfPages(file)
-for i in plt.get_fignums():
-    pdf.savefig(i)
-pdf.close()
-plt.close('all')
-
-
-fig=plt.figure(figsize=(12,8))
-
 # distinguish between zonal/nonzonal modes
 for ikx in range(phi2.shape[1]-1):
     plt.semilogy(time,phi2[:,ikx+1,0],'--')

POST_PROCESSING/post_processing/stella_data.py

@@ -8,7 +8,8 @@
 
 ####### Import variables from netcdf file #########
 #infile = input("Path to netcdf file: ")
-input_directory = '/Users/giorgiaacton/Documents/runs/ALL_RH/test_diag/'
+#input_directory = '/pitagora_work/FUPA1_STELTURB/gacton00/Paul/linear_60_higherres/'
+input_directory = '/pitagora_scratch/userexternal/gacton00/RH_noj0_higherres/'
 file_prefix = 'input'
 infile = input_directory + file_prefix + '.out.nc'
 outdir = input_directory
@@ -193,8 +194,8 @@ def read_stella_float(var):
 temperature, temperature_present = \
     read_stella_float('temperature')
 
-entropy_vs_kx, entropy_vs_kx_present = \
-    read_stella_float('entropy_vs_kx')
+free_energy_vs_kx, free_energy_vs_kx_present = \
+    read_stella_float('free_energy_vs_kx')
 
 ncfile.close()
 

POST_PROCESSING/post_processing/stella_plots.py

@@ -51,7 +51,9 @@ def logyplot_1d(x,y,xlab,title='',ylab=''):
         plt.title(title)
     return fig
 
-def plot_2d(z,xin,yin,zmin,zmax,xlab='',ylab='',title='',cmp='RdBu'):
+def plot_2d(z,xin,yin,zmin,zmax,xlab='',ylab='',title='',cmp='RdBu', outfile=''):
+
+    z = np.real(z)
 
     fig = plt.figure(figsize=(12,8))
     x,y = np.meshgrid(xin,yin)
@@ -63,8 +65,9 @@ def plot_2d(z,xin,yin,zmin,zmax,xlab='',ylab='',title='',cmp='RdBu'):
     plt.xlabel(xlab)
     plt.ylabel(ylab)
     plt.title(title)
-    return fig
-
+    plt.savefig(outfile)
+    plt.close(fig)
+    
 def movie_2d(z, xin, yin, zmin, zmax, nframes, outfile,
              xlab='', ylab='', title='', step=1, cmp='RdBu',
              fps=10, interval=50, abs_flag=True):
@@ -80,8 +83,11 @@ def movie_2d(z, xin, yin, zmin, zmax, nframes, outfile,
 
     if abs_flag == True :
         print('abs flag is True')
-        z_even = np.sign(np.real(z_even_full)) * np.abs(z_even_full)
-        z_odd  = np.sign(np.real(z_odd_full)) * np.abs(z_odd_full)
+        #z_even = np.sign(np.real(z_even_full)) * np.abs(z_even_full)
+        #z_odd  = np.sign(np.real(z_odd_full)) * np.abs(z_odd_full)
+        z_even = np.abs(z_even_full)
+        z_odd  = np.abs(z_odd_full)
+        zmin[:] = 0.0
     else:
         z_even = np.real(z_even_full)
         z_odd  = np.real(z_odd_full) 
@@ -104,9 +110,9 @@ def movie_2d(z, xin, yin, zmin, zmax, nframes, outfile,
     # Create initial images once (frame 0)
     i0 = frames[0]
     im1 = ax1.imshow(z_even[i0], origin='lower', aspect='auto',
-                     extent=extent, cmap='RdBu', animated=False)
+                     extent=extent, cmap=cmp, animated=False)
     im2 = ax2.imshow(z_odd[i0], origin='lower', aspect='auto',
-                     extent=extent, cmap='RdBu', animated=False)
+                     extent=extent, cmap=cmp, animated=False)
 
     # Set initial clim from provided arrays/scalars:
     try: