Add facility for printing the averaged fluxes to screen and give the …

…user the ability to specify at command line the minimum t (in vref/lref) for the start of the time averaging.

xarray_post_processing/plot_stella.py

@@ -16,16 +16,20 @@
 from utils import plot_fields_nc
 from utils import plot_fields_fluxes_spectra_nc
 from utils import plot_fields_kyspectra_with_time_nc
+from utils import print_average_fluxes
 
 # define the command line inputs
 parser = argparse.ArgumentParser()
 ncfilehelpstr = "the file_name of the stella netCDF4 file_name.out.nc file to be analysed"
 parser.add_argument("ncfile", help=ncfilehelpstr)
+parser.add_argument("-t","--tmin", help="minimum value for time averaging window", default=-1.0, type=float)
 # read the command line inputs
 args = parser.parse_args()
+tmin = args.tmin
+
+
 # get the current working directory to construct an absolute path
 workdir = os.path.abspath(os.getcwd())
-
 filename = workdir + "/" + args.ncfile
 
 stelladata = xr.open_dataset(filename+".out.nc")
@@ -36,6 +40,8 @@
 charge, mass, dens, temp, tprim, fprim, vnew, typeint, typestring = species_data(stelladata,filename)
 stelladata.close()
 
+# print information about the average fluxes
+print_average_fluxes(filename, pflx, vflx, qflx, time, typestring, tmin=tmin)
 # plot the heat, particle and momentum fluxes
 plot_fluxes_nc(filename, pflx, vflx, qflx, time, typestring)
 # plot the fields

xarray_post_processing/utils.py

@@ -187,15 +187,50 @@ def plot_1d_semilog_list_pdf (xlist,ylist,marker_list,xlab, pdf,
     plt.close (fig)
     return
 
-def timeavg(ft,time,axis=-1):
+def timeavg(ft,time,axis=-1,tmin=-1.0):
     ntime = time.size
-    it_min = ntime//2 + 1
+    zero = 1.0e-8
+    if (not tmin < -zero):
+        if tmin < time[0]:
+            it_min = 0
+        elif abs(tmin - time[0]) < zero:
+            it_min = 0
+        elif tmin > time[ntime-1]:
+            it_min = ntime - 1
+        else: #find it given tmin
+            for it in range(0,ntime-1):
+                if (tmin - time[it])*(time[it+1]-tmin) > 0.0 or abs(tmin - time[it]) < zero:
+                    it_min = it
+                    break
+    else:
+        # use half of time window
+        it_min = ntime//2 + 1
     it_max = ntime - 1
     time_interval = time[it_max] - time[it_min]
     favg = simps(ft[it_min:it_max+1],x=time[it_min:it_max+1], axis=axis) / time_interval
     return favg
 
-
+def print_average_fluxes(filename, pflx, vflx, qflx, time, typestring, tmin=-1.0):
+    nspec = len(typestring)
+    print("")
+    print("time averaged fluxes")
+    print("####################")
+    print("pflx/pflx_GBref ")
+    pflx_average = timeavg(pflx,time,axis=0,tmin=tmin)
+    for ispec in range(0,nspec):
+        print("pflx "+typestring[ispec]," = ",pflx_average[ispec])
+    print("####################")
+    print("vflx/vflx_GBref ")
+    vflx_average = timeavg(vflx,time,axis=0,tmin=tmin)
+    for ispec in range(0,nspec):
+        print("vflx "+typestring[ispec]," = ",vflx_average[ispec])
+    print("####################")
+    print("qflx/qflx_GBref ")
+    qflx_average = timeavg(qflx,time,axis=0,tmin=tmin)
+    for ispec in range(0,nspec):
+        print("qflx "+typestring[ispec]," = ",qflx_average[ispec])
+    print("####################\n")
+    return pflx_average      
 
 def plot_fluxes_nc(filename, pflx, vflx, qflx, time, typestring):
     nspec = len(typestring)