Merge pull request #273 from trontrytel/thesis_plots

Thesis plots

Author: pdziekan

src/particles_impl_chem_henry.ipp

@@ -191,9 +191,9 @@ namespace libcloudphxx
 
           // implicit solution to the eq. 8.22 from chapter 8.4.2 
           // in Peter Warneck Chemistry of the Natural Atmosphere  
-          return
-            (
-              ( m_old 
+          real_t mass_helper =     
+          (
+            ( m_old 
                 + 
                 (dt * si::seconds) * (V * si::cubic_metres) 
                 * common::henry::mass_trans(
@@ -210,7 +210,9 @@ namespace libcloudphxx
                  * common::henry::mass_trans(rw2, (D * si::metres * si::metres / si::seconds), 
                                             acc_coeff * si::seconds / si::seconds, T, (M_gas * si::kilograms / si::moles)) 
                  / Henry / common::moist_air::kaBoNA<real_t>() / T)
-              ) / si::kilograms;
+          ) / si::kilograms ;
+
+          return mass_helper;// > 0 ? mass_helper : 0;
         }
       };
     };

src/particles_impl_chem_react.ipp

@@ -19,10 +19,11 @@ namespace libcloudphxx
       { // returns the change in mass per second of each chemical compounds
         // due to oxidation reaction
         const int chem_iter;
+        const real_t dt;
 
         // ctor
-        chem_rhs_helper(const int &chem_iter) : 
-          chem_iter(chem_iter)
+        chem_rhs_helper(const int &chem_iter, const real_t dt) : 
+          chem_iter(chem_iter), dt(dt)
         {}
 
         BOOST_GPU_ENABLED
@@ -67,7 +68,14 @@ namespace libcloudphxx
             O3_react = V * m_O3 / M_O3<real_t>() / V * m_S_IV / M_SO2_H2O<real_t>() / V
                        / (real_t(1) + Kt_SO2 / conc_H + Kt_SO2 * Kt_HSO3 / conc_H / conc_H)
                        * (R_O3_k0 + R_O3_k1 * Kt_SO2 / conc_H + R_O3_k2 * Kt_SO2 * Kt_HSO3 / conc_H / conc_H);
-       
+            
+            // check if reaction rate won't take more O3 than there is
+            O3_react = (O3_react * dt * si::seconds < m_O3 / M_O3<real_t>()) ? O3_react : 
+                         m_O3 / M_O3<real_t>() / (dt * si::seconds);
+            // check if reaction rate won't take more S_IV than there is
+            O3_react = (O3_react * dt * si::seconds < m_S_IV / M_SO2_H2O<real_t>()) ? O3_react : 
+                         m_S_IV / M_SO2_H2O<real_t>() / (dt * si::seconds);
+ 
 	  // helper for H2O2 reaction
 	  quantity<divide_typeof_helper<si::amount, si::time>::type, real_t>
             H2O2_react = V * R_H2O2_k * Kt_SO2 
@@ -76,6 +84,13 @@ namespace libcloudphxx
                          / (real_t(1) + Kt_SO2 / conc_H + Kt_SO2 * Kt_HSO3 / conc_H / conc_H)
                          / (real_t(1) + R_S_H2O2_K<real_t>() * conc_H);
 
+            // check if reaction rate won't take more H2O2 than there is
+            H2O2_react = (H2O2_react * dt * si::seconds < m_H2O2 / M_H2O2<real_t>()) ? H2O2_react : 
+                           m_H2O2 / M_H2O2<real_t>() / (dt * si::seconds);
+            // check if reaction rate won't take more S_IV than there is (this silently gives precedence to O3 reaction)
+            H2O2_react = (H2O2_react * dt * si::seconds < m_S_IV / M_SO2_H2O<real_t>() - O3_react * dt * si::seconds) ? 
+                          H2O2_react : m_S_IV / M_SO2_H2O<real_t>() / (dt * si::seconds) - O3_react;
+
           switch (chem_iter)
           {
             case SO2:
@@ -104,7 +119,8 @@ namespace libcloudphxx
       // functor called by odeint
       template <typename real_t>
       struct chem_rhs
-      {
+      { 
+        const real_t dt;
         const thrust_device::vector<real_t> &V;
         const thrust_device::vector<unsigned int> &chem_flag;
 
@@ -119,12 +135,13 @@ namespace libcloudphxx
 
         // ctor
         chem_rhs(
+          const real_t dt,
           const thrust_device::vector<real_t> &V,
           const pi_t &T,
           const typename thrust_device::vector<real_t>::const_iterator &m_H,
           const thrust_device::vector<unsigned int> &chem_flag
         ) :
-          V(V), T(T), m_H(m_H), n_part(V.size()), chem_flag(chem_flag)
+          dt(dt), V(V), T(T), m_H(m_H), n_part(V.size()), chem_flag(chem_flag)
         {}
 
         void operator()(
@@ -179,7 +196,7 @@ namespace libcloudphxx
                   // output
                   dot_psi.begin() + (chem_iter - chem_rhs_beg) * n_part, 
                   // op
-                  chem_rhs_helper<real_t>(chem_iter),
+                  chem_rhs_helper<real_t>(chem_iter, dt),
                   // condition
                   thrust::identity<unsigned int>()
                 );
@@ -259,6 +276,7 @@ namespace libcloudphxx
       // do chemical reactions
       chem_stepper.do_step(
         detail::chem_rhs<real_t>(
+          dt,
           V,
           thrust::make_permutation_iterator(T.begin(), ijk.begin()), 
           chem_bgn[H], 

src/particles_pimpl_ctor.ipp

@@ -162,7 +162,7 @@ namespace libcloudphxx
         > runge_kutta_fehlberg78 or the bulirsch_stoer with a very high order. But
         > should benchmark both steppers and choose the faster one.
       */
-      boost::numeric::odeint::euler<
+      boost::numeric::odeint::runge_kutta4<
         thrust_device::vector<real_t>, // state_type
         real_t,                        // value_type
         thrust_device::vector<real_t>, // deriv_type

src/particles_step.ipp

@@ -138,15 +138,21 @@ namespace libcloudphxx
             //cleanup - TODO think of something better
             pimpl->chem_cleanup();
           }
-              
+ 
           if (opts.chem_dsc)
           { //dissociation
             pimpl->chem_dissoc();
+
+            //cleanup - TODO think of something better
+            pimpl->chem_cleanup();
           }
-            
+           
           if (opts.chem_rct)
           { //oxidation 
             pimpl->chem_react(pimpl->opts_init.dt / pimpl->opts_init.sstp_chem);
+
+            //cleanup - TODO think of something better
+            pimpl->chem_cleanup();
           }
         }
       }

tests/paper_GMD_2015/src/kin_cloud_2d_lgrngn.hpp

@@ -70,6 +70,18 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common<ct_params_t>
         rng_num++;
       }
     }
+    {
+      // rw3(rd)
+      int rng_num = 0;
+      for (auto &rng_moms : params.out_dry)
+      {
+        auto &rng(rng_moms.first);
+        prtcls->diag_dry_rng(rng.first / si::metres, rng.second / si::metres);
+        prtcls->diag_wet_mom(3);
+        this->record_aux(aux_name("rw3ofrd", rng_num, 3), prtcls->outbuf());
+        rng_num++;
+      }
+    }
   } 
 
   libcloudphxx::lgrngn::arrinfo_t<real_t> make_arrinfo(

tests/paper_GMD_2015/tests/chem_sandbox/calc_chem.cpp

@@ -14,7 +14,7 @@ int main(int ac, char** av)
 {
   if (ac != 2) error_macro("expecting one argument - CMAKE_BINARY_DIR");
 
-  string bins_wet_str, bins_dry_str;
+  string bins_wet_str, bins_dry_str, bins_wet_str_2;
 
   {
     ostringstream tmp;
@@ -32,31 +32,41 @@ int main(int ac, char** av)
     bins_wet_str = tmp.str();
   }
 
+  {
+    ostringstream tmp;
+    vector<quantity<si::length>> left_edges = bins_wet();
+    for (int i = 0; i < left_edges.size()-1; ++i)
+      tmp << float(left_edges[i] / si::metres) << ":" << float(left_edges[i + 1] / si::metres) << "|0,3;";
+    bins_wet_str_2 = tmp.str();
+  }
+
+
   for (const string &kernel : {"hall_pinsky_stratocumulus"})
-                           /*   {"hall", "hall_davis_no_waals", 
+      /*                        {"hall", "hall_davis_no_waals", 
                                "onishi_hall", "onishi_hall_davis_no_waals", 
                                "vohl_davis_no_waals", 
                                "hall_pinsky_stratocumulus"
                               }) */
   {
 
     string opts_common = 
-      //"--outfreq=11800  --nt=11800 --spinup=10000 --nx=76 --nz=76 --relax_th_rv=false --rng_seed=44 ";
-      "--outfreq=200  --nt=11800 --spinup=10000 --nx=76 --nz=76 --relax_th_rv=false --rng_seed=44 ";
+      //"--outfreq=11800  --dt=1 --nt=11800 --spinup=10000 --nx=76 --nz=76 --relax_th_rv=false --rng_seed=44 ";
+      "--outfreq=200  --dt=1 --nt=11800 --spinup=10000 --nx=76 --nz=76 --relax_th_rv=false --rng_seed=44 ";
     set<string> opts_micro({
-      "--micro=lgrngn_chem --outdir=out_"+kernel+" --backend=CUDA --adv_serial=False --sd_conc=100 "
+      "--micro=lgrngn_chem --outdir=out_"+kernel+" --backend=CUDA --adv_serial=False --sd_conc=256 "
                       "--sstp_cond=10 --coal=True --sedi=True " 
                       "--w_max=.6 "
                       "--chem_switch=True --chem_dsl=True --chem_dsc=True --chem_rho=1.8e3 --sstp_chem=10 "  
                       //chem_rct switched on afetr spinup in set_chem  
                       "--mean_rd1=0.05e-6 --sdev_rd1=1.8 --n1_stp=50e6 "
                       "--mean_rd2=0.1e-6  --sdev_rd2=1.5 --n2_stp=0 "
                       "--kernel="+kernel+" --terminal_velocity=beard77fast "
-                      "--SO2_g_0=.2e-9 --O3_g_0=50e-9 --H2O2_g_0=.5e-9 --CO2_g_0=360e-6 --NH3_g_0=.1e-9 --HNO3_g_0=.1e-9 "
+                      //"--SO2_g_0=.2e-9 --O3_g_0=25e-9 --H2O2_g_0=.4e-9 --CO2_g_0=360e-6 --NH3_g_0=.4e-9 --HNO3_g_0=.1e-9 "
+                      "--SO2_g_0=.2e-9 --O3_g_0=25e-9 --H2O2_g_0=.4e-9 --CO2_g_0=360e-6 --NH3_g_0=.1e-9 --HNO3_g_0=.1e-9 "
         " --out_wet=\""
           ".5e-6:25e-6|0,1,2,3;" // FSSP
           "25e-6:1|0,3;"         // "rain"
-          + bins_wet_str +       // aerosol spectrum (wet)
+          + bins_wet_str_2 +       // aerosol spectrum (wet)
         "\""
         " --out_dry=\""
           "0.:1.|0,1;"
@@ -67,7 +77,7 @@ int main(int ac, char** av)
         "\""
         " --out_wet_pH=\""
           + bins_wet_str +       // spectrum for S_VI and H+ (wet)
-        "\""
+        "\""   
     });
 
     for (auto &opts_m : opts_micro)

tests/paper_GMD_2015/tests/chem_sandbox/chem_plots/chem_size_distr.py

@@ -0,0 +1,144 @@
+"""
+Plots size distributions from updraft and downdraft 
+    from above cloud base and from below cloud base
+"""
+import numpy   as np
+import h5py    as h5
+import Gnuplot as gp
+import math    as mt
+
+# all the collision kernels
+kernels = {"out_hall_pinsky_stratocumulus"}
+
+ub_cover = {"out_hall_pinsky_stratocumulus" : 0}
+ua_cover = {"out_hall_pinsky_stratocumulus" : 0}
+db_cover = {"out_hall_pinsky_stratocumulus" : 0}
+da_cover = {"out_hall_pinsky_stratocumulus" : 0}
+
+# left and right edges of bins for dry and wet radius
+dry_edges = np.fromiter( (1e-6 * 10**(-3 + i * .1) for i in xrange(41)), dtype="float")
+wet_edges = np.fromiter( (1e-6 * 10**(-3 + i * .1) for i in xrange(56)), dtype="float")
+num_dry   = dry_edges.size - 1
+num_wet   = wet_edges.size - 1
+
+for kernel in kernels:
+
+    # helper arrays for storing initial and cutoff dry radii as well as final and cutoff wet radii
+    n_dry_ub = np.zeros(num_dry)
+    n_dry_ua = np.zeros(num_dry)
+    n_dry_db = np.zeros(num_dry)
+    n_dry_da = np.zeros(num_dry)
+
+    n_wet_ub = np.zeros(num_wet)
+    n_wet_ua = np.zeros(num_wet)
+    n_wet_db = np.zeros(num_wet)
+    n_wet_da = np.zeros(num_wet)
+
+    for run in {"seed_44"}:
+        # open hdf5 files with data
+        h5f_cst   = h5.File('data/case_base/' + kernel + '/const.h5', 'r')
+        h5f_ini   = h5.File('data/case_base/' + kernel + '/timestep0000000000.h5', 'r')
+        h5f_fin   = h5.File('data/case_base/' + kernel + '/timestep0000011800.h5', 'r')
+
+        x_grid = h5f_cst['X'][:-1][:-1]
+        y_grid = h5f_cst['Y'][:-1][:-1]
+        
+        idx_ub = np.where((x_grid > 1 ) & (x_grid < 37) & (y_grid > 15) & (y_grid < 45))
+        idx_ua = np.where((x_grid > 1 ) & (x_grid < 37) & (y_grid > 45) & (y_grid < 70))
+        idx_db = np.where((x_grid > 39) & (x_grid < 76) & (y_grid > 15) & (y_grid < 45))
+        idx_da = np.where((x_grid > 39) & (x_grid < 76) & (y_grid > 45) & (y_grid < 70))
+
+        # count dry cloudy and rainy grid cells
+        ub_cover[kernel] += idx_ub[0].size
+        ua_cover[kernel] += idx_ua[0].size
+        db_cover[kernel] += idx_db[0].size
+        da_cover[kernel] += idx_da[0].size
+
+        for i in range(num_dry-1): # first bin id for total conc
+            name = "rd_rng" + str(i+1).zfill(3) + "_mom0"
+            tmp_ini  = 1e-6 * h5f_ini[name][:]
+            tmp_fin  = 1e-6 * h5f_fin[name][:]
+
+            n_dry_ub[i] += tmp_fin[idx_ub].sum() / idx_ub[0].size / (mt.log(dry_edges[i+1], 10) - mt.log(dry_edges[i], 10))
+            n_dry_ua[i] += tmp_fin[idx_ua].sum() / idx_ua[0].size / (mt.log(dry_edges[i+1], 10) - mt.log(dry_edges[i], 10))
+            n_dry_db[i] += tmp_fin[idx_db].sum() / idx_db[0].size / (mt.log(dry_edges[i+1], 10) - mt.log(dry_edges[i], 10))
+            n_dry_da[i] += tmp_fin[idx_da].sum() / idx_da[0].size / (mt.log(dry_edges[i+1], 10) - mt.log(dry_edges[i], 10))
+   
+        for i in range(num_wet-2): # first two bins are for total cloud and rain water conc
+            name = "rw_rng" + str(i+2).zfill(3) + "_mom0"
+            tmp_fin  = 1e-6 * h5f_fin[name][:]
+
+            n_wet_ub[i] += tmp_fin[idx_ub].sum() / idx_ub[0].size / (mt.log(wet_edges[i+1], 10) - mt.log(wet_edges[i], 10))
+            n_wet_ua[i] += tmp_fin[idx_ua].sum() / idx_ua[0].size / (mt.log(wet_edges[i+1], 10) - mt.log(wet_edges[i], 10))
+            n_wet_db[i] += tmp_fin[idx_db].sum() / idx_db[0].size / (mt.log(wet_edges[i+1], 10) - mt.log(wet_edges[i], 10))
+            n_wet_da[i] += tmp_fin[idx_da].sum() / idx_da[0].size / (mt.log(wet_edges[i+1], 10) - mt.log(wet_edges[i], 10))
+    
+        # close hdf5 files
+        h5f_ini.close()
+        h5f_fin.close()
+
+    ub_cover[kernel] = ub_cover[kernel] / 76. / 76 * 100 
+    ua_cover[kernel] = ua_cover[kernel] / 76. / 76 * 100 
+    db_cover[kernel] = db_cover[kernel] / 76. / 76 * 100 
+    da_cover[kernel] = da_cover[kernel] / 76. / 76 * 100 
+
+    print kernel, "  ", ub_cover[kernel]
+    print kernel, "  ", ua_cover[kernel]
+    print kernel, "  ", db_cover[kernel]
+    print kernel, "  ", da_cover[kernel]
+
+    # plotting
+    ymin = 0.001
+    ymax = 1000
+    xmin = 0.001
+    xmax = 5
+
+    g = gp.Gnuplot()
+    g('reset')
+    g('set term svg dynamic enhanced font "Verdana, 14"')
+    g('set output "plots/' + kernel + '_aerosol_ud.svg" ')
+    g('set logscale xy')
+    g('set key samplen 1.2')
+    g('set xtics rotate by 65 right (.01, .1, 1, 10, 100)')
+    g('set ytics (.001, .01, .1, 1, 10, 100, 1000)')
+    g('set xlabel "particle radius [um]"')
+    g('set ylabel "dN/dlog_{10}(r) [mg^{-1} per log_{10}(size interval)]"')
+    g('set xrange [' +  str(xmin) + ':' + str(xmax) + ']')
+    g('set yrange [' +  str(ymin) + ':' + str(ymax) + ']')
+    g('set grid')
+    g('set nokey')
+    
+    plot_rd_ub  = gp.PlotItems.Data(dry_edges[:-1] * 1e6 , n_dry_ub,  with_="steps lw 4 lc rgb 'black'")
+    plot_rd_ua  = gp.PlotItems.Data(dry_edges[:-1] * 1e6 , n_dry_ua,  with_="steps lw 4 lc rgb 'green'")
+    plot_rd_db  = gp.PlotItems.Data(dry_edges[:-1] * 1e6 , n_dry_db,  with_="steps lw 4 lc rgb 'red'")
+    plot_rd_da  = gp.PlotItems.Data(dry_edges[:-1] * 1e6 , n_dry_da,  with_="steps lw 4 lc rgb 'blue'")
+
+    g.plot(plot_rd_ub, plot_rd_ua, plot_rd_db, plot_rd_da)
+
+    # plotting
+    ymin = 0.001
+    ymax = 1000
+    xmin = 0.001
+    xmax = 500
+
+    g = gp.Gnuplot()
+    g('reset')
+    g('set term svg dynamic enhanced font "Verdana, 14"')
+    g('set output "plots/' + kernel + '_rain_ud.svg" ')
+    g('set logscale xy')
+    g('set key samplen 1.2')
+    g('set xtics rotate by 65 right (.01, .1, 1, 10, 100)')
+    g('set ytics (.001, .01, .1, 1, 10, 100, 1000)')
+    g('set xlabel "particle radius [um]"')
+    g('set ylabel "dN/dlog_{10}(r) [mg^{-1} per log_{10}(size interval)]"')
+    g('set xrange [' +  str(xmin) + ':' + str(xmax) + ']')
+    g('set yrange [' +  str(ymin) + ':' + str(ymax) + ']')
+    g('set grid')
+    g('set nokey')
+
+    plot_rw_ub  = gp.PlotItems.Data(wet_edges[:-1] * 1e6 , n_wet_ub,  with_="steps lw 4 lc rgb 'black'")
+    plot_rw_ua  = gp.PlotItems.Data(wet_edges[:-1] * 1e6 , n_wet_ua,  with_="steps lw 4 lc rgb 'green'")
+    plot_rw_db  = gp.PlotItems.Data(wet_edges[:-1] * 1e6 , n_wet_db,  with_="steps lw 4 lc rgb 'red'")
+    plot_rw_da  = gp.PlotItems.Data(wet_edges[:-1] * 1e6 , n_wet_da,  with_="steps lw 4 lc rgb 'blue'")
+    
+    g.plot(plot_rw_ub, plot_rw_ua, plot_rw_db, plot_rw_da)