Merge pull request #386 from pdziekan/sgs_length_profile

SGS_mixing_length as a profile

Author: pdziekan

bindings/python/lgrngn.hpp

@@ -345,6 +345,26 @@ namespace libcloudphxx
       {
         throw std::runtime_error("w_LS profile does not feature a getter yet - TODO");
       }
+
+      template <typename real_t>
+      void set_SGS_mix_len(
+        lgr::opts_init_t<real_t> *arg,
+        const bp_array &vec
+      )
+      {
+        sanity_checks(vec);
+        arg->SGS_mix_len.resize(0);
+        for (int i = 0; i < len(vec); ++i)
+          arg->SGS_mix_len.push_back(bp::extract<real_t>(vec[i]));
+      }
+
+      template <typename real_t>
+      bp_array get_SGS_mix_len(
+        lgr::opts_init_t<real_t> *arg
+      )
+      {
+        throw std::runtime_error("SGS_mix_len profile does not feature a getter yet - TODO");
+      }
     };
   };
 };

bindings/python/lib.cpp

@@ -214,10 +214,6 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .value("beard77fast", lgr::vt_t::beard77fast)
       .value("khvorostyanov_spherical", lgr::vt_t::khvorostyanov_spherical)
       .value("khvorostyanov_nonspherical", lgr::vt_t::khvorostyanov_nonspherical);
-    bp::enum_<lgr::SGS_length_scale_t::SGS_length_scale_t>("SGS_length_scale_t") 
-      .value("vertical", lgr::SGS_length_scale_t::vertical)
-      .value("arithmetic_mean", lgr::SGS_length_scale_t::arithmetic_mean)
-      .value("geometric_mean", lgr::SGS_length_scale_t::geometric_mean);
     bp::enum_<lgr::RH_formula_t::RH_formula_t>("RH_formula_t") 
       .value("pv_cc", lgr::RH_formula_t::pv_cc)
       .value("rv_cc", lgr::RH_formula_t::rv_cc)
@@ -296,12 +292,12 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("src_sd_conc", &lgr::opts_init_t<real_t>::src_sd_conc)
       .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("SGS_length_scale", &lgr::opts_init_t<real_t>::SGS_length_scale)
       .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)
       .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>)
     ;
     bp::class_<lgr::particles_proto_t<real_t>/*, boost::noncopyable*/>("particles_proto_t")

include/libcloudph++/lgrngn/SGS_length_scale.hpp

@@ -10,7 +10,6 @@
 #include "extincl.hpp"
 #include "kernel.hpp"
 #include "terminal_velocity.hpp"
-#include "SGS_length_scale.hpp"
 #include "advection_scheme.hpp"
 #include "RH_formula.hpp"
 #include "../common/chem.hpp"
@@ -88,9 +87,6 @@ namespace libcloudphxx
       // terminal velocity formula
       vt_t::vt_t terminal_velocity;
 
-      // SGS mixing length
-      SGS_length_scale_t::SGS_length_scale_t SGS_length_scale;
-
       // super-droplet advection scheme
       as_t::as_t adve_scheme;
 
@@ -130,6 +126,9 @@ namespace libcloudphxx
       // subsidence rate profile, positive downwards [m/s]
       std::vector<real_t> w_LS;
 
+      // SGS mixing length profile [m]
+      std::vector<real_t> SGS_mix_len;
+
       real_t rd_min; // minimal dry radius of droplets (works only for init from spectrum)
 
       // ctor with defaults (C++03 compliant) ...
@@ -158,7 +157,6 @@ namespace libcloudphxx
         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),
         terminal_velocity(vt_t::undefined),
-        SGS_length_scale(SGS_length_scale_t::geometric_mean),
         kernel(kernel_t::undefined),
         adve_scheme(as_t::implicit),
         RH_formula(RH_formula_t::pv_cc),

include/libcloudph++/lgrngn/opts_init.hpp

@@ -14,8 +14,6 @@
 #include <boost/numeric/odeint/external/thrust/thrust_operations.hpp>
 #include <boost/numeric/odeint/external/thrust/thrust_resize.hpp>
 
-#include <libcloudph++/common/SGS_length_scale.hpp>
-
 #include <map>
 
 namespace libcloudphxx
@@ -150,9 +148,8 @@ namespace libcloudphxx
         eta,// dynamic viscosity 
         diss_rate; // turbulent kinetic energy dissipation rate
 
-      real_t lambda;
-
       thrust_device::vector<real_t> w_LS; // large-scale subsidence velocity profile
+      thrust_device::vector<real_t> SGS_mix_len; // SGS mixing length profile
 
       // sorting needed only for diagnostics and coalescence
       bool sorted;
@@ -293,6 +290,7 @@ namespace libcloudphxx
                         halo_size * (opts_init.nz + 1) * opts_init.ny   // 3D
         ),
         w_LS(opts_init.w_LS),
+        SGS_mix_len(opts_init.SGS_mix_len),
         adve_scheme(opts_init.adve_scheme),
         pure_const_multi (((opts_init.sd_conc) == 0) && (opts_init.sd_const_multi > 0 || opts_init.dry_sizes.size() > 0)) // coal prob can be greater than one only in sd_conc simulations
       {
@@ -311,29 +309,6 @@ namespace libcloudphxx
 #endif
         *increase_sstp_coal = false;
 
-        switch (opts_init.SGS_length_scale)
-        {
-          case SGS_length_scale_t::vertical:
-            lambda =  
-              n_dims == 1 ? common::SGS_length_scale::vertical(opts_init.dx * si::metres)                                                      / si::metres: // 1D
-              n_dims == 2 ? common::SGS_length_scale::vertical(opts_init.dx * si::metres, opts_init.dz * si::metres)                           / si::metres: // 2D
-                            common::SGS_length_scale::vertical(opts_init.dx * si::metres, opts_init.dy * si::metres, opts_init.dz * si::metres)/ si::metres; // 3D
-            break;
-          case SGS_length_scale_t::arithmetic_mean:
-            lambda =  
-              n_dims == 1 ? common::SGS_length_scale::arithmetic_mean(opts_init.dx * si::metres)                                                      / si::metres: // 1D
-              n_dims == 2 ? common::SGS_length_scale::arithmetic_mean(opts_init.dx * si::metres, opts_init.dz * si::metres)                           / si::metres: // 2D
-                            common::SGS_length_scale::arithmetic_mean(opts_init.dx * si::metres, opts_init.dy * si::metres, opts_init.dz * si::metres)/ si::metres; // 3D
-            break;
-          case SGS_length_scale_t::geometric_mean:
-            lambda =  
-              n_dims == 1 ? common::SGS_length_scale::geometric_mean(opts_init.dx * si::metres)                                                      / si::metres: // 1D
-              n_dims == 2 ? common::SGS_length_scale::geometric_mean(opts_init.dx * si::metres, opts_init.dz * si::metres)                           / si::metres: // 2D
-                            common::SGS_length_scale::geometric_mean(opts_init.dx * si::metres, opts_init.dy * si::metres, opts_init.dz * si::metres)/ si::metres; // 3D
-            break;
-          default: assert(false && "unrecognized value of opts_init.SGS_length_scale"); 
-        }
-
         // initialising host temporary arrays
         {
           thrust_size_t n_grid;

src/impl/particles_impl.ipp

@@ -16,24 +16,32 @@ namespace libcloudphxx
       template<class real_t>
       struct common__turbulence__tke
       {
-        const quantity<si::length, real_t> lambda;
-        common__turbulence__tke(const real_t &lambda):
-          lambda(lambda * si::metres){}
-
         BOOST_GPU_ENABLED
-        real_t operator()(const real_t &diss_rate)
+        real_t operator()(const real_t &diss_rate, const real_t &lambda)
         {
+          assert(lambda > 0);
           return common::GA17_turbulence::tke(
             diss_rate * si::metres * si::metres / si::seconds / si::seconds / si::seconds,
-            lambda) / si::metres / si::metres * si::seconds * si::seconds;
+            lambda * si::metres) / si::metres / si::metres * si::seconds * si::seconds;
         }
       };
     };
     // calc the SGS TKE, in place of dissipation rate
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::hskpng_tke()
     {   
-      thrust::transform(diss_rate.begin(), diss_rate.end(), diss_rate.begin(), detail::common__turbulence__tke<real_t>(lambda));
+      namespace arg = thrust::placeholders;
+
+      thrust::transform(
+        diss_rate.begin(), diss_rate.end(),
+        thrust::make_permutation_iterator(SGS_mix_len.begin(),   // profile of the SGS mixing length
+          thrust::make_transform_iterator(                       // calculate vertical index from cell index
+            thrust::make_counting_iterator<thrust_size_t>(0),
+            arg::_1 % opts_init.nz
+          )
+        ),
+        diss_rate.begin(),
+        detail::common__turbulence__tke<real_t>());
     }
   };
 };

src/impl/particles_impl_hskpng_tke.ipp

@@ -16,16 +16,13 @@ namespace libcloudphxx
       template<class real_t>
       struct common__turbulence__tau
       {
-        const quantity<si::length, real_t> lambda;
-        common__turbulence__tau(const real_t &lambda):
-          lambda(lambda * si::metres){}
-
         BOOST_GPU_ENABLED
-        real_t operator()(const real_t &tke)
+        real_t operator()(const real_t &tke, const real_t &lambda)
         {
+          assert(lambda > 0);
           return common::GA17_turbulence::tau(
             tke * si::metres * si::metres / si::seconds / si::seconds,
-            lambda) / si::seconds;
+            lambda * si::metres) / si::seconds;
         }
       };
 
@@ -54,9 +51,21 @@ namespace libcloudphxx
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::hskpng_turb_vel(const bool only_vertical)
     {   
+      namespace arg = thrust::placeholders;
+
       thrust_device::vector<real_t> &tau(tmp_device_real_cell);
       thrust_device::vector<real_t> &tke(diss_rate); // should be called after hskpng_tke, which replaces diss_rate with tke
-      thrust::transform(tke.begin(), tke.end(), tau.begin(), detail::common__turbulence__tau<real_t>(lambda));
+
+      thrust::transform(tke.begin(), tke.end(),
+        thrust::make_permutation_iterator(SGS_mix_len.begin(),   // profile of the SGS mixing length
+          thrust::make_transform_iterator(                       // calculate vertical index from cell index
+            thrust::make_counting_iterator<thrust_size_t>(0),
+            arg::_1 % opts_init.nz
+          )
+        ),
+        tau.begin(),
+        detail::common__turbulence__tau<real_t>()
+      );
 
       thrust_device::vector<real_t> &r_normal(tmp_device_real_part);
       thrust_device::vector<real_t> * vel_turbs_vctrs_a[] = {&up, &wp, &vp};

src/impl/particles_impl_hskpng_turb_vel.ipp

@@ -74,39 +74,47 @@ namespace libcloudphxx
       if(opts_init.sd_const_multi > 0 && opts_init.src_switch)
         throw std::runtime_error("aerosol source and constant multiplicity option are not compatible");
 
-        if (n_dims > 0)
-        {
-          if (!(opts_init.x0 >= 0 && opts_init.x0 < m1(opts_init.nx) * opts_init.dx))
-            throw std::runtime_error("!(x0 >= 0 & x0 < min(1,nx)*dz)");
-          if (!(opts_init.y0 >= 0 && opts_init.y0 < m1(opts_init.ny) * opts_init.dy))
-            throw std::runtime_error("!(y0 >= 0 & y0 < min(1,ny)*dy)");
-          if (!(opts_init.z0 >= 0 && opts_init.z0 < m1(opts_init.nz) * opts_init.dz))
-            throw std::runtime_error("!(z0 >= 0 & z0 < min(1,nz)*dz)");
-          // check temporarily disabled since dewv_id is not passed anymore, TODO: fix it
-//        if (!(opts_init.x1 > opts_init.x0 && opts_init.x1 <= m1(opts_init.nx) * opts_init.dx) && dev_id == -1) // only for single device runs, since on multi_CUDA x1 is not yet adjusted to local domain
-//            throw std::runtime_error("!(x1 > x0 & x1 <= min(1,nx)*dx)");
-          if (!(opts_init.y1 > opts_init.y0 && opts_init.y1 <= m1(opts_init.ny) * opts_init.dy))
-            throw std::runtime_error("!(y1 > y0 & y1 <= min(1,ny)*dy)");
-          if (!(opts_init.z1 > opts_init.z0 && opts_init.z1 <= m1(opts_init.nz) * opts_init.dz))
-            throw std::runtime_error("!(z1 > z0 & z1 <= min(1,nz)*dz)");
-        }
-
-        if (opts_init.dt == 0) throw std::runtime_error("please specify opts_init.dt");
-        if (opts_init.sd_conc * opts_init.sd_const_multi != 0) throw std::runtime_error("specify either opts_init.sd_conc or opts_init.sd_const_multi, not both");
-        if (opts_init.sd_conc == 0 && opts_init.sd_const_multi == 0 && opts_init.dry_sizes.size() == 0) throw std::runtime_error("please specify opts_init.sd_conc, opts_init.sd_const_multi or opts_init.dry_sizes");
-        if (opts_init.coal_switch)
-        {
-          if(opts_init.terminal_velocity == vt_t::undefined) throw std::runtime_error("please specify opts_init.terminal_velocity or turn off opts_init.coal_switch");
-          if(opts_init.kernel == kernel_t::undefined) throw std::runtime_error("please specify opts_init.kernel");
-        }
-        if (opts_init.sedi_switch)
-          if(opts_init.terminal_velocity == vt_t::undefined) throw std::runtime_error("please specify opts_init.terminal_velocity or turn off opts_init.sedi_switch");
-        if (opts_init.sedi_switch && opts_init.nz == 0)
-          throw std::runtime_error("opts_init.sedi_switch can be True only if n_dims > 1");
-        if (opts_init.subs_switch && opts_init.nz == 0)
-          throw std::runtime_error("opts_init.subs_switch can be True only if n_dims > 1");
-        if (opts_init.subs_switch && opts_init.nz != w_LS.size())
-          throw std::runtime_error("opts_init.subs_switch == True, but subsidence velocity profile size != nz");
+      if (n_dims > 0)
+      {
+        if (!(opts_init.x0 >= 0 && opts_init.x0 < m1(opts_init.nx) * opts_init.dx))
+          throw std::runtime_error("!(x0 >= 0 & x0 < min(1,nx)*dz)");
+        if (!(opts_init.y0 >= 0 && opts_init.y0 < m1(opts_init.ny) * opts_init.dy))
+          throw std::runtime_error("!(y0 >= 0 & y0 < min(1,ny)*dy)");
+        if (!(opts_init.z0 >= 0 && opts_init.z0 < m1(opts_init.nz) * opts_init.dz))
+          throw std::runtime_error("!(z0 >= 0 & z0 < min(1,nz)*dz)");
+        // check temporarily disabled since dewv_id is not passed anymore, TODO: fix it
+//      if (!(opts_init.x1 > opts_init.x0 && opts_init.x1 <= m1(opts_init.nx) * opts_init.dx) && dev_id == -1) // only for single device runs, since on multi_CUDA x1 is not yet adjusted to local domain
+//          throw std::runtime_error("!(x1 > x0 & x1 <= min(1,nx)*dx)");
+        if (!(opts_init.y1 > opts_init.y0 && opts_init.y1 <= m1(opts_init.ny) * opts_init.dy))
+          throw std::runtime_error("!(y1 > y0 & y1 <= min(1,ny)*dy)");
+        if (!(opts_init.z1 > opts_init.z0 && opts_init.z1 <= m1(opts_init.nz) * opts_init.dz))
+          throw std::runtime_error("!(z1 > z0 & z1 <= min(1,nz)*dz)");
+      }
+
+      if (opts_init.dt == 0) throw std::runtime_error("please specify opts_init.dt");
+      if (opts_init.sd_conc * opts_init.sd_const_multi != 0) throw std::runtime_error("specify either opts_init.sd_conc or opts_init.sd_const_multi, not both");
+      if (opts_init.sd_conc == 0 && opts_init.sd_const_multi == 0 && opts_init.dry_sizes.size() == 0) throw std::runtime_error("please specify opts_init.sd_conc, opts_init.sd_const_multi or opts_init.dry_sizes");
+      if (opts_init.coal_switch)
+      {
+        if(opts_init.terminal_velocity == vt_t::undefined) throw std::runtime_error("please specify opts_init.terminal_velocity or turn off opts_init.coal_switch");
+        if(opts_init.kernel == kernel_t::undefined) throw std::runtime_error("please specify opts_init.kernel");
+      }
+      if (opts_init.sedi_switch)
+        if(opts_init.terminal_velocity == vt_t::undefined) throw std::runtime_error("please specify opts_init.terminal_velocity or turn off opts_init.sedi_switch");
+      if (opts_init.sedi_switch && opts_init.nz == 0)
+        throw std::runtime_error("opts_init.sedi_switch can be True only if n_dims > 1");
+      if (opts_init.subs_switch && opts_init.nz == 0)
+        throw std::runtime_error("opts_init.subs_switch can be True only if n_dims > 1");
+      if (opts_init.turb_adve_switch && opts_init.nz == 0)
+        throw std::runtime_error("opts_init.turb_adve_switch can be True only if n_dims > 1");
+      if (opts_init.turb_cond_switch && opts_init.nz == 0)
+        throw std::runtime_error("opts_init.turb_cond_switch can be True only if n_dims > 1");
+      if (opts_init.subs_switch && opts_init.nz != w_LS.size())
+        throw std::runtime_error("opts_init.subs_switch == True, but subsidence velocity profile size != nz");
+      if ((opts_init.turb_adve_switch || opts_init.turb_cond_switch) && opts_init.nz != SGS_mix_len.size())
+        throw std::runtime_error("at least one of opts_init.turb_adve_switch, opts_init.turb_cond_switch is true, but SGS mixing length profile size != nz");
+      if(opts_init.SGS_mix_len.size() > 0 && *std::min(opts_init.SGS_mix_len.begin(), opts_init.SGS_mix_len.end()) <= 0)
+        throw std::runtime_error("SGS_mix_len <= 0");
     }
   };
 };

src/impl/particles_impl_init_sanity_check.ipp

@@ -1,5 +1,5 @@
 # non-pytest tests
-foreach(test api_blk_1m api_blk_2m api_lgrngn api_common segfault_20150216 col_kernels terminal_velocities SGS_length_scales SD_removal uniform_init source chem_coal sstp_cond multiple_kappas adve_scheme lgrngn_subsidence sat_adj_blk_1m )
+foreach(test api_blk_1m api_blk_2m api_lgrngn api_common segfault_20150216 col_kernels terminal_velocities SD_removal uniform_init source chem_coal sstp_cond multiple_kappas adve_scheme lgrngn_subsidence sat_adj_blk_1m )
   #TODO: indicate that tests depend on the lib
   add_test(
     NAME ${test}