Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

adding calc/plot scripts for the 8icmw_case1_1800s_sdm_vs_2dbin example

  • Loading branch information...
commit d5133ef12d6cd1b4c9321868469b22e6a0ca666b 1 parent 3fd9766
authored November 20, 2012
16  examples/8icmw_case1_1800s_sdm_vs_2dbin/bins.hpp
... ...
@@ -0,0 +1,16 @@
  1
+#pragma once
  2
+// Zach's mass-doubling bin layout
  3
+
  4
+vector<quantity<si::length, real_t>> bins()
  5
+{
  6
+  const int n_bins = 48;
  7
+  vector<quantity<si::length, real_t>> bins(n_bins);
  8
+  quantity<si::length, real_t> r = real_t(.5e-8) * si::metres;
  9
+  for (int i = 0; i < n_bins; ++i)
  10
+  {
  11
+    bins[i] = r;
  12
+    quantity<si::volume, real_t> v = real_t(4./3) * phc::pi<real_t>() * pow<3>(r);
  13
+    r = root<3>(real_t(3./4) / phc::pi<real_t>() * 2 * v);
  14
+  }
  15
+  return bins;
  16
+}
93  examples/8icmw_case1_1800s_sdm_vs_2dbin/calc.cpp
... ...
@@ -0,0 +1,93 @@
  1
+/* @file
  2
+ *  @author Sylwester Arabas <slayoo@igf.fuw.edu.pl>
  3
+ *  @author Ania Jaruga <ajaruga@igf.fuw.edu.pl>
  4
+ *  @copyright University of Warsaw
  5
+ *  @date March 2012
  6
+ *  @section LICENSE
  7
+ *    GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/)
  8
+ */
  9
+
  10
+#include <vector>
  11
+using std::vector;
  12
+
  13
+#include <string>
  14
+using std::string;
  15
+
  16
+#include <iostream>
  17
+using std::ostringstream;
  18
+using std::cerr;
  19
+using std::endl;
  20
+
  21
+#include <boost/units/systems/si.hpp>
  22
+#include <boost/units/io.hpp>
  23
+namespace si = boost::units::si;
  24
+using boost::units::quantity;
  25
+using boost::units::divide_typeof_helper;
  26
+using boost::units::detail::get_value;
  27
+
  28
+#include <boost/filesystem.hpp>
  29
+
  30
+#include "../../src/cmn/cmn_error.hpp"
  31
+#include "../../src/cmn/cmn_units.hpp"
  32
+#include "../../src/phc/phc.hpp"
  33
+#include "../../src/wkc/wkc_icmw8_case1.hpp"
  34
+
  35
+typedef double real_t;
  36
+
  37
+#include "bins.hpp"
  38
+
  39
+int main(int argc, char **argv)
  40
+{
  41
+  string dir = argc > 1 ? argv[1] : "tmp";
  42
+  notice_macro("creating " << dir << "...");
  43
+  boost::filesystem::create_directory(dir);
  44
+
  45
+  notice_macro("creating input files ...")
  46
+  string opts = wkc::icmw8_case1::create_files(
  47
+    wkc::icmw8_case1::sdm,
  48
+    dir + "/rho.nc", 
  49
+    dir + "/ini.nc",
  50
+    50, // nx
  51
+    50, // ny
  52
+    real_t(30) * si::metres, // dx
  53
+    real_t(30) * si::metres  // dy
  54
+  );
  55
+
  56
+  notice_macro("calling the solver ...")
  57
+  ostringstream cmd;
  58
+  cmd
  59
+    << "../../icicle " << opts 
  60
+    << " --adv mpdata"
  61
+      << " --adv.mpdata.fct " << true
  62
+      << " --adv.mpdata.iord " << 2
  63
+      << " --adv.mpdata.third_order " << false
  64
+    << " --dt " << real_t(1.)
  65
+    << " --nt " << real_t(1800) // 1800
  66
+    << " --nout " << real_t(60) // 60
  67
+    << " --out netcdf" 
  68
+    << " --out.netcdf.file " << dir << "/out.nc"
  69
+    << " --slv serial"
  70
+
  71
+    << " --eqs.todo_sdm.adve.sstp " << 1
  72
+    << " --eqs.todo_sdm.sedi.sstp " << 0
  73
+    << " --eqs.todo_sdm.cond.sstp " << 2
  74
+    << " --eqs.todo_sdm.chem.sstp " << 0
  75
+    << " --eqs.todo_sdm.coal.sstp " << 0
  76
+
  77
+    << " --eqs.todo_sdm.adve " << true
  78
+    << " --eqs.todo_sdm.cond " << true
  79
+    << " --eqs.todo_sdm.coal " << false
  80
+    << " --eqs.todo_sdm.sedi " << false
  81
+    << " --eqs.todo_sdm.chem " << false
  82
+    << " --eqs.todo_sdm.sd_conc_mean " << 128
  83
+
  84
+    // Zach's output setting (mass-doubling layout)
  85
+    << " --eqs.todo_sdm.out_m0 \"";
  86
+  vector<quantity<si::length, real_t>> left_edges = bins();
  87
+  for (int i = 0; i < left_edges.size()-1; ++i) 
  88
+    cmd << real_t(left_edges[i] / si::metres) << ":" << real_t(left_edges[i + 1] / si::metres) << ";";
  89
+  cmd << "\"";
  90
+    
  91
+  if (EXIT_SUCCESS != system(cmd.str().c_str()))
  92
+    error_macro("model run failed: " << cmd.str())
  93
+}
282  examples/8icmw_case1_1800s_sdm_vs_2dbin/plot.cpp
... ...
@@ -0,0 +1,282 @@
  1
+/* @file
  2
+ *  @author Sylwester Arabas <slayoo@igf.fuw.edu.pl>
  3
+ *  @copyright University of Warsaw
  4
+ *  @date March 2012
  5
+ *  @section LICENSE
  6
+ *    GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/)
  7
+ */
  8
+
  9
+#include "../../src/cmn/cmn_error.hpp"
  10
+#include "../../src/cmn/cmn_units.hpp"
  11
+#include "../../src/phc/phc.hpp"
  12
+
  13
+#include <string>
  14
+using std::string;
  15
+
  16
+#include <list>
  17
+using std::list;
  18
+
  19
+#include <vector>
  20
+using std::vector;
  21
+
  22
+#include <netcdf>
  23
+using netCDF::NcFile;
  24
+typedef vector<size_t> start;
  25
+typedef vector<size_t> count;
  26
+
  27
+#include <iostream>
  28
+using std::endl;
  29
+using std::cerr;
  30
+using std::ostringstream;
  31
+
  32
+#include <boost/format.hpp>
  33
+using boost::format;
  34
+
  35
+#include <blitz/array.h>
  36
+using blitz::Array;
  37
+using blitz::Range;
  38
+using blitz::firstDim;
  39
+using blitz::secondDim;
  40
+
  41
+#define GNUPLOT_ENABLE_BLITZ
  42
+#include <gnuplot-iostream/gnuplot-iostream.h>
  43
+
  44
+typedef double real_t;
  45
+
  46
+#include "bins.hpp"
  47
+
  48
+std::string zeropad(int n)
  49
+{
  50
+  std::ostringstream tmp;
  51
+  tmp << std::setw(5) << std::setfill('0') << n;
  52
+  return tmp.str();
  53
+}
  54
+
  55
+int main(int argc, char **argv)
  56
+{
  57
+
  58
+  string dir= argc > 1 ? argv[1] : "tmp";
  59
+ 
  60
+  notice_macro("opening super-droplet netCDF file")
  61
+  NcFile nf(dir+"/out.nc", NcFile::read);
  62
+
  63
+  notice_macro("opening super-droplet netCDF file")
  64
+  NcFile nf2("/users/arabas/devel/test2.nc", NcFile::read);
  65
+
  66
+  notice_macro("reading dt_out")
  67
+  quantity<si::time, real_t> dt_out;
  68
+  {
  69
+    real_t tmp;
  70
+    nf.getVar("dt_out").getVar(&tmp);
  71
+    dt_out = tmp * si::seconds;
  72
+  }
  73
+
  74
+  notice_macro("reading nx, ny")
  75
+  size_t 
  76
+    nt = nf.getDim("time").getSize(),
  77
+    nx = nf.getDim("X").getSize(),
  78
+    ny = nf.getDim("Y").getSize();
  79
+
  80
+  notice_macro("reading dx, dy")
  81
+  real_t dx, dy;
  82
+  {
  83
+    real_t tmp;
  84
+    nf.getVar("X").getVar(start({1}), count({1}), &dx);
  85
+    nf.getVar("X").getVar(start({0}), count({1}), &tmp);
  86
+    dx -= tmp;
  87
+    nf.getVar("Y").getVar(start({1}), count({1}), &dy);
  88
+    nf.getVar("Y").getVar(start({0}), count({1}), &tmp);
  89
+    dy -= tmp;
  90
+  }
  91
+
  92
+  string dxdy;
  93
+  {
  94
+    ostringstream tmp;
  95
+    tmp << " dx=" << dx/1000;
  96
+    tmp << " dy=" << dy/1000;
  97
+    tmp << " origin=(" << dx/2000 << "," << dy/2000 << ",0)";
  98
+    dxdy = tmp.str();
  99
+  }
  100
+
  101
+  Array<real_t, 2> rhod(nx, ny);
  102
+  {
  103
+    NcFile nfini(dir+"/ini.nc", NcFile::read);
  104
+    for (int i = 0; i < nx; ++i)
  105
+    {
  106
+      Range j(0, ny-1);
  107
+      assert(rhod(i,j).isStorageContiguous());
  108
+      nfini.getVar("rhod").getVar(start({0,0,0}), count({1,ny,1}), rhod(i,j).data());
  109
+    }
  110
+  }
  111
+
  112
+  notice_macro("setting-up plot parameters")
  113
+  Gnuplot gp;
  114
+
  115
+  for (size_t t = 0 ; t < nt; ++t) for (string &ext : list<string>({"eps","png"}))
  116
+  {
  117
+    notice_macro("generating frame at t=" << t)
  118
+    gp << "reset" << endl;
  119
+    // progressive-rock connoisseur palette ;)
  120
+    gp << "set palette defined (0 '#000000', 1 '#993399', 2 '#00CCFF', 3 '#66CC00', 4 '#FFFF00', 5 '#FC8727', 6 '#FD0000') maxcolors 20" << endl; 
  121
+    gp << "set view map" << endl;
  122
+    gp << "set tics out scale .25" << endl;
  123
+
  124
+    gp << "set ylabel 'Y [km]'" << endl;
  125
+    gp << "set yrange [" << 0 << ":" << ny * dy/1000 << "]" << endl;
  126
+    gp << "set ytics .5" << endl;
  127
+
  128
+    gp << "set contour base" << endl;
  129
+    gp << "set nosurface" << endl;
  130
+    gp << "set cntrparam levels 0" << endl;
  131
+    gp << "set nokey" << endl;
  132
+
  133
+    gp << "set label 't = " << int(real_t(t) * dt_out / si::seconds) << " s' at screen .51,.99 center" << endl;
  134
+    gp << "set label 'icicle/sdm' at screen .25,.99 center" << endl;
  135
+    gp << "set label 'Zach''s 2D-bin' at screen .75,.99 center" << endl;
  136
+
  137
+    int cols = 2, rows = 5;
  138
+
  139
+    if (ext == "png")
  140
+      gp << "set term png enhanced size " << cols * 500 << "," << rows * 500 << endl;
  141
+    else if (ext == "eps")
  142
+      gp << "set term postscript size " << 12 * cols<< "cm," << 12 * rows << "cm solid enhanced color" << endl;
  143
+    else assert(false);
  144
+
  145
+    gp << "set output '" << dir << "/test_" << zeropad(t) << "." << ext << "'" << endl;
  146
+    gp << "set multiplot layout " << rows << "," << cols << endl;
  147
+
  148
+    //// spectrum plot
  149
+    {
  150
+      vector<quantity<si::length, real_t>> left_edges = bins();
  151
+      int ns = left_edges.size() - 1;
  152
+
  153
+      gp << "set title 'mean particle concentration [1/cm^3]'" << endl;
  154
+      gp << "set logscale cb" << endl;
  155
+      gp << "set cbrange [1:200]" << endl;
  156
+      gp << "set xrange [0:" << ns << "]" << endl;
  157
+      gp << "set xtics (";
  158
+      for (int i = 0; i < left_edges.size()-1; i+=9)
  159
+        gp << (i!=0?",":"") << "\"" << format("%4.3g") % real_t(left_edges[i] / si::metres / real_t(1e-6)) << "\" " << i;
  160
+      gp << ")" << endl;
  161
+      gp << "set xlabel 'particle (wet) radius [{/Symbol m}m]'" << endl;
  162
+
  163
+      Array<real_t, 2> tmps(ns, ny), tmp(nx, ny);
  164
+
  165
+      {
  166
+        // super-droplets
  167
+        for (int i = 0; i < ns; ++i)
  168
+        {
  169
+          ostringstream name;
  170
+          name << "m_0_" << i;
  171
+          nf.getVar(name.str()).getVar(start({t,0,0,0}), count({1,nx,ny,1}), tmp.data());
  172
+          for (int y = 0; y < ny; ++y) 
  173
+            tmps(i, y) = sum(tmp(blitz::Range::all(), y)) / nx;
  174
+        }
  175
+        tmps /= 1e6; // 1/m3 -> 1/cm3
  176
+        gp << "splot '-' binary" << gp.binfmt(tmps)
  177
+          << " dx=1 dy=" << dy/1000 << " origin=(.5," << dy/2000 << ",0)"
  178
+          << " using 1 with image notitle" << endl;
  179
+        gp.sendBinary(tmps);
  180
+      }
  181
+      // 2D-bin
  182
+      { 
  183
+        for (int i = 0; i < ns; ++i)
  184
+        {
  185
+          nf2.getVar("QCDIST").getVar(start({t,0,0,i}), count({1,ny,nx,1}), tmp.data());
  186
+          for (int y = 0; y < ny; ++y) 
  187
+            tmps(i, y) = sum(tmp(y, blitz::Range::all())) / nx;
  188
+        }
  189
+        nf2.getVar("GAC").getVar(start({t,0,0}), count({1,ny,nx}), tmp.data());
  190
+        tmp.transposeSelf(secondDim, firstDim);
  191
+        tmps *= tmp; // kg-1 -> m-3
  192
+        tmps /= 1e6; // 1/m3 -> 1/cm3
  193
+        gp << "splot '-' binary" << gp.binfmt(tmps)
  194
+          << " dx=1 dy=" << dy/1000 << " origin=(.5," << dy/2000 << ",0)"
  195
+          << " using 1 with image notitle" << endl;
  196
+        gp.sendBinary(tmps);
  197
+      }
  198
+    }
  199
+
  200
+    gp << "set xlabel 'X [km]'" << endl;
  201
+    gp << "set xrange [" << 0 << ":" << nx * dx/1000 << "]" << endl;
  202
+    gp << "set xtics .5" << endl;
  203
+    gp << "unset logscale cb" << endl;
  204
+
  205
+    //// water vapour mixing ratio
  206
+    {
  207
+      Array<real_t, 2> tmp(nx, ny);
  208
+
  209
+      gp << "set title 'water vapour mixing ratio [g/kg]'" << endl;
  210
+      gp << "set cbrange [6:8]" << endl;
  211
+      // super-droplets
  212
+      nf.getVar("rhod_rv").getVar(start({t,0,0,0}), count({1,nx,ny,1}), tmp.data()); 
  213
+      tmp /= rhod;
  214
+      gp << "splot '-' binary" << gp.binfmt(tmp) << dxdy << " using ($1*1000) with image notitle" << endl;
  215
+      gp.sendBinary(tmp);
  216
+      // 2D-bin
  217
+      nf2.getVar("QV").getVar(start({t,0,0}), count({1,ny,nx}), tmp.data()); 
  218
+      tmp.transposeSelf(secondDim, firstDim);
  219
+      gp << "splot '-' binary" << gp.binfmt(tmp) << dxdy << " using ($1*1000) with image notitle" << endl;
  220
+      gp.sendBinary(tmp);
  221
+    }
  222
+
  223
+
  224
+    //// potential temperature
  225
+    {
  226
+      Array<real_t, 2> tmp(nx, ny);
  227
+
  228
+      gp << "set title 'potential temperature [K]'" << endl;
  229
+      gp << "set cbrange [288:293]" << endl;
  230
+      // super-droplet
  231
+      nf.getVar("rhod_th").getVar(start({t,0,0,0}), count({1,nx,ny,1}), tmp.data()); 
  232
+      tmp /= rhod;
  233
+      gp << "splot '-' binary" << gp.binfmt(tmp) << dxdy << " with image notitle" << endl;
  234
+      gp.sendBinary(tmp);
  235
+      // 2D-bin
  236
+      nf2.getVar("THETA").getVar(start({t,0,0}), count({1,ny,nx}), tmp.data()); 
  237
+      tmp.transposeSelf(secondDim, firstDim);
  238
+      gp << "splot '-' binary" << gp.binfmt(tmp) << dxdy << " with image notitle" << endl;
  239
+      gp.sendBinary(tmp);
  240
+    }
  241
+
  242
+ 
  243
+    //// super-droplet concentration
  244
+    {
  245
+      Array<real_t, 2> tmp(nx, ny);
  246
+
  247
+      gp << "set title 'super-droplet conc. [1/dx/dy/dz]'" << endl;
  248
+      gp << "set cbrange [0:512]" << endl;
  249
+      nf.getVar("sd_conc").getVar(start({t,0,0,0}), count({1,nx,ny,1}), tmp.data()); 
  250
+      gp << "splot '-' binary" << gp.binfmt(tmp) << dxdy << " using 1 with image notitle";
  251
+      gp << endl;
  252
+      gp.sendBinary(tmp);
  253
+      // 2D-bin
  254
+      gp << "splot 0" << endl;
  255
+    }
  256
+
  257
+
  258
+    // effective radius
  259
+    {
  260
+      Array<real_t, 2> tmp0(nx, ny), tmp1(nx, ny);
  261
+
  262
+      gp << "set title 'effective radius [{/Symbol m}m] (FSSP range)'" << endl;
  263
+      gp << "set cbrange [1:30]" << endl;
  264
+      nf.getVar("m_3").getVar(start({t,0,0,0}), count({1,nx,ny,1}), tmp0.data()); 
  265
+      nf.getVar("m_2").getVar(start({t,0,0,0}), count({1,nx,ny,1}), tmp1.data()); 
  266
+      tmp0 /= tmp1;
  267
+      tmp0 *= 1e6;
  268
+      gp << "splot '-' binary" << gp.binfmt(tmp0) << dxdy << " using 1 with image notitle";
  269
+      gp << endl;
  270
+      gp.sendBinary(tmp0);
  271
+      // 2D-bin
  272
+      gp << "splot 0" << endl;
  273
+    }
  274
+
  275
+    gp << "unset multiplot" << endl;
  276
+    gp << "unset label" << endl;
  277
+  }
  278
+
  279
+  string cmd="convert -monitor -delay 10 -loop 1 " + dir + "/test_*.png " + dir + "/todo.gif 1>&2";
  280
+  int status = system(cmd.c_str());
  281
+  notice_macro("done.")
  282
+}

0 notes on commit d5133ef

Please sign in to comment.
Something went wrong with that request. Please try again.