Diagnose mean distance before coalescence

include/libcloudph++/lgrngn/opts_init.hpp

@@ -183,6 +183,9 @@ namespace libcloudphxx
       // relaxation time scale [s]
       real_t rlx_timescale;
 
+      // pre-collision statistics are diagnosed up to precoal_stats_tmax [s] before collision
+      real_t precoal_stats_tmax;
+
       // -- ctors ---
 
       // ctor with defaults (C++03 compliant) ...
@@ -237,7 +240,8 @@ namespace libcloudphxx
         open_side_walls(false),
         periodic_topbot_walls(false),
         variable_dt_switch(false),
-        rng_seed_init_switch(false)
+        rng_seed_init_switch(false),
+        precoal_stats_tmax(3600)
       {}
 
       // dtor (just to silence -Winline warnings)

include/libcloudph++/lgrngn/particles.hpp

@@ -69,7 +69,8 @@ namespace libcloudphxx
 
       // returns accumulated rainfall
       virtual void step_async(
-        const opts_t<real_t> &
+        const opts_t<real_t> &,
+        const real_t &
       ) { 
         assert(false); 
       }  
@@ -108,6 +109,8 @@ namespace libcloudphxx
       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 void store_ijk(const real_t&)                                     { assert(false); }
+      virtual std::vector<real_t> diag_precoal_distance()                       { 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; }
 
@@ -167,7 +170,8 @@ namespace libcloudphxx
       );
 
       void step_async(
-        const opts_t<real_t> &
+        const opts_t<real_t> &,
+        const real_t &
       );
 
       // diagnostic methods
@@ -197,6 +201,8 @@ namespace libcloudphxx
       void diag_precip_rate();
       void diag_max_rw();
       void diag_vel_div();
+      void store_ijk(const real_t&);
+      std::vector<real_t> diag_precoal_distance();
       std::map<libcloudphxx::common::output_t, real_t> diag_puddle();
       real_t *outbuf();
 
@@ -263,7 +269,8 @@ namespace libcloudphxx
       );
 
       void step_async(
-        const opts_t<real_t> &
+        const opts_t<real_t> &,
+        const real_t &
       );
 
       // diagnostic methods
@@ -294,6 +301,8 @@ namespace libcloudphxx
       void diag_precip_rate();
       void diag_max_rw();
       void diag_vel_div();
+      void store_ijk(const real_t&);
+      std::vector<real_t> diag_precoal_distance();
       std::map<libcloudphxx::common::output_t, real_t> diag_puddle();
 
       struct impl;

models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp

@@ -252,6 +252,7 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common<ct_params_t>
         using libcloudphxx::lgrngn::multi_CUDA;
 
 #if defined(STD_FUTURE_WORKS)
+/*
         if (params.async)
         {
           assert(!ftr.valid());
@@ -261,19 +262,22 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common<ct_params_t>
               std::launch::async, 
               &particles_t<real_t, multi_CUDA>::step_async, 
               dynamic_cast<particles_t<real_t, multi_CUDA>*>(prtcls.get()),
-              params.cloudph_opts
+              params.cloudph_opts,
+              this->timestep * params.dt
             );
           else if(params.backend == CUDA)
             ftr = std::async(
               std::launch::async, 
               &particles_t<real_t, CUDA>::step_async, 
               dynamic_cast<particles_t<real_t, CUDA>*>(prtcls.get()),
-              params.cloudph_opts
+              params.cloudph_opts,
+              this->timestep * params.dt
             );
           assert(ftr.valid());
         } else 
+        */
 #endif
-          prtcls->step_async(params.cloudph_opts);
+          prtcls->step_async(params.cloudph_opts, this->timestep * params.dt);
       }
 
       // performing diagnostics

models/kinematic_2D/src/kin_cloud_2d_lgrngn_chem.hpp

@@ -279,6 +279,7 @@ class kin_cloud_2d_lgrngn_chem : public kin_cloud_2d_lgrngn<ct_params_t>
         using libcloudphxx::lgrngn::CUDA;
         using libcloudphxx::lgrngn::multi_CUDA;
 
+/*
 #if defined(STD_FUTURE_WORKS)
         if (parent_t::params.async)
         {
@@ -302,6 +303,7 @@ class kin_cloud_2d_lgrngn_chem : public kin_cloud_2d_lgrngn<ct_params_t>
         } else 
 #endif
           parent_t::prtcls->step_async(parent_t::params.cloudph_opts);
+          */
       }
 
       // performing diagnostics

src/detail/config.hpp

@@ -32,13 +32,19 @@ namespace libcloudphxx
 
         const real_t rlx_conc_tolerance = 0.1; // tolerance of the relaxation scheme; new SD will be created if missing_conc/expected_conc > tolerance
 
+        // number of bins to store data about distance (and other attributes?) before collision
+        const int precoal_stats_bins = 100;
+
         // ctor
         config():
           vt0_ln_r_min(log(5e-7)),
           vt0_ln_r_max(log(3e-3)),  // Beard 1977 is defined on 1um - 6mm diameter range
           eps_tolerance(sizeof(real_t) * 8 / 4)
           {}
       };
+
+      // just some constant, not related to config but had to put them somewhere
+      enum { invalid = -1, no_initial_value = -44 };
     };
   };
 };

src/impl/particles_impl.ipp

@@ -23,12 +23,6 @@ namespace libcloudphxx
 {
   namespace lgrngn
   {
-    namespace detail
-    {   
-      enum { invalid = -1 };
-
-    };  
-
     // pimpl stuff 
     template <typename real_t, backend_t device>
     struct particles_t<real_t, device>::impl
@@ -91,6 +85,7 @@ namespace libcloudphxx
         sstp_tmp_p, // ditto for pressure
         incloud_time; // time this SD has been within a cloud
 
+
       // dry radii distribution characteristics
       real_t log_rd_min, // logarithm of the lower bound of the distr
              log_rd_max, // logarithm of the upper bound of the distr
@@ -230,6 +225,18 @@ namespace libcloudphxx
       // to simplify foreach calls
       const thrust::counting_iterator<thrust_size_t> zero;
 
+      // --- stuff for diagnosing distance (and other attributes?) before collision ---
+      // cell number of each SD at moments in time when store_ijk was called,
+      // TODO: with MPI, ijk_history needs to account for domains 
+      // NOTE: ijk_history would better be thrust_size_t type, not n_t, but then we couldnt use distmem_n_vctrs 
+      //       to manage it and to add distmem_size_vctrs would require adding more GPU/MPI communications as there currently are none for thrust_size_t
+      //std::vector<thrust_device::vector<thrust_size_t>> ijk_history;
+//      std::vector<thrust_device::vector<n_t>> ijk_history;
+      std::vector<std::unique_ptr<thrust_device::vector<n_t>>> ijk_history;
+      std::vector<real_t> ijk_history_time; // time at which history was saved
+      std::vector<real_t> precoal_distance;
+      std::vector<n_t> precoal_distance_count;
+
       // -- distributed memory stuff --
       // TODO: move to a separate struct?
 
@@ -270,11 +277,19 @@ namespace libcloudphxx
       // ids of sds to be copied with distmem
       thrust_device::vector<thrust_size_t> &lft_id, &rgt_id;
 
-      // real_t vectors copied in distributed memory case
-      std::set<thrust_device::vector<real_t>*> distmem_real_vctrs;
+      // --- containters with vector pointers to help resize and copy vectors ---
+
+      // vectors copied between distributed memories (MPI, multi_CUDA), these are SD attributes
+      std::set<std::pair<thrust_device::vector<real_t>*, real_t>>         distmem_real_vctrs; // pair of vector and its initial value
+      std::set<thrust_device::vector<n_t>*>                               distmem_n_vctrs;
+      //std::set<thrust_device::vector<thrust_size_t>*>                     distmem_size_vctrs; 
 
+      // vetors that are not in distmem_real_vctrs that need to be resized when the number of SDs changes, these are helper variables
+      std::set<thrust_device::vector<real_t>*>         resize_real_vctrs;
+//      std::set<thrust_device::vector<n_t>*>            resize_n_vctrs;
+      std::set<thrust_device::vector<thrust_size_t>*>  resize_size_vctrs;
 
-      // methods
+      // --- methods ---
 
       // fills u01 with n random real numbers uniformly distributed in range [0,1)
       void rand_u01(thrust_size_t n) { rng.generate_n(u01, n); }
@@ -312,7 +327,7 @@ namespace libcloudphxx
         rng(_opts_init.rng_seed),
         src_stp_ctr(0),
         rlx_stp_ctr(0),
-	bcond(bcond),
+        bcond(bcond),
         n_x_bfr(0),
         n_cell_bfr(0),
         mpi_rank(mpi_rank),
@@ -384,40 +399,71 @@ namespace libcloudphxx
         }
 
         // initializing distmem_real_vctrs - list of real_t vectors with properties of SDs that have to be copied/removed/recycled when a SD is copied/removed/recycled
-        // TODO: add to that list vectors of other types (e.g integer pimpl->n)
         // NOTE: this does not include chemical stuff due to the way chem vctrs are organized! multi_CUDA / MPI does not work with chemistry as of now
-        typedef thrust_device::vector<real_t>* ptr_t;
-        ptr_t arr[] = {&rd3, &rw2, &kpa, &vt};
-        distmem_real_vctrs = std::set<ptr_t>(arr, arr + sizeof(arr) / sizeof(ptr_t) );
+        distmem_real_vctrs.insert({&rd3, detail::no_initial_value});
+        distmem_real_vctrs.insert({&rw2, detail::no_initial_value});
+        distmem_real_vctrs.insert({&kpa, detail::no_initial_value});
 
-        if (opts_init.nx != 0)  distmem_real_vctrs.insert(&x);
-        if (opts_init.ny != 0)  distmem_real_vctrs.insert(&y);
-        if (opts_init.nz != 0)  distmem_real_vctrs.insert(&z);
+        distmem_real_vctrs.insert({&vt,  detail::invalid});
+
+        if (opts_init.nx != 0)  distmem_real_vctrs.insert({&x, detail::no_initial_value});
+        if (opts_init.ny != 0)  distmem_real_vctrs.insert({&y, detail::no_initial_value});
+        if (opts_init.nz != 0)  distmem_real_vctrs.insert({&z, detail::no_initial_value});
 
         if(allow_sstp_cond && opts_init.exact_sstp_cond)
         {
-           distmem_real_vctrs.insert(&sstp_tmp_rv);
-           distmem_real_vctrs.insert(&sstp_tmp_th);
-           distmem_real_vctrs.insert(&sstp_tmp_rh);
+           distmem_real_vctrs.insert({&sstp_tmp_rv, detail::no_initial_value});
+           distmem_real_vctrs.insert({&sstp_tmp_th, detail::no_initial_value});
+           distmem_real_vctrs.insert({&sstp_tmp_rh, detail::no_initial_value});
            // sstp_tmp_p needs to be added if a constant pressure profile is used, but this is only known after init - see particles_init
         }
 
         if(opts_init.turb_adve_switch)
         {
-          if(opts_init.nx != 0) distmem_real_vctrs.insert(&up);
-          if(opts_init.ny != 0) distmem_real_vctrs.insert(&vp);
-          if(opts_init.nz != 0) distmem_real_vctrs.insert(&wp);
+          if(opts_init.nx != 0) distmem_real_vctrs.insert({&up, 0});
+          if(opts_init.ny != 0) distmem_real_vctrs.insert({&vp, 0});
+          if(opts_init.nz != 0) distmem_real_vctrs.insert({&wp, 0});
         }
 
         if(opts_init.turb_cond_switch)
         {
-          distmem_real_vctrs.insert(&wp);
-          distmem_real_vctrs.insert(&ssp);
-          distmem_real_vctrs.insert(&dot_ssp);
+          distmem_real_vctrs.insert({&wp, 0});
+          distmem_real_vctrs.insert({&ssp, 0});
+          distmem_real_vctrs.insert({&dot_ssp, 0});
         }
 
         if(opts_init.diag_incloud_time)
-          distmem_real_vctrs.insert(&incloud_time);
+          distmem_real_vctrs.insert({&incloud_time, detail::no_initial_value});
+
+        precoal_distance.resize(config.precoal_stats_bins);
+        precoal_distance_count.resize(config.precoal_stats_bins);
+        thrust::fill(precoal_distance.begin(), precoal_distance.end(), real_t(0));
+        thrust::fill(precoal_distance_count.begin(), precoal_distance_count.end(), n_t(0));
+
+        // initializing distmem_n_vctrs - list of n_t vectors with properties of SDs that have to be copied/removed/recycled when a SD is copied/removed/recycled
+        distmem_n_vctrs.insert(&n);
+
+        // real vctrs that need to be resized but do need to be copied in distmem
+        resize_real_vctrs.insert(&tmp_device_real_part);
+        if(opts_init.chem_switch || allow_sstp_cond || n_dims >= 2)
+          resize_real_vctrs.insert(&tmp_device_real_part1);
+        if((allow_sstp_cond && opts_init.exact_sstp_cond) || n_dims==3 || opts_init.turb_cond_switch)
+          resize_real_vctrs.insert(&tmp_device_real_part2);
+        if(allow_sstp_cond && opts_init.exact_sstp_cond)
+        {
+          resize_real_vctrs.insert(&tmp_device_real_part3);
+          resize_real_vctrs.insert(&tmp_device_real_part4);
+          if(const_p)
+            resize_real_vctrs.insert(&tmp_device_real_part5);
+        }
+
+        resize_size_vctrs.insert(&ijk);
+        resize_size_vctrs.insert(&sorted_ijk);
+        resize_size_vctrs.insert(&sorted_id);
+        resize_size_vctrs.insert(&tmp_device_size_part);
+        if (opts_init.nx != 0) resize_size_vctrs.insert(&i);
+        if (opts_init.ny != 0) resize_size_vctrs.insert(&j);
+        if (opts_init.nz != 0) resize_size_vctrs.insert(&k);
       }
 
       void sanity_checks();
@@ -565,7 +611,7 @@ namespace libcloudphxx
       void update_pstate(thrust_device::vector<real_t> &, thrust_device::vector<real_t> &);
       void update_incloud_time(const real_t &dt);
 
-      void coal(const real_t &dt, const bool &turb_coal);
+      void coal(const real_t &dt, const bool &turb_coal, const real_t &time);
 
       void chem_vol_ante();
       void chem_flag_ante();

src/impl/particles_impl_bcnd.ipp

@@ -124,7 +124,7 @@ namespace libcloudphxx
               arg::_1 >= opts_init.x1
             ) - rgt_id.begin();
 
-            const int no_of_n_vctrs_copied(int(1));
+            const int no_of_n_vctrs_copied(distmem_n_vctrs.size());
             const int no_of_real_vctrs_copied(distmem_real_vctrs.size());
 
             if(lft_count*no_of_n_vctrs_copied > in_n_bfr.size() || rgt_count*no_of_n_vctrs_copied  > in_n_bfr.size())