Merge pull request #427 from pdziekan/relax_ccn_profile_and_accr_acnv_diag

Relax ccn profile and accr acnv diag

Author: pdziekan

bindings/python/lgrngn.hpp

@@ -319,6 +319,35 @@ namespace libcloudphxx
         }
       }
 
+      template <typename real_t>
+      void set_rdd( // rlx_dry_distros
+        lgr::opts_init_t<real_t> *arg,
+        const bp::dict &kappa_func
+      )
+      {
+        arg->rlx_dry_distros.clear();
+        if(len(kappa_func.keys()) == 0)
+          return;
+
+        // loop over kappas
+        for (int i = 0; i < len(kappa_func.keys()); ++i)
+        {
+          std::cerr << i << std::endl;
+          const real_t kappa = bp::extract<real_t>(kappa_func.keys()[i]);
+          std::cerr << kappa << std::endl;
+          const bp::list nlnrd_kparange_zrange_list = bp::extract<bp::list>(kappa_func.values()[i]);
+          assert(len(nlnrd_kparange_zrange_list) == 3);
+          auto nlnrd = std::make_shared<detail::pyunary<real_t>>(nlnrd_kparange_zrange_list[0]);
+          const bp::list kparange_list = bp::extract<bp::list>(nlnrd_kparange_zrange_list[1]);
+          std::pair<real_t, real_t> kparange{bp::extract<real_t>(kparange_list[0]), bp::extract<real_t>(kparange_list[1])};
+          std::cerr << kparange.first << " " << kparange.second << std::endl;
+          const bp::list zrange_list = bp::extract<bp::list>(nlnrd_kparange_zrange_list[2]);
+          std::pair<real_t, real_t> zrange{bp::extract<real_t>(zrange_list[0]), bp::extract<real_t>(zrange_list[1])};
+          std::cerr << zrange.first << " " << zrange.second << std::endl;
+          arg->rlx_dry_distros[kappa] = std::make_tuple(nlnrd, kparange, zrange);
+        }
+      }
+
       template <typename real_t>
       void get_ds(
         lgr::opts_init_t<real_t> *arg
@@ -351,6 +380,14 @@ namespace libcloudphxx
         throw std::runtime_error("source_dry_distros does not feature a getter yet - TODO");
       }
 
+      template <typename real_t>
+      void get_rdd(
+        lgr::opts_init_t<real_t> *arg
+      )
+      {
+        throw std::runtime_error("relax_dry_distros does not feature a getter yet - TODO");
+      }
+
       template <typename real_t>
       void set_kp(
         lgr::opts_init_t<real_t> *arg,

bindings/python/lib.cpp

@@ -243,6 +243,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("coal", &lgr::opts_t<real_t>::coal)
       .def_readwrite("src", &lgr::opts_t<real_t>::src)
       .def_readwrite("rcyc", &lgr::opts_t<real_t>::rcyc)
+      .def_readwrite("rlx", &lgr::opts_t<real_t>::rlx)
       .def_readwrite("chem_dsl", &lgr::opts_t<real_t>::chem_dsl)
       .def_readwrite("chem_dsc", &lgr::opts_t<real_t>::chem_dsc)
       .def_readwrite("chem_rct", &lgr::opts_t<real_t>::chem_rct)
@@ -257,6 +258,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .add_property("dry_sizes", &lgrngn::get_ds<real_t>, &lgrngn::set_ds<real_t>)
       .add_property("src_dry_distros", &lgrngn::get_sdd<real_t>, &lgrngn::set_sdd<real_t>)
       .add_property("src_dry_sizes", &lgrngn::get_ds<real_t>, &lgrngn::set_sds<real_t>)
+      .add_property("rlx_dry_distros", &lgrngn::get_rdd<real_t>, &lgrngn::set_rdd<real_t>)
       .def_readwrite("nx", &lgr::opts_init_t<real_t>::nx)
       .def_readwrite("ny", &lgr::opts_init_t<real_t>::ny)
       .def_readwrite("nz", &lgr::opts_init_t<real_t>::nz)
@@ -285,6 +287,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("sedi_switch", &lgr::opts_init_t<real_t>::sedi_switch)
       .def_readwrite("subs_switch", &lgr::opts_init_t<real_t>::subs_switch)
       .def_readwrite("src_switch", &lgr::opts_init_t<real_t>::src_switch)
+      .def_readwrite("rlx_switch", &lgr::opts_init_t<real_t>::rlx_switch)
       .def_readwrite("turb_adve_switch", &lgr::opts_init_t<real_t>::turb_adve_switch)
       .def_readwrite("turb_cond_switch", &lgr::opts_init_t<real_t>::turb_cond_switch)
       .def_readwrite("turb_coal_switch", &lgr::opts_init_t<real_t>::turb_coal_switch)
@@ -293,20 +296,25 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("sstp_coal", &lgr::opts_init_t<real_t>::sstp_coal)
       .def_readwrite("sstp_chem", &lgr::opts_init_t<real_t>::sstp_chem)
       .def_readwrite("supstp_src", &lgr::opts_init_t<real_t>::supstp_src)
+      .def_readwrite("supstp_rlx", &lgr::opts_init_t<real_t>::supstp_rlx)
       .def_readwrite("kernel", &lgr::opts_init_t<real_t>::kernel)
       .def_readwrite("adve_scheme", &lgr::opts_init_t<real_t>::adve_scheme)
       .def_readwrite("sd_conc", &lgr::opts_init_t<real_t>::sd_conc)
       .def_readwrite("sd_conc_large_tail", &lgr::opts_init_t<real_t>::sd_conc_large_tail)
       .def_readwrite("aerosol_independent_of_rhod", &lgr::opts_init_t<real_t>::aerosol_independent_of_rhod)
       .def_readwrite("sd_const_multi", &lgr::opts_init_t<real_t>::sd_const_multi)
       .def_readwrite("src_sd_conc", &lgr::opts_init_t<real_t>::src_sd_conc)
+      .def_readwrite("rlx_bins", &lgr::opts_init_t<real_t>::rlx_bins)
+      .def_readwrite("rlx_sd_per_bin", &lgr::opts_init_t<real_t>::rlx_sd_per_bin)
+      .def_readwrite("rlx_timescale", &lgr::opts_init_t<real_t>::rlx_timescale)
       .def_readwrite("n_sd_max", &lgr::opts_init_t<real_t>::n_sd_max)
       .def_readwrite("terminal_velocity", &lgr::opts_init_t<real_t>::terminal_velocity)
       .def_readwrite("RH_formula", &lgr::opts_init_t<real_t>::RH_formula)
       .def_readwrite("chem_rho", &lgr::opts_init_t<real_t>::chem_rho)
       .def_readwrite("RH_max", &lgr::opts_init_t<real_t>::RH_max)
       .def_readwrite("rng_seed", &lgr::opts_init_t<real_t>::rng_seed)
       .def_readwrite("rng_seed_init", &lgr::opts_init_t<real_t>::rng_seed_init)
+      .def_readwrite("rng_seed_init_switch", &lgr::opts_init_t<real_t>::rng_seed_init_switch)
       .def_readwrite("diag_incloud_time", &lgr::opts_init_t<real_t>::diag_incloud_time)
       .add_property("w_LS", &lgrngn::get_w_LS<real_t>, &lgrngn::set_w_LS<real_t>)
       .add_property("SGS_mix_len", &lgrngn::get_SGS_mix_len<real_t>, &lgrngn::set_SGS_mix_len<real_t>)

include/libcloudph++/blk_1m/extincl.hpp

@@ -2,7 +2,10 @@
 
 #include <algorithm>
 
-#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+#include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
+#include <boost/numeric/odeint/algebra/default_operations.hpp>
+#include <boost/numeric/odeint/util/resizer.hpp>
 
 #include "../common/const_cp.hpp"
 #include "../common/theta_dry.hpp"

include/libcloudph++/lgrngn/opts.hpp

@@ -19,7 +19,7 @@ namespace libcloudphxx
     struct opts_t 
     {
       // process toggling
-      bool adve, sedi, subs, cond, coal, src, rcyc, turb_adve, turb_cond, turb_coal;
+      bool adve, sedi, subs, cond, coal, src, rlx, rcyc, turb_adve, turb_cond, turb_coal;
 
       // RH limit for drop growth
       real_t RH_max;       
@@ -33,7 +33,7 @@ namespace libcloudphxx
 
       // ctor with defaults (C++03 compliant) ...
       opts_t() : 
-        adve(true), sedi(true), subs(false), cond(true), coal(true), src(false), rcyc(false),
+        adve(true), sedi(true), subs(false), cond(true), coal(true), src(false), rlx(false), rcyc(false),
         chem_dsl(false), chem_dsc(false), chem_rct(false),
         turb_adve(false), turb_cond(false), turb_coal(false),
         RH_max(44), // :) (anything greater than 1.1 would be enough

include/libcloudph++/lgrngn/opts_init.hpp

@@ -48,9 +48,6 @@ namespace libcloudphxx
 
       // no. of substeps 
       int sstp_cond, sstp_coal; 
-  
-      // timestep interval at which source will be applied
-      int supstp_src;
 
       // Lagrangian domain extents
       real_t x0, y0, z0, x1, y1, z1;
@@ -74,19 +71,6 @@ namespace libcloudphxx
       // should be enough to store particles from sources
       unsigned long long n_sd_max; 
 
-      // source distro per unit time
-      dry_distros_t src_dry_distros;
-
-      // number of SDs created from src_dry_distros per cell per source iteration
-      unsigned long long src_sd_conc;
-
-      // dry sizes of droplets added from the source, STP_concentration created per unit time instead of the STP_concentration 
-      dry_sizes_t src_dry_sizes;
-
-      // box in which aerosol from source will be created
-      // will be rounded to cell number - cells are supposed to be uniform
-      real_t src_x0, src_y0, src_z0, src_x1, src_y1, src_z1;
-
       // coalescence Kernel type
       kernel_t::kernel_t kernel;
 
@@ -108,6 +92,7 @@ namespace libcloudphxx
            sedi_switch,  // if false no sedimentation throughout the whole simulation
            subs_switch,  // if false no subsidence throughout the whole simulation
            src_switch,   // if false no source throughout the whole simulation
+           rlx_switch,   // if false no relaxation throughout the whole simulation
            turb_adve_switch,   // if true, turbulent motion of SDs is modeled
            turb_cond_switch,   // if true, turbulent condensation of SDs is modeled
            turb_coal_switch,   // if true, turbulent coalescence kernels can be used
@@ -126,6 +111,9 @@ namespace libcloudphxx
       int rng_seed,
           rng_seed_init; // seed used to init SD (positions and dry sizes)
 
+      // should the separate rng seed for initialization be used?
+      bool rng_seed_init_switch;
+
       // no of GPUs per MPI node to use, 0 for all available
       int dev_count; 
 
@@ -146,6 +134,54 @@ namespace libcloudphxx
            periodic_topbot_walls; // if true, top and bot walls are periodic. Open otherwise
 
 
+      // --- aerosol source stuff ---
+
+      // source distro per unit time
+      dry_distros_t src_dry_distros;
+
+      // number of SDs created from src_dry_distros per cell per source iteration
+      unsigned long long src_sd_conc;
+
+      // dry sizes of droplets added from the source, STP_concentration created per unit time instead of the STP_concentration 
+      dry_sizes_t src_dry_sizes;
+
+      // box in which aerosol from source will be created
+      // will be rounded to cell number - cells are supposed to be uniform
+      real_t src_x0, src_y0, src_z0, src_x1, src_y1, src_z1;
+  
+      // timestep interval at which source will be applied
+      int supstp_src;
+
+
+      // --- aerosol relaxation stuff ---
+      // initial dry sizes of aerosol
+      // defined with a distribution
+      // uses shared_ptr to make opts_init copyable
+      typedef std::unordered_map<
+        real_t,                // kappa
+        std::tuple<
+          std::shared_ptr<unary_function<real_t>>, // n(ln(rd)) @ STP; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true
+          std::pair<real_t, real_t>, // kappa range of CCN considered to belong to this distribution, ranges of different members of the map need to be exclusive (TODO: add a check of this)
+          std::pair<real_t, real_t>  // range of altitudes at which this relaxation acts
+        >
+      > rlx_dry_distros_t;
+
+      rlx_dry_distros_t rlx_dry_distros;
+
+      // number of bins into which the relaxation distro is divided
+      unsigned long long rlx_bins;
+
+      // number of SD created per bin
+      real_t rlx_sd_per_bin; // floating, because it gets divided by the number of GPUs per node * number of nodes
+  
+      // timestep interval at which relaxation will be applied
+      int supstp_rlx;
+
+      // relaxation time scale [s]
+      real_t rlx_timescale;
+
+      // -- ctors ---
+
       // ctor with defaults (C++03 compliant) ...
       opts_init_t() : 
         nx(0), ny(0), nz(0),
@@ -158,20 +194,20 @@ namespace libcloudphxx
         sd_const_multi(0),
         dt(0),   
         sstp_cond(1), sstp_coal(1), sstp_chem(1),         
-        supstp_src(1),
         chem_switch(false),  // chemical reactions turned off by default
         sedi_switch(true),  // sedimentation turned on by default
         subs_switch(false),  // subsidence turned off by default
         coal_switch(true),  // coalescence turned on by default
         src_switch(false),  // source turned off by default
+        rlx_switch(false), 
         exact_sstp_cond(false),
         turb_cond_switch(false),
         turb_adve_switch(false),
         turb_coal_switch(false),
         RH_max(.95), // value seggested in Lebo and Seinfeld 2011
         chem_rho(0), // dry particle density  //TODO add checking if the user gave a different value (np w init)  (was 1.8e-3)
         rng_seed(44),
-        rng_seed_init(rng_seed),
+        rng_seed_init(44),
         terminal_velocity(vt_t::undefined),
         kernel(kernel_t::undefined),
         adve_scheme(as_t::implicit),
@@ -186,12 +222,18 @@ namespace libcloudphxx
         src_y1(0),
         src_z0(0),
         src_z1(0),
+        supstp_src(1),
+        rlx_bins(0),
+        rlx_timescale(1),
+        rlx_sd_per_bin(0),
+        supstp_rlx(1),
         rd_min(0.),
         diag_incloud_time(false),
         no_ccn_at_init(false),
         open_side_walls(false),
         periodic_topbot_walls(false),
-        variable_dt_switch(false)
+        variable_dt_switch(false),
+        rng_seed_init_switch(false)
       {}
 
       // dtor (just to silence -Winline warnings)

include/libcloudph++/lgrngn/particles.hpp

@@ -110,6 +110,8 @@ namespace libcloudphxx
       virtual void diag_vel_div()                                               { assert(false); }
       virtual std::map<libcloudphxx::common::output_t, real_t> diag_puddle()    { assert(false); return std::map<libcloudphxx::common::output_t, real_t>(); }
       virtual real_t *outbuf()                                                  { assert(false); return NULL; }
+      virtual void diag_accr25()                                   { assert(false); }
+      virtual void diag_acnv25()                                   { assert(false); }
 
       // storing a pointer to opts_init (e.g. for interrogatin about
       // dimensions in Python bindings)
@@ -199,6 +201,8 @@ namespace libcloudphxx
       void diag_vel_div();
       std::map<libcloudphxx::common::output_t, real_t> diag_puddle();
       real_t *outbuf();
+      virtual void diag_accr25();
+      virtual void diag_acnv25();
 
       struct impl;
       std::unique_ptr<impl> pimpl;
@@ -286,6 +290,8 @@ namespace libcloudphxx
       void diag_incloud_time_mom(const int&);
       void diag_wet_mass_dens(const real_t&, const real_t&);
       real_t *outbuf();
+      virtual void diag_accr25();
+      virtual void diag_acnv25();
 
       void diag_chem(const enum common::chem::chem_species_t&);
       void diag_rw_ge_rc();

src/detail/config.hpp

@@ -30,6 +30,8 @@ namespace libcloudphxx
 
         const real_t bcond_tolerance = 2e-4; // [m]; error tolerance for position near bcond after distmem copy  
 
+        const real_t rlx_conc_tolerance = 0.1; // tolerance of the relaxation scheme; new SD will be created if missing_conc/expected_conc > tolerance
+
         // ctor
         config():
           vt0_ln_r_min(log(5e-7)),

src/detail/distmem_opts.hpp

@@ -46,6 +46,8 @@ namespace libcloudphxx
 
         opts_init.n_sd_max = opts_init.n_sd_max / size + 1;
 
+        opts_init.rlx_sd_per_bin /= size;
+
         return n_x_bfr;
       }
     }

src/detail/urand.hpp

@@ -66,7 +66,7 @@ namespace libcloudphxx
           const thrust_size_t n
         ) {
           // note: generate_n copies the third argument!!!
-          std::generate_n(u01.begin(), n, fnctr_u01({engine, dist_u01})); 
+          std::generate_n(u01.begin(), n, fnctr_u01({engine, dist_u01})); // [0,1) range 
         }
 
         void generate_normal_n(
@@ -131,20 +131,25 @@ namespace libcloudphxx
           const thrust_size_t n
         )
         {
-          int status = curandGenerateUniform(gen, thrust::raw_pointer_cast(v.data()), n);
+          int status = curandGenerateUniform(gen, thrust::raw_pointer_cast(v.data()), n); // (0,1] range
           assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/);
           _unused(status);
-
+          // shift into the expected [0,1) range
+          namespace arg = thrust::placeholders;
+          thrust::transform(v.begin(), v.begin() + n, v.begin(), float(1) - arg::_1);
         }
 
         void generate_n(
           thrust_device::vector<double> &v, 
           const thrust_size_t n
         )
         {
-          int status = curandGenerateUniformDouble(gen, thrust::raw_pointer_cast(v.data()), n);
+          int status = curandGenerateUniformDouble(gen, thrust::raw_pointer_cast(v.data()), n); // (0,1] range
           assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/);
           _unused(status);
+          // shift into the expected [0,1) range
+          namespace arg = thrust::placeholders;
+          thrust::transform(v.begin(), v.begin() + n, v.begin(), double(1) - arg::_1);
         }
 
         void generate_normal_n(