Merge pull request #327 from pdziekan/delayed_advection

Delayed advection

bindings/python/lgrngn.hpp

@@ -140,6 +140,72 @@ namespace libcloudphxx
 	);
       }
 
+      // 
+      template <typename real_t>
+      void sync_in(
+	lgr::particles_proto_t<real_t> *arg,
+	const bp_array &th,
+	const bp_array &rv,
+	const bp_array &rhod,
+	const bp_array &Cx,
+	const bp_array &Cy,
+	const bp_array &Cz,
+        bp::dict &ambient_chem
+      )
+      {
+        typedef std::map<enum lgr::chem_species_t, lgr::arrinfo_t<real_t> > map_t;
+        map_t map;
+
+        for (int i = 0; i < len(ambient_chem.keys()); ++i)
+          map.insert(typename map_t::value_type(
+            bp::extract<enum lgr::chem_species_t>(ambient_chem.keys()[i]),
+            np2ai<real_t>(bp::extract<bp_array>(ambient_chem.values()[i]), sz(*arg))
+          ));
+
+	lgr::arrinfo_t<real_t>
+	  np2ai_th(np2ai<real_t>(th, sz(*arg))),
+	  np2ai_rv(np2ai<real_t>(rv, sz(*arg)));
+	arg->sync_in(
+	  np2ai_th,
+	  np2ai_rv,
+	  np2ai<real_t>(rhod, sz(*arg)),
+	  np2ai<real_t>(Cx, sz(*arg)),
+	  np2ai<real_t>(Cy, sz(*arg)),
+	  np2ai<real_t>(Cz, sz(*arg)),
+          map
+	);
+      }
+
+      // 
+      template <typename real_t>
+      void step_cond(
+	lgr::particles_proto_t<real_t> *arg,
+	const lgr::opts_t<real_t> &opts,
+	const bp_array &th,
+	const bp_array &rv,
+        bp::dict &ambient_chem
+      )
+      {
+        typedef std::map<enum lgr::chem_species_t, lgr::arrinfo_t<real_t> > map_t;
+        map_t map;
+
+        for (int i = 0; i < len(ambient_chem.keys()); ++i)
+          map.insert(typename map_t::value_type(
+            bp::extract<enum lgr::chem_species_t>(ambient_chem.keys()[i]),
+            np2ai<real_t>(bp::extract<bp_array>(ambient_chem.values()[i]), sz(*arg))
+          ));
+
+	lgr::arrinfo_t<real_t>
+	  np2ai_th(np2ai<real_t>(th, sz(*arg))),
+	  np2ai_rv(np2ai<real_t>(rv, sz(*arg)));
+	arg->step_cond(
+	  opts, 
+	  np2ai_th,
+	  np2ai_rv,
+          map
+	);
+      }
+
       template <typename real_t>
       bp::dict diag_puddle(lgr::particles_proto_t<real_t> *arg)
       {

bindings/python/lib.cpp

@@ -311,6 +311,20 @@ BOOST_PYTHON_MODULE(libcloudphxx)
         bp::arg("Cz")  = BP_ARR_FROM_BP_OBJ,
         bp::arg("ambient_chem") = bp::dict()
       ))
+      .def("sync_in",    &lgrngn::sync_in<real_t>, (
+        bp::arg("th")  = BP_ARR_FROM_BP_OBJ,
+        bp::arg("rv")  = BP_ARR_FROM_BP_OBJ,
+        bp::arg("rhod")= BP_ARR_FROM_BP_OBJ,
+        bp::arg("Cx")  = BP_ARR_FROM_BP_OBJ,
+        bp::arg("Cy")  = BP_ARR_FROM_BP_OBJ,
+        bp::arg("Cz")  = BP_ARR_FROM_BP_OBJ,
+        bp::arg("ambient_chem") = bp::dict()
+      ))
+      .def("step_cond",    &lgrngn::step_cond<real_t>, (
+        bp::arg("th")  = BP_ARR_FROM_BP_OBJ,
+        bp::arg("rv")  = BP_ARR_FROM_BP_OBJ,
+        bp::arg("ambient_chem") = bp::dict()
+      ))
       .def("step_async",   &lgr::particles_proto_t<real_t>::step_async)
       .def("diag_sd_conc", &lgr::particles_proto_t<real_t>::diag_sd_conc)
       .def("diag_all",     &lgr::particles_proto_t<real_t>::diag_all)

include/libcloudph++/lgrngn/particles.hpp

@@ -44,6 +44,27 @@ namespace libcloudphxx
         assert(false); 
       }  
 
+      virtual void sync_in(
+        arrinfo_t<real_t> th,
+        arrinfo_t<real_t> rv,
+        const arrinfo_t<real_t> rhod      = arrinfo_t<real_t>(),
+        const arrinfo_t<real_t> courant_x = arrinfo_t<real_t>(),
+        const arrinfo_t<real_t> courant_y = arrinfo_t<real_t>(),
+        const arrinfo_t<real_t> courant_z = arrinfo_t<real_t>(),
+        std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem = std::map<enum chem_species_t, arrinfo_t<real_t> >()
+      ) { 
+        assert(false); 
+      }  
+
+      virtual void step_cond(
+        const opts_t<real_t> &,
+        arrinfo_t<real_t> th,
+        arrinfo_t<real_t> rv,
+        std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem = std::map<enum chem_species_t, arrinfo_t<real_t> >()
+      ) { 
+        assert(false); 
+      }  
+
       // returns accumulated rainfall
       virtual void step_async(
         const opts_t<real_t> &
@@ -109,6 +130,23 @@ namespace libcloudphxx
         std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem
       );
 
+      void sync_in(
+        arrinfo_t<real_t> th,
+        arrinfo_t<real_t> rv,
+        const arrinfo_t<real_t> rhod,
+        const arrinfo_t<real_t> courant_x,
+        const arrinfo_t<real_t> courant_y,
+        const arrinfo_t<real_t> courant_z,
+        std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem
+      );
+
+      void step_cond(
+        const opts_t<real_t> &,
+        arrinfo_t<real_t> th,
+        arrinfo_t<real_t> rv,
+        std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem = std::map<enum chem_species_t, arrinfo_t<real_t> >()
+      );
+
       void step_async(
         const opts_t<real_t> &
       );
@@ -185,6 +223,24 @@ namespace libcloudphxx
         const arrinfo_t<real_t> courant_z,
         std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem 
       );
+
+      void sync_in(
+        arrinfo_t<real_t> th,
+        arrinfo_t<real_t> rv,
+        const arrinfo_t<real_t> rhod     ,
+        const arrinfo_t<real_t> courant_x,
+        const arrinfo_t<real_t> courant_y,
+        const arrinfo_t<real_t> courant_z,
+        std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem 
+      );
+
+      void step_cond(
+        const opts_t<real_t> &,
+        arrinfo_t<real_t> th,
+        arrinfo_t<real_t> rv,
+        std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem = std::map<enum chem_species_t, arrinfo_t<real_t> >()
+      );
+
       void step_async(
         const opts_t<real_t> &
       );

src/impl/particles_impl.ipp

@@ -34,7 +34,10 @@ namespace libcloudphxx
       typedef unsigned long long n_t; // thrust_size_t?
 
       // order of operation flags
-      bool init_called, should_now_run_async, selected_before_counting;
+      bool init_called, should_now_run_async, selected_before_counting, should_now_run_cond;
+
+      // did density vary in this step
+      bool var_rho;
 
       // member fields
       opts_init_t<real_t> opts_init; // a copy
@@ -236,6 +239,8 @@ namespace libcloudphxx
         init_called(false),
         should_now_run_async(false),
         selected_before_counting(false),
+        should_now_run_cond(false),
+        var_rho(false),
 	opts_init(_opts_init),
 	n_dims( // 0, 1, 2 or 3
           opts_init.nx/m1(opts_init.nx) + 
@@ -446,10 +451,10 @@ namespace libcloudphxx
 
       void src(const real_t &dt);
 
-      void sstp_step(const int &step, const bool &var_rho);
+      void sstp_step(const int &step);
       void sstp_step_exact(const int &step);
       void sstp_save();
-      void sstp_step_chem(const int &step, const bool &var_rho);
+      void sstp_step_chem(const int &step);
       void sstp_save_chem();
 
       void step_finalize(const opts_t<real_t>&);

src/impl/particles_impl_sstp.ipp

@@ -46,8 +46,7 @@ namespace libcloudphxx
 
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::sstp_step(
-      const int &step,
-      const bool &var_rho // if rho varied and need to be updated
+      const int &step
     )
     {   
       if (opts_init.sstp_cond == 1) return;

src/impl/particles_impl_sstp_chem.ipp

@@ -51,8 +51,7 @@ namespace libcloudphxx
 
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::sstp_step_chem(
-      const int &step,
-      const bool &var_rho // if rho varied and need to be updated
+      const int &step
     )
     {   
       if (opts_init.sstp_chem == 1) return;

src/particles_multi_gpu_step.ipp

@@ -28,6 +28,31 @@ namespace libcloudphxx
       pimpl->mcuda_run(&particles_t<real_t, CUDA>::step_sync, opts, th, rv, rhod, courant_1, courant_2, courant_3, ambient_chem);
     }
 
+    template <typename real_t>
+    void particles_t<real_t, multi_CUDA>::sync_in(
+      arrinfo_t<real_t> th,
+      arrinfo_t<real_t> rv,
+      const arrinfo_t<real_t> rhod,
+      const arrinfo_t<real_t> courant_1,
+      const arrinfo_t<real_t> courant_2,
+      const arrinfo_t<real_t> courant_3,
+      std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem
+    )
+    {
+      pimpl->mcuda_run(&particles_t<real_t, CUDA>::sync_in, th, rv, rhod, courant_1, courant_2, courant_3, ambient_chem);
+    }
+
+    template <typename real_t>
+    void particles_t<real_t, multi_CUDA>::step_cond(
+      const opts_t<real_t> &opts,
+      arrinfo_t<real_t> th,
+      arrinfo_t<real_t> rv,
+      std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem
+    )
+    {
+      pimpl->mcuda_run(&particles_t<real_t, CUDA>::step_cond, opts, th, rv, ambient_chem);
+    }
+
     template <typename real_t>
     void particles_t<real_t, multi_CUDA>::step_async(
       const opts_t<real_t> &opts

src/particles_step.ipp

@@ -22,6 +22,21 @@ namespace libcloudphxx
       const arrinfo_t<real_t> courant_z, // defaults to NULL-NULL pair (e.g. kinematic model)
       std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem
     )
+    {
+      sync_in(th, rv, rhod, courant_x, courant_y, courant_z, ambient_chem);
+      step_cond(opts, th, rv, ambient_chem);
+    }
+
+    template <typename real_t, backend_t device>
+    void particles_t<real_t, device>::sync_in(
+      arrinfo_t<real_t> th,
+      arrinfo_t<real_t> rv,
+      const arrinfo_t<real_t> rhod,      // defaults to NULL-NULL pair (e.g. kinematic or boussinesq model)
+      const arrinfo_t<real_t> courant_x, // defaults to NULL-NULL pair (e.g. kinematic model)
+      const arrinfo_t<real_t> courant_y, // defaults to NULL-NULL pair (e.g. kinematic model)
+      const arrinfo_t<real_t> courant_z, // defaults to NULL-NULL pair (e.g. kinematic model)
+      std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem
+    )
     {
       // sanity checks
       if (!pimpl->init_called)
@@ -67,6 +82,9 @@ namespace libcloudphxx
         if (!courant_z.is_null()) pimpl->init_e2l(courant_z, &pimpl->courant_z, 0, 0, 1, pimpl->n_x_bfr * max(1, pimpl->opts_init.ny) - pimpl->halo_z);
       }
 
+      // did rhod change
+      pimpl->var_rho = !rhod.is_null();
+
       // syncing in Eulerian fields (if not null)
       pimpl->sync(th,             pimpl->th);
       pimpl->sync(rv,             pimpl->rv);
@@ -98,6 +116,23 @@ namespace libcloudphxx
           );
         }
       }
+
+      pimpl->should_now_run_cond = true;
+    }
+
+    template <typename real_t, backend_t device>
+    void particles_t<real_t, device>::step_cond(
+      const opts_t<real_t> &opts,
+      arrinfo_t<real_t> th,                                            // for sync-out
+      arrinfo_t<real_t> rv,                                            // for sync-out
+      std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem   // for sync-out
+    )
+    {
+      //sanity checks
+      if (!pimpl->should_now_run_cond)
+        throw std::runtime_error("please call sync_in() before calling step_cond()");
+
+      pimpl->should_now_run_cond = false;
 
       // condensation/evaporation 
       // if const_p == True, pressure is not substepped
@@ -122,7 +157,7 @@ namespace libcloudphxx
         {
           for (int step = 0; step < pimpl->opts_init.sstp_cond; ++step) 
           {   
-            pimpl->sstp_step(step, !rhod.is_null());
+            pimpl->sstp_step(step);
             pimpl->hskpng_Tpr(); 
             pimpl->cond(pimpl->opts_init.dt / pimpl->opts_init.sstp_cond, opts.RH_max);
           }
@@ -146,7 +181,7 @@ namespace libcloudphxx
           if (opts.chem_dsl)
           {
             //adjust trace gases to substepping
-            pimpl->sstp_step_chem(step, !rhod.is_null());
+            pimpl->sstp_step_chem(step);
 
             //dissolving trace gases (Henrys law)
             pimpl->chem_henry(pimpl->opts_init.dt / pimpl->opts_init.sstp_chem);

tests/python/unit/api_lgrngn.py

@@ -130,6 +130,41 @@ def lognormal(lnr):
 print frombuffer(prtcls.outbuf())
 assert frombuffer(prtcls.outbuf()) == opts_init.sd_conc # parcel set-up
 
+# ----------
+# 0D (parcel) with explicit calls to sync_in and step_cond 
+print "0D"
+rhod = arr_t([  1.])
+th   = arr_t([300.])
+rv   = arr_t([  0.01])
+
+prtcls = lgrngn.factory(backend, opts_init)
+prtcls.init(th, rv, rhod)
+try: 
+  prtcls.init(th, rv, rhod)
+  raise Exception("multiple init call not reported!")
+except:
+  pass
+try:
+  prtcls.step_cond(opts, th, rv)
+  raise Exception("sync_in/cond order mismatch not reported!")
+except:
+  pass
+prtcls.sync_in(th, rv, rhod)
+prtcls.step_cond(opts, th, rv)
+prtcls.step_async(opts)
+prtcls.step_sync(opts, th, rv)
+prtcls.diag_dry_rng(0.,1.)
+prtcls.diag_wet_rng(0.,1.)
+prtcls.diag_dry_mom(1)
+prtcls.diag_wet_mom(1)
+prtcls.diag_kappa_mom(1)
+puddle = prtcls.diag_puddle()
+print 'puddle: ', puddle
+#prtcls.diag_chem(lgrngn.chem_species_t.OH)
+prtcls.diag_all()
+prtcls.diag_sd_conc()
+print frombuffer(prtcls.outbuf())
+assert frombuffer(prtcls.outbuf()) == opts_init.sd_conc # parcel set-up
 
 
 # ----------