Synth turb sgs model

bindings/python/lib.cpp

@@ -225,6 +225,9 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .value("implicit", lgr::as_t::implicit)
       .value("euler", lgr::as_t::euler)
       .value("pred_corr", lgr::as_t::pred_corr);
+    bp::enum_<lgr::sgs_adve_t>("sgs_adve_t") 
+      .value("ST_periodic", lgr::sgs_adve_t::ST_periodic)
+      .value("GA17", lgr::sgs_adve_t::GA17);
     bp::enum_<lgr::src_t>("src_t") 
       .value("off", lgr::src_t::off)
       .value("simple", lgr::src_t::simple)
@@ -253,7 +256,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("chem_dsc", &lgr::opts_t<real_t>::chem_dsc)
       .def_readwrite("chem_rct", &lgr::opts_t<real_t>::chem_rct)
       .def_readwrite("RH_max", &lgr::opts_t<real_t>::RH_max)
-      .def_readwrite("turb_adve", &lgr::opts_t<real_t>::turb_adve)
+      .def_readwrite("sgs_adve", &lgr::opts_t<real_t>::sgs_adve)
       .def_readwrite("turb_cond", &lgr::opts_t<real_t>::turb_cond)
       .def_readwrite("turb_coal", &lgr::opts_t<real_t>::turb_coal)
       .def_readwrite("dt", &lgr::opts_t<real_t>::dt)
@@ -292,7 +295,6 @@ 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("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)
       .def_readwrite("exact_sstp_cond", &lgr::opts_init_t<real_t>::exact_sstp_cond)
@@ -302,6 +304,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .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("sgs_adve", &lgr::opts_init_t<real_t>::sgs_adve)
       .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)
@@ -321,10 +324,17 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .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)
+      .def_readwrite("init_pair_separation", &lgr::opts_init_t<real_t>::init_pair_separation)
       .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>)
       .add_property("kernel_parameters", &lgrngn::get_kp<real_t>, &lgrngn::set_kp<real_t>)
       .def_readwrite("variable_dt_switch", &lgr::opts_init_t<real_t>::variable_dt_switch)
+      .def_readwrite("ST_eps", &lgr::opts_init_t<real_t>::ST_eps)
+      .def_readwrite("ST_Nmodes", &lgr::opts_init_t<real_t>::ST_Nmodes)
+      .def_readwrite("ST_Nwaves_max", &lgr::opts_init_t<real_t>::ST_Nwaves_max)
+      .def_readwrite("ST_Lmax", &lgr::opts_init_t<real_t>::ST_Lmax)
+      .def_readwrite("ST_Lmin", &lgr::opts_init_t<real_t>::ST_Lmin)
+      .def_readwrite("periodic_topbot_walls", &lgr::opts_init_t<real_t>::periodic_topbot_walls)
     ;
     bp::class_<lgr::particles_proto_t<real_t>/*, boost::noncopyable*/>("particles_proto_t")
       .add_property("opts_init", &lgrngn::get_oi<real_t>)
@@ -381,6 +391,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def("diag_wet_mom", &lgr::particles_proto_t<real_t>::diag_wet_mom)
       .def("diag_kappa_mom",    &lgr::particles_proto_t<real_t>::diag_kappa_mom)
       .def("diag_incloud_time_mom",    &lgr::particles_proto_t<real_t>::diag_incloud_time_mom)
+      .def("diag_pair_separation_mean",    &lgr::particles_proto_t<real_t>::diag_pair_separation_mean)
       .def("diag_wet_mass_dens", &lgr::particles_proto_t<real_t>::diag_wet_mass_dens)
       .def("diag_chem",    &lgr::particles_proto_t<real_t>::diag_chem)
       .def("diag_precip_rate",    &lgr::particles_proto_t<real_t>::diag_precip_rate)

include/libcloudph++/common/GA17_turbulence.hpp

@@ -80,7 +80,7 @@ namespace libcloudphxx
 
       template <typename real_t>
       BOOST_GPU_ENABLED
-      quantity<si::velocity, real_t> update_turb_vel(
+      quantity<si::velocity, real_t> update_sgs_vel(
         const quantity<si::velocity, real_t> &wp,
         const quantity<si::time, real_t> &tau,
         const quantity<si::time, real_t> &dt,

include/libcloudph++/lgrngn/opts.hpp

@@ -19,7 +19,7 @@ namespace libcloudphxx
     struct opts_t 
     {
       // process toggling
-      bool adve, sedi, subs, cond, coal, src, rlx, rcyc, turb_adve, turb_cond, turb_coal;
+      bool adve, sedi, subs, cond, coal, src, rlx, rcyc, sgs_adve, turb_cond, turb_coal;
 
       // RH limit for drop growth
       real_t RH_max;       
@@ -35,7 +35,7 @@ namespace libcloudphxx
       opts_t() : 
         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),
+        sgs_adve(false), turb_cond(false), turb_coal(false),
         RH_max(44), // :) (anything greater than 1.1 would be enough
         dt(-1) // negative means that we do not override dt in this step
       {

include/libcloudph++/lgrngn/opts_init.hpp

@@ -12,6 +12,7 @@
 #include "terminal_velocity.hpp"
 #include "advection_scheme.hpp"
 #include "RH_formula.hpp"
+#include "sgs_adve.hpp"
 #include "ccn_source.hpp"
 #include "../common/chem.hpp"
 
@@ -93,14 +94,25 @@ namespace libcloudphxx
            sedi_switch,  // if false no sedimentation throughout the whole simulation
            subs_switch,  // if false no subsidence 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
            exact_sstp_cond;    // if true, use per-particle sstp_cond logic, if false, use per-cell
+
+      sgs_adve_t sgs_adve; // model of SGS advection to be used. 
+      // parameters of the synthetic turbulence model (if used)
+      int ST_Nmodes,      // number of Fourier modes
+          ST_Nwaves_max;  // maximum number of wave numbers per mode
+      real_t ST_eps,      // TKE dissipation rate [m^2 / s^3], constant and uniform
+             ST_Lmax,     // maximum wave length of the ST modes [m]
+             ST_Lmin;     // minimum wave length of the ST modes [m]
 
       int sstp_chem;
       real_t chem_rho;
 
+      // initial distance between every second SD and the preceeding one [m],
+      // used in the pair separation test of the SGS advection model
+      real_t init_pair_separation;
+
       // do we want to track the time SDs spend inside clouds
       bool diag_incloud_time;
 
@@ -205,13 +217,13 @@ namespace libcloudphxx
         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(44),
         terminal_velocity(vt_t::undefined),
+        sgs_adve(sgs_adve_t::undefined),
         kernel(kernel_t::undefined),
         adve_scheme(as_t::implicit),
         RH_formula(RH_formula_t::pv_cc),
@@ -237,6 +249,12 @@ namespace libcloudphxx
         open_side_walls(false),
         periodic_topbot_walls(false),
         variable_dt_switch(false),
+        init_pair_separation(-1),
+        ST_Nmodes(10),
+        ST_Nwaves_max(10),
+        ST_eps(-1),
+        ST_Lmax(100),
+        ST_Lmin(1e-3),
         rng_seed_init_switch(false)
       {}
 

include/libcloudph++/lgrngn/particles.hpp

@@ -105,12 +105,19 @@ namespace libcloudphxx
       virtual void diag_up_mom(const int&)                                      { assert(false); }
       virtual void diag_vp_mom(const int&)                                      { assert(false); }
       virtual void diag_wp_mom(const int&)                                      { assert(false); }
+      virtual std::pair<real_t, real_t> diag_up_minmax()                        { assert(false); }
+      virtual std::pair<real_t, real_t> diag_vp_minmax()                        { assert(false); }
+      virtual std::pair<real_t, real_t> diag_wp_minmax()                        { assert(false); }
       virtual void diag_incloud_time_mom(const int&)                            { assert(false); } // requires opts_init.diag_incloud_time==true
       virtual void diag_max_rw()                                                { assert(false); }
       virtual void diag_vel_div()                                               { assert(false); }
+      virtual real_t diag_pair_separation_mean()                                { assert(false); }
+      virtual int  diag_sstp_coal()                                             { 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 remove_wet_rng(const real_t&, const real_t&)                 { assert(false); }
+
       // storing a pointer to opts_init (e.g. for interrogatin about
       // dimensions in Python bindings)
       opts_init_t<real_t> *opts_init;
@@ -187,7 +194,11 @@ namespace libcloudphxx
       void diag_up_mom(const int&);
       void diag_vp_mom(const int&);
       void diag_wp_mom(const int&);
+      std::pair<real_t, real_t> diag_up_minmax();
+      std::pair<real_t, real_t> diag_vp_minmax();
+      std::pair<real_t, real_t> diag_wp_minmax();
       void diag_incloud_time_mom(const int &k);
+      real_t diag_pair_separation_mean();
       void diag_wet_mass_dens(const real_t&, const real_t&);
 
       void diag_chem(const enum common::chem::chem_species_t&);
@@ -197,9 +208,12 @@ namespace libcloudphxx
       void diag_precip_rate();
       void diag_max_rw();
       void diag_vel_div();
+      int  diag_sstp_coal();
       std::map<libcloudphxx::common::output_t, real_t> diag_puddle();
       real_t *outbuf();
 
+      void remove_wet_rng(const real_t&, const real_t&);
+
       struct impl;
       std::unique_ptr<impl> pimpl;
 
@@ -283,18 +297,24 @@ namespace libcloudphxx
       void diag_up_mom(const int&);
       void diag_vp_mom(const int&);
       void diag_wp_mom(const int&);
+      std::pair<real_t, real_t> diag_up_minmax();
+      std::pair<real_t, real_t> diag_vp_minmax();
+      std::pair<real_t, real_t> diag_wp_minmax();
       void diag_incloud_time_mom(const int&);
+      real_t diag_pair_separation_mean();
       void diag_wet_mass_dens(const real_t&, const real_t&);
-      real_t *outbuf();
-
       void diag_chem(const enum common::chem::chem_species_t&);
       void diag_rw_ge_rc();
       void diag_RH_ge_Sc();
       void diag_all();
       void diag_precip_rate();
       void diag_max_rw();
       void diag_vel_div();
+      int  diag_sstp_coal();
       std::map<libcloudphxx::common::output_t, real_t> diag_puddle();
+      real_t *outbuf();
+
+      void remove_wet_rng(const real_t&, const real_t&);
 
       struct impl;
       std::unique_ptr<impl> pimpl;

include/libcloudph++/lgrngn/sgs_adve.hpp

@@ -0,0 +1,13 @@
+#pragma once 
+
+namespace libcloudphxx
+{
+  namespace lgrngn
+  {
+//<listing>
+    enum class sgs_adve_t { undefined, GA17, ST_periodic }; 
+//</listing>
+    // GA17 - SGS velocity modeled with a Langevin equation as in Grabowski & Abade 2017. TKE dissipation rate taken from the diss_rate array; characteristic length scale taken from the SGS_mix_len profile.
+    // ST_periodic - synthetic turbulence: SGS velocity is a sum of Fourier modes, periodic within a box. constant and uniform diss_rate; constant and uniform minimum and maximum wavelengths
+  };
+};

src/detail/kernel_definitions/hall_davis_no_waals_efficiencies_like_in_EFM.hpp

@@ -0,0 +1,110 @@
+#pragma once
+namespace libcloudphxx
+{
+  namespace lgrngn
+  {
+    namespace detail
+    {
+      template<class real_t>
+      real_t hall_davis_no_waals_r_max()
+      {
+        return 300.e-6; // not used??
+      }
+
+      const int hall_davis_no_waals_n_rad = 15;
+      const int hall_davis_no_waals_n_rat = 21;
+
+      template<class real_t>
+      void hall_davis_no_waals_efficiencies(std::vector<real_t> &vec)
+      {
+        //NVCC doesn't support extended initializer list. First make an array, then copy it to a vector.
+        const real_t arr[hall_davis_no_waals_n_rat][hall_davis_no_waals_n_rad] = 
+        {
+          {0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0},
+          {0, 0, 0, 0, 0,
+          0, 0, 0.001, 0.084, 0.05,
+          0.2543, 0.5, 0.7807, 0.87, 0.97},
+          {0, 0, 0, 0, 0,
+          0, 0, 0.07, 0.4, 0.43,
+          0.58, 0.79, 0.93, 0.96, 1},
+          {0, 0, 0, 0, 0,
+          0, 0, 0.28, 0.62, 0.64,
+          0.7629, 0.91, 0.9687, 0.98, 1},
+          {0, 0.01215, 0, 0, 0,
+          0, 0, 0.5, 0.7, 0.77,
+          0.84, 0.95, 0.95, 1, 1},
+          {0.01122, 0.01215, 0.01742, 0.0209, 0.02843,
+          0, 0.1357, 0.62, 0.79, 0.84,
+          0.8829, 0.95, 1, 1, 1},
+          {0.01122, 0.01215, 0.01742, 0.0261, 0.0341,
+          0.08008, 0.1807, 0.68, 0.83, 0.87,
+          0.9, 1, 1, 1, 1},
+          {0.01122, 0.01641, 0.02219, 0.0261, 0.03975,
+          0.09509, 0.2707, 0.74, 0.864, 0.89,
+          0.9229, 1, 1, 1, 1},
+          {0.01122, 0.01641, 0.02219, 0.03131, 0.04541,
+          0.1101, 0.3157, 0.78, 0.88, 0.9,
+          0.94, 1, 1, 1, 1},
+          {0.01503, 0.02067, 0.02696, 0.03651, 0.05107,
+          0.1251, 0.4057, 0.8, 0.9, 0.91,
+          0.95, 1, 1, 1, 1},
+          {0.01503, 0.02067, 0.02696, 0.03822, 0.05673,
+          0.143, 0.4507, 0.8, 0.9, 0.91,
+          0.95, 1, 1, 1, 1},
+          {0.01503, 0.02067, 0.02696, 0.03822, 0.05381,
+          0.1336, 0.39, 0.8, 0.9, 0.91,
+          0.95, 1, 1, 1, 1},
+          {0.02, 0.02264, 0.02696, 0.03644, 0.05089,
+          0.1242, 0.3597, 0.78, 0.9, 0.91,
+          0.95, 1, 1, 1, 1},
+          {0.02, 0.02264, 0.02696, 0.03466, 0.04797,
+          0.1148, 0.299, 0.77, 0.888, 0.91,
+          0.95, 1, 1, 1, 1},
+          {0.02, 0.0246, 0.02696, 0.03288, 0.04504,
+          0.09603, 0.2687, 0.76, 0.88, 0.92,
+          0.95, 1, 1, 1, 1},
+          {0.02497, 0.0246, 0.02696, 0.03288, 0.04212,
+          0.08665, 0.2081, 0.77, 0.882, 0.93,
+          0.9743, 1, 1, 1, 1},
+          {0.02497, 0.0246, 0.02696, 0.0311, 0.0392,
+          0.07726, 0.1778, 0.77, 0.89, 0.95,
+          1, 1, 1, 1, 1},
+          {0.02497, 0.02657, 0.02696, 0.02932, 0.03627,
+          0.06788, 0.1171, 0.78, 0.938, 1,
+          1.023, 1, 1, 1, 1},
+          {0.02497, 0.02657, 0.02696, 0.02843, 0.03335,
+          0.0538, 0.08679, 0.79, 1.01, 1.03,
+          1.04, 1, 1, 1, 1},
+          {0.02696, 0.02696, 0.02696, 0.02843, 0.03335,
+          0.0538, 0.08679, 0.95, 1.5, 1.7,
+          2.543, 1, 1, 1, 1},
+          {0.02696, 0.02696, 0.02696, 0.02843, 0.03335,
+          0.0538, 0.08679, 1.4, 2.3, 3,
+          4, 1, 1, 1, 1}
+        };
+
+        vec = std::vector<real_t>(&(arr[0][0]), &(arr[hall_davis_no_waals_n_rat][hall_davis_no_waals_n_rad]));
+      }
+
+      template<class real_t>
+      void hall_davis_no_waals_radii(std::vector<real_t> &vec)
+      {
+        //NVCC doesn't support extended initializer list. First make an array, then copy it to a vector.
+        const real_t arr[hall_davis_no_waals_n_rad] = {6.e-6,8.e-6,10.e-6,15.e-6,20.e-6,25.e-6,30.e-6,40.e-6,50.e-6,60.e-6,70.e-6,100.e-6,150.e-6,200.e-6,300.e-6}; 
+
+        vec = std::vector<real_t>(arr, arr + hall_davis_no_waals_n_rad );
+      }
+
+      template<class real_t>
+      void hall_davis_no_waals_ratios(std::vector<real_t> &vec)
+      {
+        //NVCC doesn't support extended initializer list. First make an array, then copy it to a vector.
+        const real_t arr[hall_davis_no_waals_n_rat] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5, 0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0};
+
+        vec = std::vector<real_t>(arr, arr + hall_davis_no_waals_n_rat );
+      }
+    }
+  }
+}