From d1f1961427a48b59cc6cb51f84166ffd3cbe3547 Mon Sep 17 00:00:00 2001
From: Galina Avdeeva <gavdeeva@login34.chn.perlmutter.nersc.gov>
Date: Fri, 20 Sep 2024 13:01:24 -0700
Subject: [PATCH] modify script cgyro_converge

---
 cgyro/bin/cgyro_converge | 100 ++++++++++++++++++++++++++++++---------
 1 file changed, 78 insertions(+), 22 deletions(-)

diff --git a/cgyro/bin/cgyro_converge b/cgyro/bin/cgyro_converge
index 66fd27478..b83827c48 100755
--- a/cgyro/bin/cgyro_converge
+++ b/cgyro/bin/cgyro_converge
@@ -1,36 +1,92 @@
 #!/usr/bin/env python
 
 import numpy as np
+import os
+import matplotlib as mpl
+import matplotlib.pyplot as plt
 
 from pygacode.cgyro import data
 from pygacode.gacodefuncs import *
 
-# read minimal data for setup (look in current directory)
-sim = data.cgyrodata('./',fast=True,silent=True)
 
-# copy time vector and number of species
-t = sim.t
-ns = sim.n_species
+directory="/global/cfs/cdirs/cgyrodb/"
 
-# read bin.cgyro.ky_flux
-sim.getflux()
+fout = "./cases_results.txt"
+fo = open(fout, "w")
 
-# select field=0 (phi), moment=1 (Q), species (0), and sum over kx
-field=0 ; moment=1 ; species=0
-y = np.sum(sim.ky_flux[species,moment,field,:,:],axis=0)
 
-# Simulation length (max time)
-tmax = t[-1]
+for path, folders, files in os.walk(directory):
+    # Open file#
+    for filename in files:
+      if 'out.cgyro.info' in filename:
+                  
+        fo.write(f"#-----{path}-----#")
+        # read minimal data for setup (look in current directory)
+        sim = data.cgyrodata(f'{path}/',fast=True,silent=True)
 
-# size and number of time windows
-dt = 50.0
-nwin = int(tmax/dt)
-t0 = tmax-nwin*dt
+        # copy time vector and number of species
+        t = sim.t
+        ns = sim.n_species
+        linestyles=['-','--',':','-.']
+        # read bin.cgyro.ky_flux
+        sim.getflux()
 
-# 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)
-   ave = time_average(y[:],t,imin,imax)
-   print(i,ave)
+
+        # 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
+        #nwin = int(tmax/dt)
+        nwin=8
+        dt=tmax/nwin
+        t0 = tmax-nwin*dt
+        
+
+        #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))
+        convergance = True
+
+        # select field=0 (phi), moment=1 (Q), species (0), and sum over kx
+        for species_index in range(ns):
+          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])
+          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)
+            
+            
+            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])
+            
+
+
+
+        
+        
+        # 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:
+            convergance = False
+
+        if convergance:
+            fo.write("TRUE \n")
+        else:
+            fo.write("FALSE \n")
+        #plt.legend()
+        #plt.show()
+