[obsolete?] refined arrays plotting

UWLCM_plotters/include/Plotter3d.hpp

@@ -18,7 +18,7 @@ class Plotter_t<3> : public PlotterCommon
   using parent_t = PlotterCommon;
   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;
+}
+
   }
 };
 

UWLCM_plotters/include/PlotterCommon.hpp

@@ -14,7 +14,7 @@ class PlotterCommon
   std::map<std::string, double> map;
   std::map<std::string, arr_prof_t> map_prof;
   blitz::Array<float, 1> timesteps;
-  double CellVol, DomainSurf, DomainVol;
+  double CellVol, DomainSurf, DomainVol, CellVol_ref;
 
   protected:
   H5::H5File h5f;
@@ -61,6 +61,7 @@ class PlotterCommon
     auto attr = h5g.openAttribute(attr_name);
     notice_macro(std::string("about to read attribute value"))
     attr.read(attr.getDataType(), &ret);
+    notice_macro(std::string("attribute value read: ") + std::to_string(ret))
     return ret;
   }
 
@@ -120,6 +121,20 @@ class PlotterCommon
       h5s.getSimpleExtentDims(&n, NULL);
       map_prof.emplace("rhod", arr_prof_t(n));
       h5d.read(map_prof["rhod"].data(), H5::PredType::NATIVE_FLOAT);
+
+      // read dry air density profile on refined grid
+      try
+      {
+        h5load(file + "/const.h5", "refined rhod");
+        h5s.getSimpleExtentDims(&n, NULL);
+        map_prof.emplace("refined rhod", arr_prof_t(n));
+        h5d.read(map_prof["refined rhod"].data(), H5::PredType::NATIVE_FLOAT);
+      }
+      catch(...) // for pre-refinement simulations, use rhod as refined rhod
+      {
+        map_prof.emplace("refined rhod", map_prof["rhod"].copy());
+	std::cerr << "refined rhod as copy of rhod: " << map_prof["refined rhod"];
+      }
     }
   }
 };

UWLCM_plotters/include/PlotterMicro.hpp

@@ -20,57 +20,81 @@ class PlotterMicro_t : public Plotter_t<NDims>
   const double L_evap = 2264.76e3; // latent heat of evaporation [J/kg]
 
   public:
+  void multiply_by_rhod(arr_t &arr)
+  {
+    if(arr.extent(NDims-1) == this->map_prof["rhod"].extent(0))
+      arr *= this->map_prof["rhod"](this->LastIndex);
+    else if(arr.extent(NDims-1) == this->map_prof["refined rhod"].extent(0))
+      arr *= this->map_prof["refined rhod"](this->LastIndex);
+    else
+      throw std::runtime_error("multiply_by_rhod: input array is neither normal grid size nor refined grid size");
+  }
+
+  void multiply_by_CellVol(arr_t &arr)
+  {
+    if(arr.extent(NDims-1) == this->map_prof["rhod"].extent(0))
+      arr *= this->CellVol;
+    else if(arr.extent(NDims-1) == this->map_prof["refined rhod"].extent(0))
+      arr *= this->CellVol_ref;
+    else
+      throw std::runtime_error("multiply_by_CellVol: input array is neither normal grid size nor refined grid size");
+  }
+
   // functions for diagnosing fields
   //
   // aerosol droplets mixing ratio
   auto h5load_ra_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) //-> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
-      res = this->h5load_timestep("aerosol_rw_mom3", at) * 4./3. * 3.1416 * 1e3;
+      return arr_t(this->h5load_timestep("aerosol_rw_mom3", at) * 4./3. * 3.1416 * 1e3);
+
     else if(this->micro == "blk_1m")
+    {
       res = 0;
-    return blitz::safeToReturn(res + 0);
+      return res;
+     // return blitz::safeToReturn(res + 0);
+    }
   }
 
   // cloud droplets mixing ratio
   auto h5load_rc_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) //-> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
-      res = this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3;
+      return arr_t(this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3);
     else if(this->micro == "blk_1m")
-      res = this->h5load_timestep("rc", at);
+      return arr_t(this->h5load_timestep("rc", at));
     else if(this->micro == "blk_2m")
-      res = this->h5load_timestep("rc", at);
-    return blitz::safeToReturn(res + 0);
+      return arr_t(this->h5load_timestep("rc", at));
   }
 
   // rain droplets mixing ratio
   auto h5load_rr_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) //-> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
-      res = this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3;
+      return arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3);
     else if(this->micro == "blk_1m")
-      res = this->h5load_timestep("rr", at);
+      return arr_t(this->h5load_timestep("rr", at));
     else if(this->micro == "blk_2m")
-      res = this->h5load_timestep("rr", at);
-    return blitz::safeToReturn(res + 0);
+      return arr_t(this->h5load_timestep("rr", at));
   }
 
   // activated drops mixing ratio
   auto h5load_ract_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) 
   {
     if(this->micro == "lgrngn")
     {
-      res = this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3;
-      res += arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3);
+      return arr_t(
+	       arr_t(this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3) + 
+               arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3)
+	     );
     }
     else if(this->micro == "blk_1m")
     {
@@ -82,34 +106,38 @@ class PlotterMicro_t : public Plotter_t<NDims>
       res = this->h5load_timestep("rc", at);
       res += arr_t(this->h5load_timestep("rr", at));
     }
-    return blitz::safeToReturn(res + 0);
+   // return blitz::safeToReturn(res + 0);
+    return res;
   }
 
   // cloud droplets concentration [1/kg]
   auto h5load_nc_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) //-> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
-      res = this->h5load_timestep("cloud_rw_mom0", at);
+      return arr_t(this->h5load_timestep("cloud_rw_mom0", at));
     else if(this->micro == "blk_1m")
+    {
       res = 0;
+      return res;
+     // return blitz::safeToReturn(res + 0);
+    }
     else if(this->micro == "blk_2m")
-      res = this->h5load_timestep("nc", at);
-    return blitz::safeToReturn(res + 0);
+      return arr_t(this->h5load_timestep("nc", at));
   }
 
   // precipitation flux [W/m2]
   auto h5load_prflux_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  )// -> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
     {
-      res = this->h5load_timestep("precip_rate", at)
+      return arr_t(this->h5load_timestep("precip_rate", at)
               *  4./3 * 3.14 * 1e3 // to get mass
               / this->CellVol    // averaged over cell volume, TODO: make precip rate return specific moment? wouldnt need the dx and dy
-              * L_evap;
+              * L_evap);
     }
     else if(this->micro == "blk_1m")
       try
@@ -125,19 +153,21 @@ class PlotterMicro_t : public Plotter_t<NDims>
       {
         res = 0;
       }
-    return blitz::safeToReturn(res + 0);
+   // return blitz::safeToReturn(res + 0);
+     return res;
   }
 
   // RH
   auto h5load_RH_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) //-> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
-      res = this->h5load_timestep("RH", at);
+      return arr_t(this->h5load_timestep("RH", at));
     else if(this->micro == "blk_1m")
       res = 0;
-    return blitz::safeToReturn(res + 0);
+   // return blitz::safeToReturn(res + 0);
+     return res;
   }
 
 

drawbicyc/include/cases/RICO11/plots.hpp

@@ -1,14 +1,6 @@
 #pragma once
 
 const std::vector<std::string> series_rico({
-// "clb_bigrain_mean_rd",
-// "clb_bigrain_mean_kappa",
-// "clb_bigrain_mean_conc",
-// "clb_bigrain_mean_inclt",
-//, "clb_bigrain_mean_gccn_fraction"
-
-//"cl_acnv25_rico",
-//"cl_accr25_rico",
 "RH_max",
 "cloud_cover_rico",
 "min_cloud_base_rico",
@@ -17,14 +9,25 @@ const std::vector<std::string> series_rico({
 "lwp",
 "rwp",
 "surf_precip",
-"acc_precip"
-,"nc"
-,"nr"
-,"cl_nc_rico"
-,"cl_nr_rico"
-,"cloud_avg_supersat"
-,"wvarmax"
-,"cl_meanr" //TODO: zmienic maske na rico
+"acc_precip",
+"cl_nc_rico",
+"cl_nr_rico",
+"cloud_avg_supersat",
+"wvarmax",
+"cl_meanr", //TODO: zmienic maske na rico
+"sd_conc"
+
+
+// "clb_bigrain_mean_rd",
+// "clb_bigrain_mean_kappa",
+// "clb_bigrain_mean_conc",
+// "clb_bigrain_mean_inclt",
+//, "clb_bigrain_mean_gccn_fraction"
+
+//"cl_acnv25_rico",
+//"cl_accr25_rico",
+//,"nc"
+//,"nr"
 //TODO (po usprawnieniu cloud mask i ujednoliceniu tego:
 /*
 ,"cl_avg_cloud_rad"
@@ -34,7 +37,6 @@ const std::vector<std::string> series_rico({
 //,"rd_geq_0.8um_conc"
 //,"cl_rd_lt_0.8um_conc"
 //,"cl_rd_geq_0.8um_conc"
-,"sd_conc"
 
 /*
  "cloud_base",

drawbicyc/include/gnuplot_series_set_labels.hpp

@@ -151,7 +151,7 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt)
     gp << "set title 'average cloud drop conc [1/cm^3]'\n";
   else if (plt == "ntot")
     gp << "set title 'average particle conc [1/cm^3]'\n";
-  else if (plt == "cl_nc")
+  else if (plt == "cl_nc" || plt == "cl_nc_rico")
   {
     gp << "set title 'average cloud drop conc [1/cm^3] in cloudy cells'\n";
     gp << "set xlabel ''\n";
@@ -163,7 +163,7 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt)
     gp << "set xlabel ''\n";
     gp << "set ylabel ''\n";
   }
-  else if (plt == "cl_nr")
+  else if (plt == "cl_nr" || plt == "cl_nr_rico")
   {
     gp << "set title 'average rain drop conc [1/cm^3] in cloudy cells'\n";
     gp << "set xlabel ''\n";