From a54e5fa34cb6cbadfa3d3905b8432b0cd9fbc21f Mon Sep 17 00:00:00 2001
From: Galina Avdeeva <gavdeeva@login34.chn.perlmutter.nersc.gov>
Date: Fri, 20 Sep 2024 14:27:05 -0700
Subject: [PATCH] cgyro_converge: fix error in the plot of averaged windows and
 add option to swich between analysis of one case or all cases in the
 directory

---
 cgyro/bin/cgyro_converge | 66 ++++++++++++++++++++++++----------------
 1 file changed, 40 insertions(+), 26 deletions(-)

diff --git a/cgyro/bin/cgyro_converge b/cgyro/bin/cgyro_converge
index b83827c48..c40e5e840 100755
--- a/cgyro/bin/cgyro_converge
+++ b/cgyro/bin/cgyro_converge
@@ -8,19 +8,25 @@ import matplotlib.pyplot as plt
 from pygacode.cgyro import data
 from pygacode.gacodefuncs import *
 
+# --- to visualize results for one case plot=True------
+plot=True
 
-directory="/global/cfs/cdirs/cgyrodb/"
-
-fout = "./cases_results.txt"
-fo = open(fout, "w")
+if plot:
+    directory='./'
+    write_output=False
+else:
+    directory= "/global/cfs/cdirs/cgyrodb/"
+    write_output=True
+    fout = "./cases_results_20%.txt"
+    fo = open(fout, "w")
 
 
 for path, folders, files in os.walk(directory):
     # Open file#
     for filename in files:
       if 'out.cgyro.info' in filename:
-                  
-        fo.write(f"#-----{path}-----#")
+        if write_output:          
+           fo.write(f"#-----{path}-----#")
         # read minimal data for setup (look in current directory)
         sim = data.cgyrodata(f'{path}/',fast=True,silent=True)
 
@@ -34,8 +40,8 @@ for path, folders, files in os.walk(directory):
 
         # Simulation length (max time)
         tmax = t[-1]
-        #if len(t)!=tmax:
-          # fo.write("!! Problem with a time array !!")
+        
+    
 
         # size and number of time windows
         #dt = 50.0
@@ -44,13 +50,14 @@ for path, folders, files in os.walk(directory):
         dt=tmax/nwin
         t0 = tmax-nwin*dt
         
+        if plot:
+         #add plot for visualization
+         fig = plt.figure(figsize=(12,6))
 
-        #add plot for visualization
-        #fig = plt.figure(figsize=(12,6))
-
-        #plt.xlabel('time [$c_s$/a]')
-        #plt.ylabel('Flux')
-        colors = mpl.cm.rainbow(np.linspace(0, 1, nwin))
+         plt.xlabel('time [$c_s$/a]')
+         plt.ylabel('Flux')
+         colors = mpl.cm.rainbow(np.linspace(0, 1, nwin))
+        
         convergance = True
 
         # select field=0 (phi), moment=1 (Q), species (0), and sum over kx
@@ -58,20 +65,19 @@ for path, folders, files in os.walk(directory):
           field=0 ; moment=1 ; 
           species=species_index
           y = np.sum(sim.ky_flux[species,moment,field,:,:],axis=0)
-          #plt.plot(t,y, label=f"Q_{species_index}", linestyle=linestyles[species_index], color=colors[species_index])
+          if plot:
+           plt.plot(t,y, label=f"Q_{species_index}", linestyle=linestyles[species_index], color=colors[species_index])
           average_array=[]
 
          # averaging over all windows (note that time_average function is fast/optimized)
           for i in range(nwin):
             w=str(tmax-(i+1)*dt)+','+str(tmax-i*dt)
-            imin,imax=time_index(t,w)
-            
+            imin,imax=time_index(t,w)               
             
-            if imax>tmax:
-               fo.write("!! Problem with a time array !!")
             ave = time_average(y[:],t,imin,imax)
             average_array.append(ave)
-            #plt.hlines(y=ave, xmin=imin, xmax=imax, color=colors[i], lw=3, linestyle=linestyles[species_index])
+            if plot:
+              plt.hlines(y=ave, xmin=t[imin], xmax=t[imax], color=colors[i], lw=3, linestyle=linestyles[species_index])
             
 
 
@@ -79,14 +85,22 @@ for path, folders, files in os.walk(directory):
         
         
         # define the convergence by std between last free averaged regions as a pesentage to the total
-          if np.std(average_array[:3]/np.mean(average_array[:3])) > 0.1:
+          if np.std(average_array[:3]/np.mean(average_array[:3])) > 0.2:
             convergance = False
-
+       
         if convergance:
-            fo.write("TRUE \n")
+            if write_output:
+              fo.write("TRUE \n")
+            else:
+               plt.title('Converged')
         else:
-            fo.write("FALSE \n")
-        #plt.legend()
-        #plt.show()
+            if write_output:
+              fo.write("FALSE \n")
+            else:
+              plt.title("Didn't converge")
+  
+        if plot:
+          plt.legend()
+          plt.show()