Merge pull request #416 from pdziekan/variable_dt

Variable dt

Author: pdziekan

bindings/python/lib.cpp

@@ -250,6 +250,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .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)
+      .def_readwrite("dt", &lgr::opts_t<real_t>::dt)
     ;
     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>)
@@ -310,6 +311,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .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)
     ;
     bp::class_<lgr::particles_proto_t<real_t>/*, boost::noncopyable*/>("particles_proto_t")
       .add_property("opts_init", &lgrngn::get_oi<real_t>)

include/libcloudph++/lgrngn/opts.hpp

@@ -28,12 +28,16 @@ namespace libcloudphxx
       // process toggling for chemistry
       bool chem_dsl, chem_dsc, chem_rct;
 
+      // overriding dt from opts_init
+      real_t dt;
+
       // ctor with defaults (C++03 compliant) ...
       opts_t() : 
         adve(true), sedi(true), subs(false), 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
+        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

@@ -64,6 +64,9 @@ namespace libcloudphxx
       // should aerosol concentration init be independent of rhod (assumed to be in cm^{-3} and not at STP)
       bool aerosol_independent_of_rhod;
 
+      // is it allowed to change dt during simulation through opts.dt
+      bool variable_dt_switch;
+
       // or, alternatively to sd_conc_mean, multiplicity of all SDs = const
       unsigned long long  sd_const_multi;
 
@@ -187,7 +190,8 @@ namespace libcloudphxx
         diag_incloud_time(false),
         no_ccn_at_init(false),
         open_side_walls(false),
-        periodic_topbot_walls(false)
+        periodic_topbot_walls(false),
+        variable_dt_switch(false)
       {}
 
       // dtor (just to silence -Winline warnings)

src/impl/particles_impl.ipp

@@ -152,6 +152,10 @@ namespace libcloudphxx
       thrust_device::vector<real_t> w_LS; // large-scale subsidence velocity profile
       thrust_device::vector<real_t> SGS_mix_len; // SGS mixing length profile
 
+      // time steps to be used, considering that opts.dt can override opts_init.dt
+      real_t dt;
+      int sstp_cond, sstp_coal, sstp_chem;
+
       // sorting needed only for diagnostics and coalescence
       bool sorted;
 
@@ -163,6 +167,10 @@ namespace libcloudphxx
       // is a constant, external pressure profile used? (e.g. anelastic model)
       bool const_p;
 
+      // is it allowed to do substepping, if not, some memory can be saved
+      bool allow_sstp_cond,
+           allow_sstp_chem;
+
       // timestep counter
       n_t stp_ctr;
 
@@ -327,6 +335,8 @@ namespace libcloudphxx
         w_LS(_opts_init.w_LS),
         SGS_mix_len(_opts_init.SGS_mix_len),
         adve_scheme(_opts_init.adve_scheme),
+        allow_sstp_cond(_opts_init.sstp_cond > 1 || _opts_init.variable_dt_switch),
+        allow_sstp_chem(_opts_init.sstp_chem > 1 || _opts_init.variable_dt_switch),
         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
       {
 
@@ -383,7 +393,7 @@ namespace libcloudphxx
         if (opts_init.ny != 0)  distmem_real_vctrs.insert(&y);
         if (opts_init.nz != 0)  distmem_real_vctrs.insert(&z);
 
-        if(opts_init.sstp_cond > 1 && opts_init.exact_sstp_cond)
+        if(allow_sstp_cond && opts_init.exact_sstp_cond)
         {
            distmem_real_vctrs.insert(&sstp_tmp_rv);
            distmem_real_vctrs.insert(&sstp_tmp_th);
@@ -483,7 +493,7 @@ namespace libcloudphxx
       void hskpng_vterm_all();
       void hskpng_vterm_invalid();
       void hskpng_tke();
-      void hskpng_turb_vel(const bool only_vertical = false);
+      void hskpng_turb_vel(const real_t &dt, const bool only_vertical = false);
       void hskpng_turb_dot_ss();
       void hskpng_remove_n0();
       void hskpng_resize_npart();
@@ -522,12 +532,13 @@ namespace libcloudphxx
         arrinfo_t<real_t> &
       );
 
+      void adjust_timesteps(const real_t &dt);
       void adve();
-      void turb_adve();
+      void turb_adve(const real_t &dt);
       template<class adve_t>
       void adve_calc(bool, thrust_size_t = 0);
-      void sedi();
-      void subs();
+      void sedi(const real_t &dt);
+      void subs(const real_t &dt);
 
       void cond_dm3_helper();
       void cond(const real_t &dt, const real_t &RH_max, const bool turb_cond);
@@ -537,7 +548,7 @@ namespace libcloudphxx
       void update_th_rv(thrust_device::vector<real_t> &);
       void update_state(thrust_device::vector<real_t> &, thrust_device::vector<real_t> &);
       void update_pstate(thrust_device::vector<real_t> &, thrust_device::vector<real_t> &);
-      void update_incloud_time();
+      void update_incloud_time(const real_t &dt);
 
       void coal(const real_t &dt, const bool &turb_coal);
 

src/impl/particles_impl_adjust_timesteps.ipp

@@ -0,0 +1,24 @@
+// vim:filetype=cpp
+/** @file
+  * @copyright University of Warsaw
+  * @section LICENSE
+  * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/)
+  */
+
+namespace libcloudphxx
+{
+  namespace lgrngn
+  {
+    template <typename real_t, backend_t device>
+    void particles_t<real_t, device>::impl::adjust_timesteps(const real_t &_dt)
+    {
+      if(_dt > 0 && !opts_init.variable_dt_switch) throw std::runtime_error("opts.dt specified, but opts_init.variable_dt_switch is false.");
+      // dt defined in opts_init can be overriden by dt passed to this function
+      dt = _dt > 0 ? _dt : opts_init.dt;
+      // then, number of substeps is adjusted to get desired dt for specific processes
+      sstp_cond = _dt > 0 ? std::ceil(opts_init.sstp_cond * _dt / opts_init.dt) : opts_init.sstp_cond;
+      sstp_coal = _dt > 0 ? std::ceil(opts_init.sstp_coal * _dt / opts_init.dt) : opts_init.sstp_coal;
+      sstp_chem = _dt > 0 ? std::ceil(opts_init.sstp_chem * _dt / opts_init.dt) : opts_init.sstp_chem;
+    }
+  };
+};

src/impl/particles_impl_hskpng_resize.ipp

@@ -49,16 +49,16 @@ namespace libcloudphxx
         dot_ssp.resize(n_part, 0);
       }
 
-      if(opts_init.chem_switch || opts_init.sstp_cond > 1 || n_dims >= 2)
+      if(opts_init.chem_switch || allow_sstp_cond || n_dims >= 2)
       {
         tmp_device_real_part1.resize(n_part);
       }
-      if((opts_init.sstp_cond>1 && opts_init.exact_sstp_cond) || n_dims==3 || opts_init.turb_cond_switch)
+      if((allow_sstp_cond && opts_init.exact_sstp_cond) || n_dims==3 || opts_init.turb_cond_switch)
       {
         tmp_device_real_part2.resize(n_part);
       }
 
-      if(opts_init.sstp_cond>1 && opts_init.exact_sstp_cond)
+      if(allow_sstp_cond && opts_init.exact_sstp_cond)
       {
         tmp_device_real_part3.resize(n_part);
         tmp_device_real_part4.resize(n_part);  

src/impl/particles_impl_hskpng_turb_vel.ipp

@@ -49,7 +49,7 @@ namespace libcloudphxx
     };
     // calc the SGS turbulent velocity component
     template <typename real_t, backend_t device>
-    void particles_t<real_t, device>::impl::hskpng_turb_vel(const bool only_vertical)
+    void particles_t<real_t, device>::impl::hskpng_turb_vel(const real_t &dt, const bool only_vertical)
     {   
       namespace arg = thrust::placeholders;
 
@@ -86,7 +86,7 @@ namespace libcloudphxx
             r_normal.begin()
           )) + n_part,
           vel_turbs_vctrs_a[i]->begin(),
-          detail::common__turbulence__update_turb_vel<real_t>(opts_init.dt)
+          detail::common__turbulence__update_turb_vel<real_t>(dt)
         );
       }
     }

src/impl/particles_impl_init_hskpng_ncell.ipp

@@ -29,7 +29,7 @@ namespace libcloudphxx
       tmp_device_size_cell.resize(n_cell);
       tmp_host_size_cell.resize(n_cell);
       tmp_host_real_cell.resize(n_cell);
-      if(opts_init.sstp_cond > 1 && !opts_init.exact_sstp_cond)
+      if(allow_sstp_cond && !opts_init.exact_sstp_cond)
       {
         sstp_tmp_rv.resize(n_cell);
         sstp_tmp_th.resize(n_cell);

src/impl/particles_impl_reserve_hskpng_npart.ipp

@@ -56,11 +56,11 @@ namespace libcloudphxx
       {
         tmp_device_real_part1.reserve(opts_init.n_sd_max); 
       }
-      if((opts_init.sstp_cond>1 && opts_init.exact_sstp_cond) || n_dims==3 || opts_init.turb_cond_switch)
+      if((allow_sstp_cond && opts_init.exact_sstp_cond) || n_dims==3 || opts_init.turb_cond_switch)
       {
         tmp_device_real_part2.reserve(opts_init.n_sd_max); 
       }
-      if(opts_init.sstp_cond>1 && opts_init.exact_sstp_cond)
+      if(allow_sstp_cond && opts_init.exact_sstp_cond)
       {
         tmp_device_real_part3.reserve(opts_init.n_sd_max); 
         tmp_device_real_part4.reserve(opts_init.n_sd_max);  
@@ -107,7 +107,7 @@ namespace libcloudphxx
 //      }
 //      vt.resize(opts_init.n_sd_max, 0.); // so that it may be safely used in condensation before first update
 //
-//      if(opts_init.sstp_cond>1 && opts_init.exact_sstp_cond)
+//      if(allow_sstp_cond && opts_init.exact_sstp_cond)
 //      {
 //        sstp_tmp_rv.resize(opts_init.n_sd_max);
 //        sstp_tmp_th.resize(opts_init.n_sd_max);

src/impl/particles_impl_sedi.ipp

@@ -10,7 +10,7 @@ namespace libcloudphxx
   namespace lgrngn
   {
     template <typename real_t, backend_t device>
-    void particles_t<real_t, device>::impl::sedi()
+    void particles_t<real_t, device>::impl::sedi(const real_t &dt)
     {   
       namespace arg = thrust::placeholders;
 
@@ -19,7 +19,7 @@ namespace libcloudphxx
         z.begin(), z.end(),                    // position
         vt.begin(),                                                    // terminal velocity 
         z.begin(),                         // output
-        arg::_1 - opts_init.dt * arg::_2   // Euler scheme (assuming vt positive!)
+        arg::_1 - dt * arg::_2   // Euler scheme (assuming vt positive!)
       );
     }
   };