Workshop

Energy_spectrum/spectrum.py

@@ -6,9 +6,9 @@
 import matplotlib.pyplot as plt
 
 #velocities = ["u", "v", "w"]
-#velocities = ["u", "v", "w", "cloud_rw_mom3", "rv", "th", "RH", "aerosol_rw_mom3"]
+velocities = ["u", "v", "w", "cloud_rw_mom3", "rv", "th", "RH", "aerosol_rw_mom3"]
 #velocities = ["cloud_rw_mom3"]
-velocities = ["w"]
+# velocities = ["w"]
 
 time_start = int(argv[1])
 time_end = int(argv[2])
@@ -18,54 +18,54 @@
 
 directories = argv[6:len(argv):2]
 labels = argv[7:len(argv):2]
-print directories, labels
+print(directories, labels)
 
 # read in nx, ny, nz
 for directory, lab in zip(directories, labels):
-  w3d = h5py.File(directory + "/timestep" + str(time_start).zfill(10) + ".h5", "r")["u"][:,:,:]
+  w3d = h5py.File(directory + "/timestep" + str(time_start).zfill(10) + ".h5", "r")["u"][:,:]
   nx, ny, nz = w3d.shape
   Exy_avg = {}
   for vel in velocities:
     Exy_avg[vel] = np.zeros(((nx+1)/2))
-  
+
   for t in range(time_start, time_end+1, outfreq):
     filename = directory + "/timestep" + str(t).zfill(10) + ".h5"
-    print filename
-  
+    print(filename)
+
     for vel in velocities:
-  
-      w3d = h5py.File(filename, "r")[vel][:,:,:] # * 4. / 3. * 3.1416 * 1e3 
-      
+
+      w3d = h5py.File(filename, "r")[vel][:,:,:] # * 4. / 3. * 3.1416 * 1e3
+
       for lvl in range(from_lvl, to_lvl+1):
         w2d = w3d[:, :, lvl]
-        
+
         wkx = 1.0 / np.sqrt(nx - 1) * np.fft.rfft(w2d, axis = 0)
         wky = 1.0 / np.sqrt(ny - 1) * np.fft.rfft(w2d, axis = 1)
-        
+
         #Ex = 0.5 * (np.abs(wkx) ** 2)
         Ex = (np.abs(wkx) ** 2)
         Ex = np.mean(Ex, axis = 1)
         #Ey = 0.5 * (np.abs(wky) ** 2)
         Ey = (np.abs(wky) ** 2)
         Ey = np.mean(Ey, axis = 0)
-        
+
         Exy = 0.5 * (Ex + Ey)
         Exy_avg[vel] += Exy
-  
+
       K = np.fft.rfftfreq(nx - 1)
   #    plt.loglog(K, Exy)
       lmbd = 50. / K # assume dx=50m
-    
+
     if (t == time_start and lab==labels[0]):
       plt.loglog(lmbd, 2e-1* K**(-5./3.) )
-  
+
   for vel in velocities:
     Exy_avg[vel] /= (time_end - time_start) / outfreq + 1
     Exy_avg[vel] /= to_lvl+1 - from_lvl
   #crudely scale
   #  Exy_avg[vel] /= Exy_avg[vel][len(Exy_avg[vel])-1]
     plt.loglog(lmbd, Exy_avg[vel] , linewidth=2, label=lab+"_"+vel)
- 
+
 plt.xlim(10**4,10**2)
 plt.xlabel("l[m]")
 plt.ylabel("PSD")

Energy_spectrum/spectrum_2D.py

@@ -0,0 +1,79 @@
+# (C) Maciej Waruszewski
+
+import h5py
+import numpy as np
+from sys import argv
+import matplotlib.pyplot as plt
+
+# velocities = ["u", "v", "w"]
+# velocities = ["u", "v", "w", "cloud_rw_mom3", "rv", "th", "RH", "aerosol_rw_mom3"]
+velocities = ["u", "w", "cloud_rw_mom3", "rv", "th", "RH", "actrw_rw_mom3"]
+# velocities = ["cloud_rw_mom3", "actrw_rw_mom3"]
+# velocities = ["u", "w"]
+
+time_start = int(argv[1])
+time_end = int(argv[2])
+outfreq = int(argv[3])
+from_lvl = int(argv[4])
+to_lvl = int(argv[5])
+
+directories = argv[6:len(argv):2]
+labels = argv[7:len(argv):2]
+print(directories, labels)
+
+# read in nx, ny, nz
+for directory, lab in zip(directories, labels):
+  w3d = h5py.File(directory + "/timestep" + str(time_start).zfill(10) + ".h5", "r")["u"][:,:]
+  nx, nz = w3d.shape
+  Exy_avg = {}
+  for vel in velocities:
+    Exy_avg[vel] = np.zeros(int((nx+1)/2))
+
+  for t in range(time_start, time_end+1, outfreq):
+    filename = directory + "/timestep" + str(t).zfill(10) + ".h5"
+    print(filename)
+
+    for vel in velocities:
+      if (vel == "cloud_rw_mom3") | (vel == "actrw_rw_mom3"):
+        w3d = h5py.File(filename, "r")[vel][:,:] * 4. / 3. * 3.1416 * 1e3
+      else:
+        w3d = h5py.File(filename, "r")[vel][:,:] # * 4. / 3. * 3.1416 * 1e3
+
+      for lvl in range(from_lvl, to_lvl+1):
+        w2d = w3d[:, lvl]
+        wkx = 1.0 / np.sqrt(nx - 1) * np.fft.rfft(w2d, axis = 0)
+
+        # wky = 1.0 / np.sqrt(nz - 1) * np.fft.rfft(w2d, axis = 1)
+
+        #Ex = 0.5 * (np.abs(wkx) ** 2)
+        Ex = (np.abs(wkx) ** 2)
+        # Ex = np.mean(Ex)
+        #Ey = 0.5 * (np.abs(wky) ** 2)
+        # Ey = (np.abs(wky) ** 2)
+        # Ey = np.mean(Ey, axis = 0)
+
+        # Exy = 0.5 * (Ex + Ey)
+        Exy = Ex
+        Exy_avg[vel] += Exy
+
+      K = np.fft.rfftfreq(nx - 1)
+  #    plt.loglog(K, Exy)
+      lmbd = 100. / K # assume dx=50m
+
+    if (t == time_start and lab==labels[0]):
+      plt.loglog(lmbd, 2e-1* K**(-5./3.), label="-5/3" )
+
+  for vel in velocities:
+    Exy_avg[vel] /= (time_end - time_start) / outfreq + 1
+    Exy_avg[vel] /= to_lvl+1 - from_lvl
+  #crudely scale
+  #  Exy_avg[vel] /= Exy_avg[vel][len(Exy_avg[vel])-1]
+    plt.loglog(lmbd, Exy_avg[vel] , linewidth=2, label=lab+"_"+vel)
+
+plt.xlim(2*10**4,10**2)
+plt.xlabel("l[m]")
+plt.ylabel("PSD")
+# plt.legend()
+plt.grid(True, which='both', linestyle='--')
+# plt.title("Mean PSD of w 322m<z<642m @3h")
+plt.show()

Energy_spectrum/spectrum_2D_multi_many.py

@@ -0,0 +1,101 @@
+# (C) Maciej Waruszewski
+
+import h5py
+import numpy as np
+from sys import argv
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import os
+'''
+How to run
+
+python3 spectrum_2D_multi_many.py 14400 18400 240 10 30 /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Energy/ /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Dane/SD100/ "SD100" /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Dane/SD1000/ "SD1000" /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Dane/SD10000/ "SD10000"
+'''
+# velocities = ["u", "v", "w"]
+# velocities = ["u", "v", "w", "cloud_rw_mom3", "rv", "th", "RH", "aerosol_rw_mom3"]
+# velocities = ["u", "w", "cloud_rw_mom3", "rv", "th", "RH", "actrw_rw_mom3"]
+# velocities = ["cloud_rw_mom3", "actrw_rw_mom3"]
+# velocities = ["cloud_rw_mom3"]
+velocities = ["u", "w"]
+
+time_start = int(argv[1])
+time_end = int(argv[2])
+outfreq = int(argv[3])
+from_lvl = int(argv[4])
+to_lvl = int(argv[5])
+outfile = argv[6]
+directories = argv[7:len(argv):2]
+labels = argv[8:len(argv):2]
+
+def Energy(directories, labels, velocities):
+    # read in nx, ny, nz
+    for directory, lab in zip(directories, labels):
+        # path_file = []
+        path_file = os.listdir(directory)
+        Exy_mean = {}
+        for path in path_file:
+            path_to_file = os.listdir(directory+path)
+            joined = os.path.join(directory,path)
+            #print(joined.split("/", -1)[-1])
+
+            Exy_avg_many = [0 for i in range(len(path_to_file))]
+            Exy_avg = [0 for i in range(len(velocities))]
+            # Exy_mean = [0 for i in range(len(velocities))]
+            print(joined+ "/timestep" + str(time_start).zfill(10) + ".h5", "elo")
+            for file in range(len(path_to_file)):
+              w3d = h5py.File(joined+ "/timestep" + str(time_start).zfill(10) + ".h5", "r")["u"][:,:]
+              print(joined+ "/timestep" + str(time_start).zfill(10) + ".h5")
+              # w3d = h5py.File(joined+'/'+path_to_file[file] + "/timestep" + str(time_start).zfill(10) + ".h5", "r")["u"][:,:]
+              nx, nz = w3d.shape
+              for vel in range(len(velocities)):
+                Exy_avg[vel] = np.zeros(int((nx+1)/2))
+                # Exy_avg_many[file][vel] =np.zeros(int((nx+1)/2))
+
+              for t in range(time_start, time_end+1, outfreq):
+                filename = joined+"/timestep" + str(t).zfill(10) + ".h5"
+                # filename = joined+'/'+path_to_file[file] + "/timestep" + str(t).zfill(10) + ".h5"
+                # print(filename)
+
+                for vel in range(len(velocities)):
+                  if (velocities[vel] == "cloud_rw_mom3") | (velocities[vel] == "actrw_rw_mom3"):
+                    w3d = h5py.File(filename, "r")[velocities[vel]][:,:] * 4. / 3. * 3.1416 * 1e3
+                  else:
+                    w3d = h5py.File(filename, "r")[velocities[vel]][:,:] # * 4. / 3. * 3.1416 * 1e3
+
+                  for lvl in range(from_lvl, to_lvl+1):
+                    w2d = w3d[:, lvl]
+                    wkx = 1.0 / np.sqrt(nx - 1) * np.fft.rfft(w2d, axis = 0)
+                    Ex = (np.abs(wkx) ** 2)
+                    Exy = Ex
+                    Exy_avg[vel] += Exy
+                Exy_avg_many[file] = Exy_avg[-len(velocities):]
+            Exy_mean[path] = np.mean(Exy_avg_many[file],axis=0)
+
+        K = np.fft.rfftfreq(nx - 1)
+        lmbd = 100. / K # assume dx=100m
+
+        for path in path_file:
+            Exy_mean[path] /= (time_end - time_start) / outfreq + 1
+            Exy_mean[path] /= to_lvl+1 - from_lvl
+            plt.loglog(lmbd, Exy_mean[path] , linewidth=2, label=velocities+'_'+lab+"_"+path)
+'''
+>>>>>>> 9ceb2c4e37b01d11666e449561e48044c01296e6
+    plt.loglog(lmbd, 2e-1* K**(-5./3.), label="-5/3" )
+    plt.xlim(2*10**4,10**2)
+    plt.xlabel("l[m]")
+    plt.ylabel("PSD")
+    plt.legend()
+    plt.grid(True, which='both', linestyle='--')
+        # plt.title("Mean PSD of w 322m<z<642m @3h")
+<<<<<<< HEAD
+    #print(outfile + 'Energy_spc.png')
+    plt.savefig(outfile + 'Energy_spc_vel_u.png')
+
+Energy(directories, labels)
+=======
+    plt.savefig(outfile + 'Energy_spc.png')
+    plt.show()
+'''
+Energy(directories, labels, velocities[0])
+# Energy(directories, labels, velocities[1])

Energy_spectrum/spectrum_2D_multi_many_mod_2.py

@@ -0,0 +1,84 @@
+# (C) Maciej Waruszewski
+
+import h5py
+import numpy as np
+from sys import argv
+import matplotlib.pyplot as plt
+import os
+'''
+How to run
+
+python3 spectrum_2D_multi_many.py 14400 18400 240 10 30 /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Energy/ /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Dane/SD100/ "SD100" /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Dane/SD1000/ "SD1000" /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Dane/SD10000/ "SD10000"
+
+python3 spectrum_2D_multi_many_modyfikacje.py 14400 18400 240 10 30 /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Energy/ /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Dane/SD100/VF/ "SD100_VF" /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Dane/SD100/multi/ "SD100_multi" /home/piotr-pc/Piotr/WORKSHOPS/Dane_do_AF_2D/Dane/SD100/tail/ "SD100_tail"
+
+'''
+# velocities = ["u", "v", "w"]
+# velocities = ["u", "v", "w", "cloud_rw_mom3", "rv", "th", "RH", "aerosol_rw_mom3"]
+# velocities = ["u", "w", "cloud_rw_mom3", "rv", "th", "RH", "actrw_rw_mom3"]
+# velocities = ["cloud_rw_mom3", "actrw_rw_mom3"]
+# velocities = ["cloud_rw_mom3"]
+velocities = ["u", "w"]
+
+time_start = int(argv[1])
+time_end = int(argv[2])
+outfreq = int(argv[3])
+from_lvl = int(argv[4])
+to_lvl = int(argv[5])
+outfile = argv[6]
+directories = argv[7:len(argv):2]
+labels = argv[8:len(argv):2]
+
+
+def Energy(directories, labels):
+    for directory, lab in zip(directories, labels):
+        path_file = os.listdir(directory)
+        Exy_mean = {}
+        Exy_avg_many = {}
+        Exy_avg = {}
+        for path in path_file:
+            path_to_file = os.listdir(directory+path)
+            joined = os.path.join(directory,path)
+            w3d = h5py.File(joined+"/"+path+"_out_lgrngn"+ "/timestep" + str(time_start).zfill(10) + ".h5", "r")["u"][:,:]
+            nx, nz = w3d.shape
+            for vel in velocities:
+                Exy_avg[path, vel] = np.zeros(int((nx+1)/2))
+            for t in range(time_start, time_end+1, outfreq):
+                filename = joined+"/"+path+"_out_lgrngn"+"/timestep" + str(t).zfill(10) + ".h5"
+                for vel in velocities:
+                    if (vel == "cloud_rw_mom3") | (vel == "actrw_rw_mom3"):
+                        w3d = h5py.File(filename, "r")[vel][:,:] * 4. / 3. * 3.1416 * 1e3
+                    else:
+                        w3d = h5py.File(filename, "r")[vel][:,:]
+                    for lvl in range(from_lvl, to_lvl+1):
+                        w2d = w3d[:, lvl]
+                        wkx = 1.0 / np.sqrt(nx - 1) * np.fft.rfft(w2d, axis = 0)
+                        Ex = (np.abs(wkx) ** 2)
+                        Exy = Ex
+                        Exy_avg[path, vel] += Exy
+
+        K = np.fft.rfftfreq(nx - 1)
+        lmbd = 100. / K # assume dx=100m
+        A = np.zeros((len(velocities), len(Exy_avg[path, vel]) ))
+        for vel in range(len(velocities)):
+            for path in path_file:
+                Exy_avg[path,velocities[vel]] /= (time_end - time_start) / outfreq + 1
+                Exy_avg[path,velocities[vel]] /= to_lvl+1 - from_lvl
+                # Exy_mean[directory, vel] = Exy_avg[vel]
+                A[vel,:] += Exy_avg[path,velocities[vel]]
+                # print(Exy_avg[path,velocities[vel]])
+            A = A/len(path_file)
+            plt.loglog(lmbd, A[vel] , linewidth=2, label=velocities[vel]+'_'+lab)
+
+    plt.loglog(lmbd, 2e-1* K**(-5./3.), label="-5/3" )
+    plt.xlim(2*10**4,10**2)
+    plt.xlabel("l[m]")
+    plt.ylabel("PSD")
+    plt.legend()
+    plt.grid(True, which='both', linestyle='--')
+        # plt.title("Mean PSD of w 322m<z<642m @3h")
+    plt.savefig(outfile + 'Energy_spc.png')
+    plt.show()
+
+Energy(directories, labels)
+# Energy(directories, labels, velocities[1])