Merge pull request #366 from pdziekan/turbulent_coalescence

Turbulent coalescence kernels: use dissipation rate from sync_in (includes turb adve and cond)

CMakeLists.txt

@@ -262,6 +262,8 @@ if (EXISTS "${CMAKE_SOURCE_DIR}/.git")
   )
 endif()
 
+include_directories(${libcloudphxx_INCLUDE_DIRS})
+
 add_subdirectory(src) 
 enable_testing()
 add_subdirectory(tests)

bindings/python/lgrngn.hpp

@@ -113,6 +113,7 @@ namespace libcloudphxx
 	const bp_array &Cx,
 	const bp_array &Cy,
 	const bp_array &Cz,
+	const bp_array &diss_rate,
         bp::dict &ambient_chem
       )
       {
@@ -136,6 +137,7 @@ namespace libcloudphxx
 	  np2ai<real_t>(Cx, sz(*arg)),
 	  np2ai<real_t>(Cy, sz(*arg)),
 	  np2ai<real_t>(Cz, sz(*arg)),
+	  np2ai<real_t>(diss_rate, sz(*arg)),
           map
 	);
       }
@@ -150,6 +152,7 @@ namespace libcloudphxx
 	const bp_array &Cx,
 	const bp_array &Cy,
 	const bp_array &Cz,
+	const bp_array &diss_rate,
         bp::dict &ambient_chem
       )
       {
@@ -172,6 +175,7 @@ namespace libcloudphxx
 	  np2ai<real_t>(Cx, sz(*arg)),
 	  np2ai<real_t>(Cy, sz(*arg)),
 	  np2ai<real_t>(Cz, sz(*arg)),
+	  np2ai<real_t>(diss_rate, sz(*arg)),
           map
 	);
       }

bindings/python/lib.cpp

@@ -244,6 +244,9 @@ 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("turb_cond", &lgr::opts_t<real_t>::turb_cond)
+      .def_readwrite("turb_coal", &lgr::opts_t<real_t>::turb_coal)
     ;
     bp::class_<lgr::opts_init_t<real_t>>("opts_init_t")
       .add_property("dry_distros", &lgrngn::get_dd<real_t>, &lgrngn::set_dd<real_t>)
@@ -269,6 +272,9 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("coal_switch", &lgr::opts_init_t<real_t>::coal_switch)
       .def_readwrite("sedi_switch", &lgr::opts_init_t<real_t>::sedi_switch)
       .def_readwrite("src_switch", &lgr::opts_init_t<real_t>::src_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)
       .def_readwrite("sstp_cond", &lgr::opts_init_t<real_t>::sstp_cond)
       .def_readwrite("sstp_coal", &lgr::opts_init_t<real_t>::sstp_coal)
@@ -308,6 +314,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
         bp::arg("Cx")  = BP_ARR_FROM_BP_OBJ,
         bp::arg("Cy")  = BP_ARR_FROM_BP_OBJ,
         bp::arg("Cz")  = BP_ARR_FROM_BP_OBJ,
+        bp::arg("diss_rate")  = BP_ARR_FROM_BP_OBJ,
         bp::arg("ambient_chem") = bp::dict()
       ))
       .def("sync_in",    &lgrngn::sync_in<real_t>, (
@@ -317,6 +324,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
         bp::arg("Cx")  = BP_ARR_FROM_BP_OBJ,
         bp::arg("Cy")  = BP_ARR_FROM_BP_OBJ,
         bp::arg("Cz")  = BP_ARR_FROM_BP_OBJ,
+        bp::arg("diss_rate")  = BP_ARR_FROM_BP_OBJ,
         bp::arg("ambient_chem") = bp::dict()
       ))
       .def("step_cond",    &lgrngn::step_cond<real_t>, (

include/libcloudph++/common/turbulence.hpp

@@ -0,0 +1,136 @@
+#pragma once
+
+namespace libcloudphxx
+{
+  namespace common
+  {
+    // Grabowski and Abade 2017
+    namespace turbulence
+    {
+      libcloudphxx_const_derived(si::dimensionless, C_E, real_t(0.845));
+      libcloudphxx_const_derived(si::dimensionless, C_tau, real_t(1.5));
+      libcloudphxx_const_derived(si::dimensionless, cube_root_of_two_pi, pow(real_t(2) *
+#if !defined(__NVCC__)
+        pi<real_t>()
+#else
+        CUDART_PI
+#endif
+        , real_t(1./3.))
+      );
+
+#if !defined(__NVCC__)
+      using std::pow;
+      using std::sqrt;
+#endif
+
+      typedef divide_typeof_helper<
+        si::dimensionless,  
+        si::length
+      >::type one_over_length;
+      libcloudphxx_const(one_over_length, a_1, real_t(3e-4), real_t(1) / si::meters);
+
+      typedef divide_typeof_helper<
+        si::area,  
+        si::time
+       >::type area_over_time;
+      libcloudphxx_const(area_over_time, a_2, real_t(2.8e-4), si::square_meters / si::seconds);
+
+      typedef divide_typeof_helper<
+        si::area,  
+        multiply_typeof_helper<si::time, si::time>::type
+      >::type area_over_time_squared;
+
+      typedef divide_typeof_helper<
+        area_over_time_squared, 
+        si::time
+      >::type area_over_time_cubed;
+
+      typedef divide_typeof_helper<
+        si::dimensionless, 
+        si::time
+      >::type one_over_time;
+
+      typedef divide_typeof_helper<
+        si::dimensionless, 
+        si::area
+      >::type one_over_area;
+
+      template <typename real_t>
+      quantity<si::length, real_t> length_scale(
+        quantity<si::length, real_t> dx
+      )
+      { return quantity<si::length, real_t>(dx);}
+
+      template <typename real_t>
+      quantity<si::length, real_t> length_scale(
+        quantity<si::length, real_t> dx,
+        quantity<si::length, real_t> dz
+      )
+      //{ return quantity<si::length, real_t>(sqrt(dx*dz));}
+      { return quantity<si::length, real_t>(dz);}
+
+      template <typename real_t>
+      quantity<si::length, real_t> length_scale(
+        quantity<si::length, real_t> dx,
+        quantity<si::length, real_t> dy,
+        quantity<si::length, real_t> dz
+      )
+      //{ return quantity<si::length, real_t>(pow((dx/si::metres)*(dy/si::metres)*(dz/si::metres), real_t(1./3.)) * si::metres);}
+      { return quantity<si::length, real_t>(dz);}
+
+      template <typename real_t>
+      BOOST_GPU_ENABLED
+      quantity<area_over_time_squared, real_t> tke(
+        const quantity<area_over_time_cubed, real_t> &diss_rate, // dissipation rate (epsilon)
+        const quantity<si::length, real_t> &L // characteristic length-scale
+      )
+      {
+        return quantity<area_over_time_squared, real_t>(pow((L * diss_rate) / si::cubic_metres * si::seconds * si::seconds * si::seconds / C_E<real_t>(), real_t(2./3.)) * si::metres * si::metres / si::seconds / si::seconds);
+      };
+
+      template <typename real_t>
+      BOOST_GPU_ENABLED
+      quantity<si::time, real_t> tau(
+        const quantity<area_over_time_squared, real_t> &tke,
+        const quantity<si::length, real_t> &L // characteristic length-scale
+      )
+      {
+        return quantity<si::time, real_t>(L / cube_root_of_two_pi<real_t>()  * sqrt(C_tau<real_t>() / tke));
+      };
+
+      template <typename real_t>
+      BOOST_GPU_ENABLED
+      quantity<si::velocity, real_t> update_turb_vel(
+        const quantity<si::velocity, real_t> &wp,
+        const quantity<si::time, real_t> &tau,
+        const quantity<si::time, real_t> &dt,
+        const quantity<area_over_time_squared, real_t> &tke,
+        const quantity<si::dimensionless, real_t> &r_normal
+      )
+      {
+        quantity<si::dimensionless, real_t> exp_m_dt_ov_tau(exp(-dt / tau));
+        return quantity<si::velocity, real_t>(wp * exp_m_dt_ov_tau + sqrt((real_t(1) - exp_m_dt_ov_tau * exp_m_dt_ov_tau) * real_t(2./3.) * tke ) * r_normal);
+      };
+
+      template <typename real_t>
+      BOOST_GPU_ENABLED
+      quantity<si::time, real_t> tau_relax(
+        const quantity<one_over_area, real_t> &wet_mom_1_over_vol
+      )
+      {
+        return quantity<si::time, real_t>(real_t(1) / ( a_2<real_t>() * wet_mom_1_over_vol ) ); 
+      };
+
+      template <typename real_t>
+      BOOST_GPU_ENABLED
+      quantity<one_over_time, real_t> dot_turb_ss(
+        const quantity<si::dimensionless, real_t> &ssp,
+        const quantity<si::velocity, real_t> &wp,
+        const quantity<si::time, real_t> &tau_rlx
+      )
+      {
+        return quantity<one_over_time, real_t>(a_1<real_t>() *  wp - ssp / tau_rlx);
+      };
+    };
+  };
+};

include/libcloudph++/lgrngn/opts.hpp

@@ -19,7 +19,7 @@ namespace libcloudphxx
     struct opts_t 
     {
       // process toggling
-      bool adve, sedi, cond, coal, src, rcyc;
+      bool adve, sedi, cond, coal, src, rcyc, turb_adve, turb_cond, turb_coal;
 
       // RH limit for drop growth
       real_t RH_max;       
@@ -32,6 +32,7 @@ namespace libcloudphxx
       opts_t() : 
         adve(true), sedi(true), cond(true), coal(true), src(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

@@ -102,7 +102,10 @@ namespace libcloudphxx
            coal_switch,  // if false no coalescence throughout the whole simulation
            sedi_switch,  // if false no sedimentation throughout the whole simulation
            src_switch,   // if false no source throughout the whole simulation
-           exact_sstp_cond; // if true, use per-particle sstp_cond logic, if false, use per-cell
+           exact_sstp_cond, // if true, use per-particle sstp_cond logic, if false, use per-cell
+           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
 
       int sstp_chem;
       real_t chem_rho;
@@ -122,6 +125,8 @@ namespace libcloudphxx
       // large-scale horizontal wind divergence [1/s], used to calculate subsidence rate as -div_LS*z
       real_t div_LS;
 
+      real_t rd_min; // minimal dry radius of droplets (works only for init from spectrum)
+
       // ctor with defaults (C++03 compliant) ...
       opts_init_t() : 
         nx(0), ny(0), nz(0),
@@ -140,6 +145,9 @@ namespace libcloudphxx
         coal_switch(true),  // coalescence turned on by default
         src_switch(false),  // source turned off by default
         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),
@@ -152,7 +160,8 @@ namespace libcloudphxx
         n_sd_max(0),
         src_sd_conc(0),
         src_z1(0),
-        div_LS(0.)
+        div_LS(0.),
+        rd_min(0.)
       {}
 
       // dtor (just to silence -Winline warnings)

include/libcloudph++/lgrngn/particles.hpp

@@ -39,6 +39,7 @@ namespace libcloudphxx
         const arrinfo_t<real_t> courant_x = arrinfo_t<real_t>(),
         const arrinfo_t<real_t> courant_y = arrinfo_t<real_t>(),
         const arrinfo_t<real_t> courant_z = arrinfo_t<real_t>(),
+        const arrinfo_t<real_t> diss_rate = arrinfo_t<real_t>(), // TKE dissipation rate (epsilon)
         std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem = std::map<enum chem_species_t, arrinfo_t<real_t> >()
       ) { 
         assert(false); 
@@ -51,6 +52,7 @@ namespace libcloudphxx
         const arrinfo_t<real_t> courant_x = arrinfo_t<real_t>(),
         const arrinfo_t<real_t> courant_y = arrinfo_t<real_t>(),
         const arrinfo_t<real_t> courant_z = arrinfo_t<real_t>(),
+        const arrinfo_t<real_t> diss_rate = arrinfo_t<real_t>(),
         std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem = std::map<enum chem_species_t, arrinfo_t<real_t> >()
       ) { 
         assert(false); 
@@ -127,6 +129,7 @@ namespace libcloudphxx
         const arrinfo_t<real_t> courant_x,
         const arrinfo_t<real_t> courant_y,
         const arrinfo_t<real_t> courant_z,
+        const arrinfo_t<real_t> diss_rate,
         std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem
       );
 
@@ -137,6 +140,7 @@ namespace libcloudphxx
         const arrinfo_t<real_t> courant_x,
         const arrinfo_t<real_t> courant_y,
         const arrinfo_t<real_t> courant_z,
+        const arrinfo_t<real_t> diss_rate,
         std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem
       );
 
@@ -221,6 +225,7 @@ namespace libcloudphxx
         const arrinfo_t<real_t> courant_x,
         const arrinfo_t<real_t> courant_y,
         const arrinfo_t<real_t> courant_z,
+        const arrinfo_t<real_t> diss_rate,
         std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem 
       );
 
@@ -231,6 +236,7 @@ namespace libcloudphxx
         const arrinfo_t<real_t> courant_x,
         const arrinfo_t<real_t> courant_y,
         const arrinfo_t<real_t> courant_z,
+        const arrinfo_t<real_t> diss_rate,
         std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem 
       );
 

models/kinematic_2D/src/kin_cloud_2d_lgrngn_chem.hpp

@@ -259,6 +259,7 @@ class kin_cloud_2d_lgrngn_chem : public kin_cloud_2d_lgrngn<ct_params_t>
         libcloudphxx::lgrngn::arrinfo_t<real_t>(),
         libcloudphxx::lgrngn::arrinfo_t<real_t>(),
         libcloudphxx::lgrngn::arrinfo_t<real_t>(),
+        libcloudphxx::lgrngn::arrinfo_t<real_t>(),
         ambient_chem
       ); 
 

src/detail/kernel_onishi_nograv.hpp

@@ -29,6 +29,8 @@ namespace libcloudphxx
       BOOST_GPU_ENABLED
       real_t kernel_onishi_nograv(const real_t &r1, const real_t &r2, const real_t &Re_l, const real_t &eps, real_t dnu, real_t ratio_den)
       {
+        if(eps < 1e-10) return 0.;
+
 #if !defined(__NVCC__)
         using std::max;
         using std::pow;

src/detail/kernels.hpp

@@ -202,14 +202,14 @@ namespace libcloudphxx
 
       //ctor
       BOOST_GPU_ENABLED
-      kernel_onishi(thrust_device::pointer<real_t> k_params, real_t r_max) : kernel_geometric<real_t, n_t>(k_params, 2, r_max) {}
+      kernel_onishi(thrust_device::pointer<real_t> k_params, real_t r_max) : kernel_geometric<real_t, n_t>(k_params, 1, r_max) {}
       kernel_onishi() = default;
 
       BOOST_GPU_ENABLED
       virtual real_t calc(const tpl_calc_wrap<real_t,n_t> &tpl_wrap) const
       {
         enum { n_a_ix, n_b_ix, rw2_a_ix, rw2_b_ix, vt_a_ix, vt_b_ix, rd3_a_ix, rd3_b_ix };
-        enum { rhod_ix, eta_ix };
+        enum { rhod_ix, eta_ix, diss_rate_ix };
 
 #if !defined(__NVCC__)
         using std::sqrt;
@@ -218,7 +218,7 @@ namespace libcloudphxx
         real_t rwb = sqrt( thrust::get<rw2_b_ix>(tpl_wrap.get_rw()));
         real_t onishi_nograv = detail::kernel_onishi_nograv<real_t>           // value of the turbulent onishi kernel that neglects gravitational settling
         (
-          rwa, rwb, kernel_base<real_t, n_t>::k_params[1], kernel_base<real_t, n_t>::k_params[0],                              // k_params[0] - epsilon, k_params[1] - Re_lambda
+          rwa, rwb, kernel_base<real_t, n_t>::k_params[0], thrust::get<diss_rate_ix>(tpl_wrap.get_ro_calc()),                  // k_params[0] - Re_lambda
           thrust::get<eta_ix>(tpl_wrap.get_ro_calc()) / thrust::get<rhod_ix>(tpl_wrap.get_ro_calc()),                          // kinetic viscosity 
           common::moist_air::rho_w<real_t>() / si::kilograms * si::cubic_metres / thrust::get<rhod_ix>(tpl_wrap.get_ro_calc()) // ratio of water to air density
         );

src/detail/tpl_calc_wrapper.hpp

@@ -1,6 +1,8 @@
 #pragma once
 #include <thrust/tuple.h>
 
+// TODO: enums used to get from these tuples are defined separately in at least two other places, define them once here
+
 namespace libcloudphxx
 {
   namespace lgrngn
@@ -21,7 +23,8 @@ namespace libcloudphxx
 
         typedef thrust::tuple<
           real_t,                       // rhod
-          real_t                        // eta
+          real_t,                       // eta
+          real_t                        // tke dissipation rate
         > tpl_ro_calc_t;
 
         tpl_rw_t tpl_rw;