diag precoal distance

src/detail/user_params.hpp

@@ -7,8 +7,8 @@
 // note: description and default values are in uwlcm.cpp, all parameters have to be handled there
 struct user_params_t
 {
-  int nt, outfreq, spinup, rng_seed, rng_seed_init;
-  setup::real_t X, Y, Z, dt;
+  int nt, outfreq, spinup, rng_seed, rng_seed_init, precoal_start, precoal_outfreq;
+  setup::real_t X, Y, Z, dt, precoal_tmax;
   std::string outdir, model_case;
   setup::real_t sgs_delta;
   quantity<si::length, setup::real_t> mean_rd1, mean_rd2;		

src/run_hlpr.cpp

@@ -31,6 +31,7 @@
   #include "solvers/lgrngn/hook_ante_delayed_step_lgrngn.hpp"
   #include "solvers/lgrngn/hook_ante_loop_lgrngn.hpp"
   #include "solvers/lgrngn/hook_ante_step_lgrngn.hpp" 
+  #include "solvers/lgrngn/hook_post_step_lgrngn.hpp" 
   #include "solvers/lgrngn/hook_mixed_rhs_ante_step_lgrngn.hpp"
 #endif
 

src/solvers/lgrngn/hook_ante_delayed_step_lgrngn.hpp

@@ -74,18 +74,20 @@ void slvr_lgrngn<ct_params_t>::hook_ante_delayed_step()
           ftr = async_launcher(
             &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
           );
         else if(params.backend == multi_CUDA)
           ftr = async_launcher(
             &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
           );
         assert(ftr.valid());
       } else 
 #endif
-        prtcls->step_async(params.cloudph_opts);
+        prtcls->step_async(params.cloudph_opts, this->timestep * params.dt);
 
 #if defined(UWLCM_TIMING)
       tend = setup::clock::now();

src/solvers/lgrngn/hook_ante_loop_lgrngn.hpp

@@ -98,6 +98,10 @@ void slvr_lgrngn<ct_params_t>::hook_ante_loop(int nt)
     int n_sd_from_rlx_dry_distros = params.cloudph_opts_init.rlx_sd_per_bin * params.cloudph_opts_init.rlx_bins * params.cloudph_opts_init.nz * 100; // room for 100 rounds of full relaxation... 
     params.cloudph_opts_init.n_sd_max += n_sd_from_rlx_dry_distros;
 
+    // pre-collision statistics diagnosis stuff
+    params.cloudph_opts_init.precoal_stats_tmax = params.user_params.precoal_tmax > 0 ? params.user_params.precoal_tmax : 
+      (params.user_params.nt - params.user_params.precoal_start) * params.user_params.dt;
+
     prtcls.reset(libcloudphxx::lgrngn::factory<real_t>(
       (libcloudphxx::lgrngn::backend_t)params.backend, 
       params.cloudph_opts_init
@@ -196,6 +200,10 @@ void slvr_lgrngn<ct_params_t>::hook_ante_loop(int nt)
     this->record_aux_const("rd_min", "lgrngn", params.cloudph_opts_init.rd_min);  
     this->record_aux_const("rd_max", "lgrngn", params.cloudph_opts_init.rd_max);  
     this->record_aux_const("relax_ccn", "user_params", params.user_params.relax_ccn);  
+    this->record_aux_const("precoal_start", "user_params", params.user_params.precoal_start);  
+    this->record_aux_const("precoal_outfreq", "user_params", params.user_params.precoal_outfreq);  
+    this->record_aux_const("precoal_tmax", "user_params", params.user_params.precoal_tmax);  
+    this->record_aux_const("precoal_stats_tmax", "lgrngn", params.cloudph_opts_init.precoal_stats_tmax);  
   }
   this->mem->barrier();
 }

src/solvers/lgrngn/hook_post_step_lgrngn.hpp

@@ -0,0 +1,34 @@
+#pragma once
+#include "../slvr_lgrngn.hpp"
+
+template <class ct_params_t>
+void slvr_lgrngn<ct_params_t>::hook_post_step()
+{
+  parent_t::hook_post_step(); 
+
+  // storej positions for precoal diag
+  if (this->rank == 0)
+  {
+    if(params.user_params.precoal_start >= 0 && (this->timestep >= params.user_params.precoal_start) && (int(this->timestep - params.user_params.precoal_start) % int(params.user_params.precoal_outfreq) == 0))
+      prtcls->store_ijk(this->timestep * params.user_params.dt);
+  }
+  this->mem->barrier();
+
+  // write mean precoal distance at the end of the simulation to a file
+  // TODO: won't work with MPI
+  if (this->rank == 0)
+  {
+    if(params.user_params.precoal_start >= 0 && this->timestep == params.user_params.nt) 
+    {
+      std::ofstream os(this->outdir+"/precoal_mean_time.dat");
+      os << "# time before collision [s] | mean distance between droplets [m]" << std::endl;
+      std::vector<real_t> precoal_distance = prtcls->diag_precoal_distance();
+      for(int i=0; i<precoal_distance.size(); ++i)
+      {
+        os << (i+0.5) * (params.cloudph_opts_init.precoal_stats_tmax / 100.) << " " << precoal_distance.at(i) << std::endl; // assuming libcloud config.precoal_stats_bins=100
+      }
+    //  os << std::endl;
+    }
+  }
+  this->mem->barrier();
+}

src/solvers/slvr_lgrngn.hpp

@@ -123,6 +123,7 @@ class slvr_lgrngn : public std::conditional_t<ct_params_t::sgs_scheme == libmpda
 
   void hook_ante_loop(int nt);
   void hook_ante_step();
+  void hook_post_step();
   void hook_ante_delayed_step();
   void hook_mixed_rhs_ante_step();
 

src/uwlcm.cpp

@@ -25,6 +25,7 @@
   #include "solvers/lgrngn/hook_ante_delayed_step_lgrngn.hpp"
   #include "solvers/lgrngn/hook_ante_loop_lgrngn.hpp"
   #include "solvers/lgrngn/hook_ante_step_lgrngn.hpp" 
+  #include "solvers/lgrngn/hook_post_step_lgrngn.hpp" 
   #include "solvers/lgrngn/hook_mixed_rhs_ante_step_lgrngn.hpp"
 #endif
 
@@ -91,6 +92,9 @@ int main(int argc, char** argv)
       ("sgs_delta", po::value<setup::real_t>()->default_value(-1) , "subgrid-scale turbulence model length scale [m]. If negative, sgs_delta = dz")
       ("help", "produce a help message (see also --micro X --help)")
       ("relax_th_rv", po::value<bool>()->default_value(false) , "relax per-level mean theta and rv to a desired (case-specific) profile")
+      ("precoal_start", po::value<int>()->default_value(-1) , "timestep at which storing of SD positions (for precoal distance diag) starts (negative for not storing at all)")
+      ("precoal_outfreq", po::value<int>()->default_value(0) , "timestep interval at which SD positions are stored (for precoal distance diag)")
+      ("precoal_tmax", po::value<setup::real_t>()->default_value(-1) , "pre-collision statistics are diagnosed up to precoal_stats_tmax [s] before collision (negative for precoal_tmax=(nt - precoal_start)*dt)")
 
       // aerosol distribution params
       // default values are realistic params, except n1_stp=n2_stp=-1
@@ -138,6 +142,9 @@ int main(int argc, char** argv)
 
     user_params.nt = vm["nt"].as<int>(),
     user_params.spinup = vm["spinup"].as<int>();
+    user_params.precoal_start = vm["precoal_start"].as<int>();
+    user_params.precoal_outfreq = vm["precoal_outfreq"].as<int>();
+    user_params.precoal_tmax = vm["precoal_tmax"].as<setup::real_t>();
 
     // handling rng_seed
     user_params.rng_seed = vm["rng_seed"].as<int>();