Adaptive condensation substepping

CMakeLists.txt

@@ -287,6 +287,8 @@ if (CMAKE_CUDA_COMPILER)
     target_compile_options(cloudphxx_lgrngn PRIVATE  $<$<COMPILE_LANGUAGE:CUDA>: -Xcompiler -DUSE_MPI>)
   endif()
 
+#  target_compile_options(cloudphxx_lgrngn PRIVATE  $<$<COMPILE_LANGUAGE:CUDA>: --extended-lambda>)
+
   target_compile_definitions(cloudphxx_lgrngn PRIVATE CUDA_FOUND)
 
   # why was this added in the MPI PR? commented-out for now

bindings/python/lib.cpp

@@ -321,7 +321,10 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .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_mix", &lgr::opts_init_t<real_t>::sstp_cond_mix)
       .def_readwrite("sstp_cond", &lgr::opts_init_t<real_t>::sstp_cond)
+      .def_readwrite("sstp_cond_act", &lgr::opts_init_t<real_t>::sstp_cond_act)
+      .def_readwrite("rc2_T", &lgr::opts_init_t<real_t>::rc2_T)
       .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_rlx", &lgr::opts_init_t<real_t>::supstp_rlx)
@@ -352,6 +355,11 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("time_dep_ice_nucl", &lgr::opts_init_t<real_t>::time_dep_ice_nucl)
       .def_readwrite("const_p", &lgr::opts_init_t<real_t>::const_p)
       .def_readwrite("th_dry", &lgr::opts_init_t<real_t>::th_dry)
+      .def_readwrite("sstp_cond_mix", &lgr::opts_init_t<real_t>::sstp_cond_mix)
+      .def_readwrite("adaptive_sstp_cond", &lgr::opts_init_t<real_t>::adaptive_sstp_cond)
+      .def_readwrite("sstp_cond_adapt_drw2_eps", &lgr::opts_init_t<real_t>::sstp_cond_adapt_drw2_eps)
+      .def_readwrite("sstp_cond_adapt_drw2_max", &lgr::opts_init_t<real_t>::sstp_cond_adapt_drw2_max)
+
     ;
     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_init.hpp

@@ -38,8 +38,16 @@ namespace libcloudphxx
       int nx, ny, nz;
       real_t dx, dy, dz, dt;
 
-      // no. of substeps 
+      // no. of substeps for condensation/coalescence. If adaptive_sstp_cond=true, sstp_cond is the maximum no. of substeps.
       int sstp_cond, sstp_coal; 
+      // no. of condensation substeps for SDs that activate/deactivate in this timestep
+      // only work in adaptive_sstp_cond mode.
+      int sstp_cond_act; 
+      // no. of substeps for chemistry
+      int sstp_chem;
+      // Note: If dt changes during simulation (by supplying opts.dt), 
+      //       no. of substeps is adjusted accordingly to keep process timestep close to dt / sstp_cond, 
+      //       but only if initially sstp_cond/coal/chem/cond_act > 1
 
       // Lagrangian domain extents
       real_t x0, y0, z0, x1, y1, z1;
@@ -90,12 +98,16 @@ namespace libcloudphxx
            turb_coal_switch,   // if true, turbulent coalescence kernels can be used
            ice_switch,         // if true, ice is allowed
            exact_sstp_cond,    // if true, use per-particle sstp_cond logic, if false, use per-cell
+           sstp_cond_mix,      // if true, th and rv of all SDs in a cell after each timestep (instant mixing at substep timescale), else update it only after all substeps
+           adaptive_sstp_cond, // if true, use adaptive number of substeps for condensation
            time_dep_ice_nucl;  // it true, time-dependent freezing, if false, singular freezing
 
+      real_t sstp_cond_adapt_drw2_eps = 1e-4; // tolerance for adaptive substepping in condensation (drw2_err <= sstp_cond_adapt_eps * rw2)
+      real_t sstp_cond_adapt_drw2_max = 4; // tolerance for adaptive substepping in condensation (drw2 < sstp_cond_adapt_drw2_max * rw2)
+
       // ice nucleating particles type
       INP_t inp_type;
 
-      int sstp_chem;
       real_t chem_rho;
 
       // do we want to track the time SDs spend inside clouds
@@ -134,6 +146,7 @@ namespace libcloudphxx
       bool open_side_walls,       // if true, side walls are "open", i.e. SD are removed at contact. Periodic otherwise.
            periodic_topbot_walls; // if true, top and bot walls are periodic. Open otherwise
 
+      real_t rc2_T = 10; // temperature [C] at which rc2 is calculated
 
       // --- aerosol source stuff ---
 
@@ -188,7 +201,7 @@ namespace libcloudphxx
         aerosol_independent_of_rhod(false), 
         sd_const_multi(0),
         dt(0),   
-        sstp_cond(1), sstp_coal(1), sstp_chem(1),         
+        sstp_cond(1), sstp_coal(1), sstp_chem(1), sstp_cond_act(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
@@ -199,6 +212,8 @@ namespace libcloudphxx
         time_dep_ice_nucl(false),
         inp_type(INP_t::mineral),
         exact_sstp_cond(false),
+        sstp_cond_mix(true),
+        adaptive_sstp_cond(false),
         turb_cond_switch(false),
         turb_adve_switch(false),
         turb_coal_switch(false),

src/detail/tmp_vector_pool.hpp

@@ -16,15 +16,15 @@ namespace libcloudphxx
             entry(size_t n) : vec(n) {}
         };
         std::vector<entry> pool;
+        const std::string name;
+
     public:
-        tmp_vector_pool(size_t pool_size = 1): pool(pool_size, 0) {}
+        tmp_vector_pool(std::string name, size_t pool_size = 1): pool(pool_size, 0), name(name) {}
 
-        // Add a new vector to the pool
-        // void add_vector(size_t vec_size) {
-        //     pool.emplace_back(vec_size);
-        // }
-        void add_vector() {
-            pool.emplace_back(0);
+        void add_vectors(size_t no_vectors = 1) {
+            for (size_t i = 0; i < no_vectors; ++i) {
+                pool.emplace_back(0);
+            }
         }
 
         void resize(size_t n) {
@@ -41,14 +41,15 @@ namespace libcloudphxx
 
         // Acquire an available vector, returns its index
         size_t acquire() {
-//            std::cerr << "tmp_vector_pool: acquiring vector from pool of size " << pool.size() << "\n";
             for (size_t i = 0; i < pool.size(); ++i) {
                 if (!pool[i].in_use) {
                     pool[i].in_use = true;
                     return i;
                 }
             }
-            assert(false && "No available temporary vectors in pool!");
+            std::cerr << "LIBCLOUDPH++: tmp_vector_pool: No available temporary vectors in pool! Pool size: " << pool.size() << ". Pool name: " << name << " \n";
+            // assert(false && "No available temporary vectors in pool!");
+            assert(false);
             return size_t(-1);
         }
 

src/impl/particles_impl_sstp_chem.ipp → src/impl/chemistry/particles_impl_sstp_chem.ipp

@@ -33,7 +33,7 @@ namespace libcloudphxx
     }
 
     template <typename real_t, backend_t device>
-    void particles_t<real_t, device>::impl::init_sstp_chem()
+    void particles_t<real_t, device>::impl::init_percell_sstp_chem()
     {   
       if (!allow_sstp_chem) return;
 
@@ -50,7 +50,7 @@ namespace libcloudphxx
     }
 
     template <typename real_t, backend_t device>
-    void particles_t<real_t, device>::impl::sstp_step_chem(
+    void particles_t<real_t, device>::impl::sstp_percell_step_chem(
       const int &step
     )
     {   

src/impl/particles_impl_coal.ipp → src/impl/coalescence/particles_impl_coal.ipp

@@ -30,6 +30,19 @@ namespace libcloudphxx
         }
       };
 
+      struct invalidator // set value of some parameter to invalid in the SD that had its parameters changed (other than multiplicity)
+      {
+        template<class tpl_t>
+        BOOST_GPU_ENABLED
+        void operator()(tpl_t tpl)
+        {
+          if(thrust::get<2>(tpl) <= 0) return; // do nothing if no collisions or first one passed was a SD with an uneven number in the cell
+          thrust::get<3>(tpl) == na_ge_nb ?    // does the first SD of the pair have greater multiplicity?
+            thrust::get<1>(tpl) = detail::invalid: 
+            thrust::get<0>(tpl) = detail::invalid;
+        }
+      };
+
       struct summator
       {
         template<class tpl_t>
@@ -510,6 +523,26 @@ namespace libcloudphxx
         );
         nancheck(incloud_time, "incloud_time - post coalescence");
       }
+
+      // invalidate rc2 due to changed rd3 and kappa
+      if(opts_init.sstp_cond_act > 1 && allow_sstp_cond) // rc2 only used in adaptive activation substepping
+      {
+        thrust::for_each(
+          thrust::make_zip_iterator(thrust::make_tuple(
+            thrust::make_permutation_iterator(rc2.begin(), sorted_id.begin()),   
+            thrust::make_permutation_iterator(rc2.begin(), sorted_id.begin())+1,
+            col.begin(),                                                        
+            col.begin()+1
+          )),
+          thrust::make_zip_iterator(thrust::make_tuple(
+            thrust::make_permutation_iterator(rc2.begin(), sorted_id.begin()),   
+            thrust::make_permutation_iterator(rc2.begin(), sorted_id.begin())+1,
+            col.begin(),                                                        
+            col.begin()+1
+          )) + n_part -1,
+          detail::invalidator()
+        );
+      }
     }
   };  
 };

src/impl/common/calc_liq_ice_content_change.ipp

@@ -0,0 +1,50 @@
+// 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::calc_liq_ice_content_change()
+    {   
+      // TODO: assert that save_liq_ice_content_before_change was called beforehand?
+      thrust_device::vector<real_t> &drw_mom3 = drw_mom3_gp->get();
+
+      moms_all();
+      moms_calc(rw2.begin(), real_t(3./2.));
+
+      // drw_mom3 = -rw_mom3_ante + rw_mom3_post
+      thrust::transform(
+        count_mom.begin(), count_mom.begin() + count_n,                    // input - 1st arg
+        thrust::make_permutation_iterator(drw_mom3.begin(), count_ijk.begin()), // 2nd arg
+        thrust::make_permutation_iterator(drw_mom3.begin(), count_ijk.begin()), // output
+        thrust::plus<real_t>()
+      );
+
+      if(opts_init.ice_switch)
+      {
+        thrust_device::vector<real_t> &d_ice_mass = d_ice_mass_gp->get();
+
+        moms_gt0(ice_a.begin()); // choose ice particles (ice_a>0)
+        moms_calc(thrust::make_transform_iterator(
+          thrust::make_zip_iterator(thrust::make_tuple(ice_a.begin(), ice_c.begin(), ice_rho.begin())),
+          detail::ice_mass<real_t>()
+          ),
+          real_t(1)
+        );
+
+        thrust::transform(
+          count_mom.begin(), count_mom.begin() + count_n,                    // input - 1st arg
+          thrust::make_permutation_iterator(d_ice_mass.begin(), count_ijk.begin()), // 2nd arg
+          thrust::make_permutation_iterator(d_ice_mass.begin(), count_ijk.begin()), // output
+          thrust::plus<real_t>()
+        );
+      }
+    }
+  };  
+};