Refined grids plotter rework

Energy_spectrum/spectrum_refined.py

@@ -34,6 +34,7 @@
   lmbd = {}
   level_start_idx = {}
   level_end_idx = {}
+  exists = {}
 
   # read some constant parameters
   with h5py.File(directory + "/const.h5", 'r') as consth5:
@@ -58,16 +59,21 @@
   # initiliaze nx,ny,nz,dx and average energy for each variable
   for var in args.vars:
     filename = directory + "/timestep" + str(time_start_idx).zfill(10) + ".h5"
-    w3d = h5py.File(filename, "r")[var][:,:,:]
-    nx[var], ny[var], nz[var] = tuple(x for x in w3d.shape)
-    dx[var] = X / (nx[var] - 1)
-    dz[var] = Z / (nz[var] - 1) 
-    Exy_avg[var] = np.zeros(int((nx[var]-1)/2 + 1))
-    ref[var] = int(dx_adve / dx[var])
-    assert(float(args.level_start / dz[var]).is_integer())
-    assert(float(args.level_end / dz[var]).is_integer())
-    level_start_idx[var] = int(args.level_start / dz[var])
-    level_end_idx[var] = int(args.level_end / dz[var]) + 1
+
+    try:
+      w3d = h5py.File(filename, "r")[var][:,:,:] 
+      nx[var], ny[var], nz[var] = tuple(x for x in w3d.shape)
+      dx[var] = X / (nx[var] - 1)
+      dz[var] = Z / (nz[var] - 1) 
+      Exy_avg[var] = np.zeros(int((nx[var]-1)/2 + 1))
+      ref[var] = int(dx_adve / dx[var])
+      assert(float(args.level_start / dz[var]).is_integer())
+      assert(float(args.level_end / dz[var]).is_integer())
+      level_start_idx[var] = int(args.level_start / dz[var])
+      level_end_idx[var] = int(args.level_end / dz[var]) + 1
+      exists[var] = True
+    except:
+      exists[var] = False
 
 
   # time loop
@@ -78,14 +84,16 @@
     # variables loop
     for var in args.vars:
       print(var)
+      if not exists[var]:
+        continue
 
       print(nx[var],dx[var][0])
       w3d = h5py.File(filename, "r")[var][0:nx[var]-1,0:ny[var]-1,:] # * 4. / 3. * 3.1416 * 1e3 
 
       for lvl in range(level_start_idx[var], level_end_idx[var]):
         w2d = w3d[:, :, lvl]
         #print w2d
-        
+
         wkx = 1.0 / np.sqrt(nx[var] - 1) * np.fft.rfft(w2d, axis = 0)
         wky = 1.0 / np.sqrt(ny[var] - 1) * np.fft.rfft(w2d, axis = 1)
 
@@ -107,15 +115,21 @@
 #      print(K, lmbd)
 
     if (t == time_start_idx and lab==args.labels[0]):
-      plt.loglog(lmbd[var], 2e-6* K**(-5./3.) )
+      L = np.array([2e2, 2e3])
+      plt.loglog(L, 2e-5 * L**(5./3.), label = "-5/3" , color="black", ls='dotted')
+      L = np.array([5e1, 3e2])
+      plt.loglog(L, 2e-8 * L**(3.), label = "-3" , color="black", ls='dashed')
 
   for var in args.vars:
+    if not exists[var]:
+      continue
+
     Exy_avg[var] /= (time_end_idx - time_start_idx) / outfreq + 1
     Exy_avg[var] /= level_end_idx[var] - level_start_idx[var]
 
     #crudely scale
     #Exy_avg[var] /= Exy_avg[var][len(Exy_avg[var])-1]
-    Exy_avg[var] /= np.sum(Exy_avg[var])
+#    Exy_avg[var] /= np.sum(Exy_avg[var])
 
     #plot
     plt.loglog(lmbd[var], Exy_avg[var] , linewidth=2, label=lab+"_"+var)

Matplotlib_common/latex_labels.py

@@ -31,16 +31,14 @@
   "surf_precip" : 'Surf. precip. [mm/day]',
   "acc_precip" : 'Acc. precip. [mm]',
   "cl_nc" : '$N_c$ [cm$^{-3}$]',
-  "cl_nc_rico" : '$N_c$ [cm$^{-3}$]',
 #  "cl_nr" : '$N_{r>25\mu m}$ [cm$^{-3}$]',
   "cl_nr" : '$N_{r}$ [cm$^{-3}$]',
-  "cl_nr_rico" : '$N_{r}$ [cm$^{-3}$]',
   "cl_gccn_conc" : '$N_{GCCN}$ [cm$^{-3}$]',
   "thl" : r'$\theta_l$ [K]',
   "00rtot" : '$q_{t}$ [g/kg]',
   "rliq" : '$q_{l}$ [g/kg]',
   "clfrac" : 'Cloud fraction',
-  "cloud_cover_rico" : 'Cloud cover',
+  "cloud_cover" : 'Cloud cover',
   "cloud_cover_dycoms" : 'Cloud cover',
   "prflux" : 'Precip. flux [W m$^{-2}$]',
   "wvar" : r'Var$\left(w\right)$ [m$^2$ s$^{-2}$]',
@@ -53,7 +51,7 @@
   "er" : 'Entrain. rate [cm s$^{-1}$]',
   "wvarmax" : 'Max. $w$ var. [m$^{2}$ s$^{-2}$]',
   "cloud_base_dycoms" : 'Cloud base [m]',
-  "min_cloud_base_rico" : 'Lowest cloud base [m]',
+  "min_cloud_base" : 'Lowest cloud base [m]',
   "inversion_height_rico" : 'Inversion height [m]',
   "gccn_rw_cl" : '$<r>$ of GCCN droplets [um]',
   "non_gccn_rw_cl" : '$<r>$ of CCN droplets [um]',
@@ -62,13 +60,18 @@
   "clb_bigrain_mean_conc" : 'conc. of (r$>$40um) @ clbase [1/cc]',
   "clb_bigrain_mean_inclt" : 'time since act. of (r$>$40um) @ clbase [s]',
   "clb_bigrain_mean_gccn_fraction" : 'frac. of (r$>$40um) on gccn @ clbase',
-  "cl_acnv25_rico" : 'acnv. rate [g m$^{-3}$ d$^{-1}$]',
-  "cl_accr25_rico" : 'accr. rate [g m$^{-3}$ d$^{-1}$]',
-  "cl_acnv25_dycoms" : 'acnv. rate [g m$^{-3}$ d$^{-1}$]',
-  "cl_accr25_dycoms" : 'accr. rate [g m$^{-3}$ d$^{-1}$]',
-  "cloud_avg_supersat" : '$<S>$ in cloud [\%]',
-  "cl_meanr" : '$<r>$ in cloud [um]',
-  "sd_conc" : '$<SD>$'
+  "cl_acnv25" : 'acnv. rate [g m$^{-3}$ d$^{-1}$]',
+  "cl_accr25" : 'accr. rate [g m$^{-3}$ d$^{-1}$]',
+  "cl_avg_supersat" : '$<S>$ in cloud [\%]',
+  "cl_avg_th" : '$<\\theta>$ in cloud [K]',
+  "cl_avg_rv" : '$<r_v>$ in cloud [1]',
+  "cl_avg_cloud_meanr" : '$<r>$ of cloud droplets in cloud [um]',
+  "cl_avg_cloud_stddevr" : '$\sigma(r)$ of cloud droplets in cloud [um]',
+  "cl_avg_act_meanr" : '$<r>$ of activated droplets in cloud [um]',
+  "cl_avg_act_stddevr" : '$\sigma(r)$ of activated droplets in cloud [um]',
+  "sd_conc" : '$<SD>$',
+  "surf_flux_latent"   : 'latent surface flux [W/m$^2$]',
+  "surf_flux_sensible" : 'sensible surface flux [W/m$^2$]',
 }
 
 
@@ -93,4 +96,10 @@
  17 : "(r)",
  18 : "(s)",
  19 : "(t)",
+ 20 : "(u)",
+ 21 : "(v)",
+ 22 : "(w)",
+ 23 : "(x)",
+ 24 : "(y)",
+ 25 : "(z)",
 }

Matplotlib_common/plot_series.py

@@ -1,17 +1,31 @@
 import numpy as np
-from sys import argv
 from math import floor
+#from sys import argv
+import argparse
+import matplotlib.pyplot as plt
 
 from latex_labels import var_labels
 from read_UWLCM_arrays import *
 
+
 def plot_series(var_list, plot_iter, nplotx, nploty, axarr, xscaledict, yscaledict, xlimdict, ylimdict, show_bin=False, suffix='', xlabel='', ylabeldict=var_labels, file_names=[], file_labels=[], linewidth=1):
+
+  parser = argparse.ArgumentParser(description='Plot UWLCM series comparison')
+  parser.add_argument("-ts", "--time_start", type=float, required=False, help="start of the plotted period [s] (override default)")
+  parser.add_argument("-te", "--time_end", type=float, required=False, help="end of the plotted period [s] (override default)")
+  parser.add_argument("-d", "--dirs", action="extend", nargs="+", type=str, help="list of directories with the data", required=True)
+  parser.add_argument("-l", "--labels", action="extend", nargs="+", type=str, help="list of labels of the data (same order as --dirs)", required=True)
+  args, extra = parser.parse_known_args()
+
+  if args.time_start is not None and args.time_end is not None:
+    xlimdict = {x: (args.time_start, args.time_end) for x in xlimdict}
+
+
   # if file names are not defined, read them and labels from command line
   if len(file_names)==0:
-    file_no = np.arange(1, len(argv)-1 , 2)
-    for no in file_no:
-      file_names.append(argv[no] + suffix)
-      file_labels.append(argv[no+1])
+    for directory, lab in zip(args.dirs, args.labels):
+      file_names.append(directory + suffix)
+      file_labels.append(lab)
 
   for var in var_list:
     label_counter=0
@@ -20,6 +34,11 @@ def plot_series(var_list, plot_iter, nplotx, nploty, axarr, xscaledict, yscaledi
       series_file = open(file_name, "r")
       my_times = read_my_var(series_file, "position")
       my_res = read_my_var(series_file, var)
+      if len(my_res) == 0: # file does not contain this type of plot
+        print("skipping from: " + str(label_counter))
+        label_counter+=1
+        print("skipping to: " + str(label_counter))
+        continue
 
       # rescale time to hours
       my_times = my_times / 3600.
@@ -32,14 +51,15 @@ def plot_series(var_list, plot_iter, nplotx, nploty, axarr, xscaledict, yscaledi
 
       linestyles = ['--', '-.', ':']
       dashList = [(3,1),(1,1),(4,1,1,1),(4,2)] 
-      colorList = ['red', 'blue', 'green']
+  #    colorList = ['red', 'blue', 'green']
+      colorList = plt.rcParams['axes.prop_cycle'].by_key()['color'] # default prop cycle colors
 
       # x label only on he lowest row
       xlabel_used = xlabel
       if plot_iter < nploty:
         xlabel_used = ''
 
-      plot_my_array(axarr, plot_iter, my_times, my_res, nploty, xlabel=xlabel_used, ylabel=ylabeldict[var], varlabel=file_labels[label_counter], dashes = dashList[label_counter % len(dashList)], xscale=xscaledict[var], yscale=yscaledict[var], xlim=xlimdict[var], ylim=ylimdict[var], linewidth=linewidth)#, color = colorList[int(floor(label_counter / len(dashList)))])
+      plot_my_array(axarr, plot_iter, my_times, my_res, nploty, xlabel=xlabel_used, ylabel=ylabeldict[var], varlabel=file_labels[label_counter], dashes = dashList[label_counter % len(dashList)], xscale=xscaledict[var], yscale=yscaledict[var], xlim=xlimdict[var], ylim=ylimdict[var], linewidth=linewidth, color = colorList[label_counter % len(colorList)])
       label_counter+=1
     plot_iter = plot_iter + 1
   return plot_iter

Matplotlib_common/read_UWLCM_arrays.py

@@ -18,6 +18,8 @@ def read_my_var(file_obj, var_name):
     arr_name = file_obj.readline() 
     if(str(arr_name).strip() == str(var_name).strip()):
       return read_my_array(file_obj)
+    elif(len(arr_name)==0):
+      return np.empty(0)
 
 #def plot_my_array(axarr, plot_iter, time, val, nploty, xlabel=None, ylabel=None, varlabel=None , linestyle='--', dashes=(5,2), xlim=None, ylim=None, xscale="linear", yscale="linear"):
 #  x = int(plot_iter / nploty)

README.md

@@ -1,2 +1,28 @@
 # UWLCM_plotting
 a set of scripts for plotting output of the UWLCM model
+
+CONTENTS
+
+plotting:
+
+- /drawbicyc - C++ program for processing UWLCM output; can calculate and plot series, profiles and snapshots for multiple parameters; also has C++ programs for averaging and comparing series or profiles from multiple runs, although this is probably better done in Python (e.g. /cases)
+
+- /cases - Python scripts for comparing series and profiles calculated using drawbicyc; specific to each LES case
+
+- /Energy_spectrum - Python scripts for plotting and comparing Fourier spectra of variables; works directly on UWLCM output
+
+- /histograms - Python scripts for plotting and comparing histograms with spatial and/or temporal distribution of variables; also capable of plotting scatter plots of correlation between two variables; works directly on UWLCM output
+
+- /NC_vs_AF - Python scripts for plotting scatter plots of number concentration of cloud droplets vs adiabatic fraction, or vs liquid water content; works directly on UWLCM output; NOTE: redundant to scatter plots in /histograms?
+
+- /papers - Python scripts for making figures used in our papers
+
+- /Size_spectra - Python scripts for comparing droplet size distributions; works directly on UWLCM output 
+
+
+helpers:
+
+- /Matplotlib_common - Python scripts for reading drawbicyc output; also plotting comparison of series and profiles from drawbicyc
+
+- /UWLCM_plotters - C++ classes for reading UWLCM output; used in drawbicyc
+

UWLCM_plotters/include/Plotter2d.hpp

@@ -1,13 +1,13 @@
 #pragma once
 #include "common.hpp"
-#include "PlotterCommon.hpp"
+#include "PlotterH5.hpp"
 
 template<int NDims>
-class Plotter_t : public PlotterCommon {};
+class Plotter_t : public PlotterH5 {};
 
 // 2d version
 template<>
-class Plotter_t<2> : public PlotterCommon 
+class Plotter_t<2> : public PlotterH5 
 {
   public:
   static const int n_dims = 2;
@@ -18,10 +18,10 @@ class Plotter_t<2> : public PlotterCommon
   arr_t dv;
 
   protected:
-  using parent_t = PlotterCommon;
+  using parent_t = PlotterH5;
   hsize_t n[2];
   enum {x, z};
-  arr_t tmp, tmp_srfc;
+  arr_t tmp, tmp_srfc, tmp_ref;
 
   public:
 
@@ -202,6 +202,7 @@ class Plotter_t<2> : public PlotterCommon
     // other dataset are of the size x*z, resize tmp
     tmp.resize(n[0]-1, n[1]-1);
     tmp_srfc.resize(n[0]-1, 1);
+    tmp_ref.resize(tmp.shape());
   }
 };
 

UWLCM_plotters/include/Plotter3d.hpp

@@ -4,7 +4,7 @@
 
 // 3d version
 template<>
-class Plotter_t<3> : public PlotterCommon 
+class Plotter_t<3> : public PlotterH5 
 {
   public:
   static const int n_dims = 3;
@@ -15,10 +15,10 @@ class Plotter_t<3> : public PlotterCommon
   arr_t dv;
 
   protected:
-  using parent_t = PlotterCommon;
+  using parent_t = PlotterH5;
   hsize_t n[3];
   enum {x, y, z};
-  arr_t tmp, tmp_srfc;
+  arr_t tmp, tmp_srfc, tmp_ref;
   blitz::Range yrange;
 
   public:
@@ -36,15 +36,17 @@ class Plotter_t<3> : public PlotterCommon
       cnt[3] = { n[x],  n[y],  srfc ? 1 : n[z] }, 
       off[3] = { 0,     0,     0    };
     this->h5s.selectHyperslab( H5S_SELECT_SET, cnt, off);
+
+    arr_t &arr(tmp.extent(0) == n[x] ? tmp : tmp_ref); // crude check if it is refined or normal data; assume surface data is never refined
 
     hsize_t ext[3] = {
-      hsize_t(tmp.extent(0)), 
-      hsize_t(tmp.extent(1)), 
-      hsize_t(srfc ? tmp_srfc.extent(2) : tmp.extent(2)) 
+      hsize_t(arr.extent(0)), 
+      hsize_t(arr.extent(1)), 
+      hsize_t(srfc ? tmp_srfc.extent(2) : arr.extent(2)) 
     };
-    this->h5d.read(srfc ? tmp_srfc.data() : tmp.data(), H5::PredType::NATIVE_FLOAT, H5::DataSpace(3, ext), h5s);
+    this->h5d.read(srfc ? tmp_srfc.data() : arr.data(), H5::PredType::NATIVE_FLOAT, H5::DataSpace(3, ext), h5s);
 
-    return blitz::safeToReturn((srfc ? tmp_srfc : tmp) + 0);
+    return blitz::safeToReturn((srfc ? tmp_srfc : arr) + 0);
   }
 
   auto h5load_timestep(
@@ -235,6 +237,44 @@ class Plotter_t<3> : public PlotterCommon
     // other dataset are of the size x*z, resize tmp
     tmp.resize(n[0]-1, n[1]-1, n[2]-1);
     tmp_srfc.resize(n[0]-1, n[1]-1, 1);
+
+    // init refined data
+    try
+    {
+      this->h5f.openDataSet("X refined").getSpace().getSimpleExtentDims(n, NULL); 
+      this->map["refined x"] = n[0]-1;
+      this->map["refined y"] = n[1]-1;
+      this->map["refined z"] = n[2]-1;
+      tmp_ref.resize(n[0], n[1], n[2]);
+      h5load(file + "/const.h5", "X refined");
+      this->map["refined dx"] = tmp_ref(1,0,0) - tmp_ref(0,0,0);
+      h5load(file + "/const.h5", "Y refined");
+      this->map["refined dy"] = tmp_ref(0,1,0) - tmp_ref(0,0,0);
+      h5load(file + "/const.h5", "Z refined");
+      this->map["refined dz"] = tmp_ref(0,0,1) - tmp_ref(0,0,0);
+      this->CellVol_ref = this->map["refined dx"] * this->map["refined dy"] * this->map["refined dz"];
+      tmp_ref.resize(n[0]-1, n[1]-1, n[2]-1);
+    }
+    catch(...) // for pre-refinement simulations that didnt store refined stuff
+    {
+      this->map["refined x"] = this->map["x"]; 
+      this->map["refined y"] = this->map["y"]; 
+      this->map["refined z"] = this->map["z"]; 
+      this->map["refined dx"] = this->map["dx"]; 
+      this->map["refined dy"] = this->map["dy"]; 
+      this->map["refined dz"] = this->map["dz"]; 
+      this->CellVol_ref = this->CellVol;
+      tmp_ref.resize(tmp.shape());
+    }
+
+    for (auto const& x : this->map)
+{
+    std::cout << x.first  // string (key)
+              << ':'
+              << x.second // string's value
+              << std::endl;
+}
+
   }
 };