wip on removal of reindex in force calculations

drawbicyc/PlotterMicro.hpp

@@ -20,6 +20,18 @@ class PlotterMicro_t : public Plotter_t<NDims>
   public:
   // functions for diagnosing fields
   //
+  // aerosol droplets mixing ratio
+  auto h5load_ra_timestep(
+    int at
+  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  {
+    if(this->micro == "lgrngn")
+      res = 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);
+  }
+
   // cloud droplets mixing ratio
   auto h5load_rc_timestep(
     int at
@@ -40,7 +52,7 @@ class PlotterMicro_t : public Plotter_t<NDims>
     if(this->micro == "lgrngn")
       res = this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3;
     else if(this->micro == "blk_1m")
-      res = this->h5load_timestep("rc", at);
+      res = this->h5load_timestep("rr", at);
     return blitz::safeToReturn(res + 0);
   }
 
@@ -59,6 +71,48 @@ class PlotterMicro_t : public Plotter_t<NDims>
     return blitz::safeToReturn(res + 0);
   }
 
+  // cloud droplets concentration [1/kg]
+  auto h5load_nc_timestep(
+    int at
+  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  {
+    if(this->micro == "lgrngn")
+      res = this->h5load_timestep("cloud_rw_mom0", at);
+    else if(this->micro == "blk_1m")
+      res = 0;
+    return blitz::safeToReturn(res + 0);
+  }
+
+  // precipitation flux [W/m2]
+  auto h5load_prflux_timestep(
+    int at
+  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  {
+    if(this->micro == "lgrngn")
+    {
+      res = 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
+              * 2264.76e3;         // latent heat of evaporation [J/kg]
+    }
+    else if(this->micro == "blk_1m")
+      res = 0;
+    return blitz::safeToReturn(res + 0);
+  }
+
+  // RH
+  auto h5load_RH_timestep(
+    int at
+  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  {
+    if(this->micro == "lgrngn")
+      res = this->h5load_timestep("RH", at);
+    else if(this->micro == "blk_1m")
+      res = 0;
+    return blitz::safeToReturn(res + 0);
+  }
+
+
   // functions for diagnosing statistics
 
   // helper function that calculates staistics (mean and std_dev) of a field only in cloudy cells

drawbicyc/drawbicyc.cpp

@@ -56,6 +56,7 @@ int main(int argc, char** argv)
 
   if(NDims == 2)
   {
+
     if(flag_series)   plot_series(PlotterMicro_t<2>(h5, micro), Plots(type));
     if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), Plots(type));
     if(flag_fields)   plot_fields(PlotterMicro_t<2>(h5, micro), Plots(type));
@@ -68,8 +69,8 @@ int main(int argc, char** argv)
     if(flag_fields)   plot_fields(PlotterMicro_t<3>(h5, micro), Plots(type));
     if(flag_qv_qc_2_6_10_min)   plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro), Plots(type));
     if(flag_lgrngn_spec) {
-      plot_lgrngn_spec_positions(PlotterMicro_t<3>(h5, "lgrngn"), Plots(type), 12);
-      plot_lgrngn_spec(PlotterMicro_t<3>(h5, "lgrngn"), Plots(type), 12);
+      plot_lgrngn_spec_positions(PlotterMicro_t<3>(h5, "lgrngn"), Plots(type));
+      plot_lgrngn_spec(PlotterMicro_t<3>(h5, "lgrngn"), Plots(type));
     }
   }
   else

drawbicyc/plot_fields.hpp

@@ -31,7 +31,7 @@ void plot_fields(Plotter_t plotter, Plots plots)
     {
       std::cout << at * n["outfreq"] << " : " << plt << std::endl;
       Gnuplot gp;
-      init(gp, plotter.file + ".plot/" + plt + "/" + zeropad(at * n["outfreq"]) + ".svg", 1, 1, n); 
+      init(gp, plotter.file + ".plot/" + plt + "/" + zeropad(at * n["outfreq"]) + ".svg", 1, 1, n, 1, 0.666666); 
 
       if (plt == "rl")
       {
@@ -205,6 +205,26 @@ void plot_fields(Plotter_t plotter, Plots plots)
         }
         catch(...){}
       }   
+      else if (plt == "lib_pres")
+      {   
+        try{
+	std::string title = "libcloud pressure [Pa]"; 
+	gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
+        auto tmp = plotter.h5load_timestep("libcloud_pressure", at * n["outfreq"]);
+        plotter.plot(gp, tmp);
+        }
+        catch(...){}
+      }   
+      else if (plt == "lib_temp")
+      {   
+        try{
+	std::string title = "libcloud temperature [K]"; 
+	gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
+        auto tmp = plotter.h5load_timestep("libcloud_temperature", at * n["outfreq"]);
+        plotter.plot(gp, tmp);
+        }
+        catch(...){}
+      }   
       else if (plt == "supersat")
       {   
         try{

drawbicyc/plot_lgrngn_spec.hpp

@@ -4,21 +4,33 @@
 #include "bins.hpp"
 #include <map>
 
+namespace
+{
+  po::variables_map vm;
+};
 
 //plot spectra positions
 template<class Plotter_t>
-void plot_lgrngn_spec_positions(Plotter_t plotter, Plots plots, int at)
+void plot_lgrngn_spec_positions(Plotter_t plotter, Plots plots)
 {
+  // read opts
+  po::options_description opts("spectra plotting options");
+  opts.add_options()
+    ("spectra_step", po::value<int>()->required() , "time step number at which spectra are plotted")
+  ;
+  handle_opts(opts, vm);
+  const int spectra_step = vm["spectra_step"].as<int>();
+
   auto& n = plotter.map;
   Gnuplot gp;
 
   try{
   // cloud water content
-  auto tmp = plotter.h5load_ract_timestep(at * n["outfreq"]) * 1e3;
+  auto tmp = plotter.h5load_ract_timestep(spectra_step * n["outfreq"]) * 1e3;
 
   std::string title = "cloud water mixing ratio [g/kg]";
-  gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
-  init(gp, plotter.file + "_spectra_positions_at" + zeropad(at) + ".svg", 1, 1, n, 2.); 
+  gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(spectra_step) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
+  init(gp, plotter.file + "_spectra_positions_at" + zeropad(spectra_step) + ".svg", 1, 1, n, 2.); 
 
   {   
     for (auto &fcs : focus_3d)
@@ -60,8 +72,10 @@ void plot_lgrngn_spec_positions(Plotter_t plotter, Plots plots, int at)
 
 //plot spectra
 template<class Plotter_t>
-void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at)
+void plot_lgrngn_spec(Plotter_t plotter, Plots plots)
 {
+  const int spectra_step = vm["spectra_step"].as<int>();
+
   const double r_c_adiab = 7.8; // adiabatic water content [g/kg] coming from Wojtek
   auto& n = plotter.map;
   for(auto elem : n)
@@ -71,38 +85,45 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at)
   Gnuplot gp;
 //  int off = 2; // TODO!!!
   int off = 0; // TODO!!!
-  string file = plotter.file + "_spectra_at" + zeropad(at) + ".svg";
+  string file = plotter.file + "_spectra_at" + zeropad(spectra_step) + ".svg";
 
-  int hor = min<int>(focus_3d.size(), 2);
-  int ver = double(focus_3d.size()) / 2. + 0.99999;
+  int hor = min<int>(focus_3d.size(), 5);
+  int ver = double(focus_3d.size()) / 5. + 0.99999;
 
   init_prof(gp, file, ver, hor);
 
 
   gp << "set xrange [0:30]\n";
   gp << "set yrange [0:50]\n";
   gp << "set xlabel 'r [μm]'\n"; 
-  gp << "set ylabel 'n [cm^{-3} μm^{-1}]'\n"; 
 
   // read in density
   auto tmp = plotter.h5load(plotter.file + "/const.h5", "G");
   typename Plotter_t::arr_t rhod(tmp);
 
   // focus to the gridbox from where the size distribution is plotted
   char lbl = 'a';
+  bool setylabel = 1;
   for (auto &fcs : focus_3d)
   {
+    if(setylabel)
+    {
+      gp << "set ylabel 'n [cm^{-3} μm^{-1}]'\n"; 
+      setylabel = 0;
+    }
+    else
+      gp << "unset ylabel\n"; 
     const int &x = fcs[0], &y = fcs[1], &z = fcs[2];
     const blitz::RectDomain<3> focusBox({x-box_size, y-box_size, z-box_size}, {x+box_size, y+box_size, z+box_size});
 
-    gp << "set label 1 '(" << lbl << ")' at graph .1, .93 font ',20'\n";
+    gp << "set label 1 '(" << lbl << ")' at graph .1, .93 font \",15\"\n";
     lbl += 1;
 
     // calc ratio of water content to adiabatic water content
     {
-      auto tmp = plotter.h5load_ract_timestep(at * n["outfreq"]) * 1e3;
+      auto tmp = plotter.h5load_ract_timestep(spectra_step * n["outfreq"]) * 1e3;
       double ratio = mean(tmp(focusBox)) / r_c_adiab;
-      gp << "set label 4 'r_c / r_c^{adiab} = " << std::setprecision(2) << ratio << "' at graph .2, .93 font ',20'\n";
+      gp << "set label 4 'AF = " << std::fixed << std::setprecision(2) << ratio << "' at graph .2, .63 font \",15\"\n";
     }
 
     // calc mean and std dev of radius of acivated droplets in the box
@@ -111,15 +132,15 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at)
     double act_rw_std_dev = 0.; 
     try
     {
-      auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]);
+      auto tmp = plotter.h5load_timestep("actrw_rw_mom0", spectra_step * n["outfreq"]);
       typename Plotter_t::arr_t snap(tmp); 
       act_conc = mean(snap(focusBox));
     }
     catch(...) {;}
     // mean droplet radius in the box
     try
     {
-      auto tmp = plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]);
+      auto tmp = plotter.h5load_timestep("actrw_rw_mom1", spectra_step * n["outfreq"]);
       typename Plotter_t::arr_t snap(tmp); // 1st raw moment / mass [m / kg]
       if(act_conc > 0)
         act_rw_mean = mean(snap(focusBox)) / act_conc;
@@ -130,13 +151,13 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at)
     // std deviation of distribution of radius in the box
     try
     {
-      auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]);
+      auto tmp = plotter.h5load_timestep("actrw_rw_mom0", spectra_step * n["outfreq"]);
       typename Plotter_t::arr_t zeroth_raw_mom(tmp); // 0th raw moment / mass [1 / kg]
-      tmp = plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]);
+      tmp = plotter.h5load_timestep("actrw_rw_mom1", spectra_step * n["outfreq"]);
       typename Plotter_t::arr_t first_raw_mom(tmp); // 1st raw moment / mass [m / kg]
-      tmp = plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]);
+      tmp = plotter.h5load_timestep("actrw_rw_mom2", spectra_step * n["outfreq"]);
       typename Plotter_t::arr_t second_raw_mom(tmp); // 2nd raw moment / mass [m^2 / kg]
-      tmp = plotter.h5load_timestep("sd_conc", at * n["outfreq"]);
+      tmp = plotter.h5load_timestep("sd_conc", spectra_step * n["outfreq"]);
       typename Plotter_t::arr_t sd_conc(tmp); // number of SDs
       if(act_conc > 0)
       {
@@ -165,8 +186,8 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at)
     catch(...) {;}
 
 
-    gp << "set label 2 '<r>_{act} = " << std::setprecision(2) << act_rw_mean * 1e6 <<"µm' at graph .1, .83 font '15'\n";
-    gp << "set label 3 'σ(r)_{act} = " << std::setprecision(2) << act_rw_std_dev * 1e6 <<"µm' at graph .1, .73 font ',15'\n";
+    gp << "set label 2 '<r> = " << std::fixed << std::setprecision(2) << act_rw_mean * 1e6 <<" µm' at graph .2, .83 font \",15\"\n";
+    gp << "set label 3 'σ = " << std::fixed << std::setprecision(2) << act_rw_std_dev * 1e6 <<" µm' at graph .2, .73 font \",15\"\n";
 
 
     std::map<float, float> focus_w;
@@ -178,7 +199,7 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at)
     for (int i = 0; i < nsw; ++i)
     {
       const string name = "rw_rng" + zeropad(i + off) + "_mom0";
-      auto tmp_w = plotter.h5load_timestep(name, at * n["outfreq"]);
+      auto tmp_w = plotter.h5load_timestep(name, spectra_step * n["outfreq"]);
 
       focus_w[(left_edges_rw[i] + left_edges_rw[i+1]) / 2. / 1e-6 / si::metres] =
         sum((tmp_w*rhod)(focusBox))
@@ -187,7 +208,7 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at)
         / ((left_edges_rw[i+1] - left_edges_rw[i]) / 1e-6 / si::metres); // per micrometre
     }
     const string name = "rw_rng" + zeropad(nsw + off) + "_mom0";
-    auto tmp_w = plotter.h5load_timestep(name, at * n["outfreq"]);
+    auto tmp_w = plotter.h5load_timestep(name, spectra_step * n["outfreq"]);
 
     notice_macro("setting-up plot parameters");
     std::cout << "larger drops conc: " << (tmp_w * rhod)(x,y,z) * 1e-6 << endl;