[do not merge] add an option to pass RH and T directly

bindings/python/lgrngn.hpp

@@ -111,6 +111,8 @@ namespace libcloudphxx
 	const bp_array &Cx,
 	const bp_array &Cy,
 	const bp_array &Cz,
+	const bp_array &RH,
+	const bp_array &T,
         bp::dict &ambient_chem
       )
       {
@@ -134,6 +136,8 @@ namespace libcloudphxx
 	  np2ai<real_t>(Cx, sz(*arg)),
 	  np2ai<real_t>(Cy, sz(*arg)),
 	  np2ai<real_t>(Cz, sz(*arg)),
+	  np2ai<real_t>(RH, sz(*arg)),
+	  np2ai<real_t>(T, sz(*arg)),
           map
 	);
       }

bindings/python/lib.cpp

@@ -266,6 +266,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("adve_scheme", &lgr::opts_init_t<real_t>::adve_scheme)
       .def_readwrite("sd_conc", &lgr::opts_init_t<real_t>::sd_conc)
       .def_readwrite("sd_conc_large_tail", &lgr::opts_init_t<real_t>::sd_conc_large_tail)
+      .def_readwrite("aerosol_independent_of_rhod", &lgr::opts_init_t<real_t>::aerosol_independent_of_rhod)
       .def_readwrite("sd_const_multi", &lgr::opts_init_t<real_t>::sd_const_multi)
       .def_readwrite("sd_const_multi_dry_sizes", &lgr::opts_init_t<real_t>::sd_const_multi_dry_sizes)
       .def_readwrite("src_sd_conc", &lgr::opts_init_t<real_t>::src_sd_conc)
@@ -294,6 +295,8 @@ BOOST_PYTHON_MODULE(libcloudphxx)
         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("RH")  = BP_ARR_FROM_BP_OBJ,
+        bp::arg("T")  = BP_ARR_FROM_BP_OBJ,
         bp::arg("ambient_chem") = bp::dict()
       ))
       .def("step_async",   &lgr::particles_proto_t<real_t>::step_async)

include/libcloudph++/lgrngn/opts_init.hpp

@@ -28,7 +28,7 @@ namespace libcloudphxx
       // uses shared_ptr to make opts_init copyable
       typedef std::unordered_map<
         real_t,                // kappa
-        std::shared_ptr<unary_function<real_t>> // n(ln(rd)) @ STP 
+        std::shared_ptr<unary_function<real_t>> // n(ln(rd)) @ STP; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true
       > dry_distros_t;
       dry_distros_t dry_distros;
 
@@ -58,6 +58,9 @@ namespace libcloudphxx
       // should more SDs be added to better represent large tail of the distribution
       bool sd_conc_large_tail;
 
+      // should aerosol concentration init be independent of rhod (assumed to be in cm^{-3} and not at STP)
+      bool aerosol_independent_of_rhod;
+
       // or, alternatively to sd_conc_mean, multiplicity of all SDs = const
       int sd_const_multi;
 
@@ -121,6 +124,7 @@ namespace libcloudphxx
         x1(1), y1(1), z1(1),
         sd_conc(0), 
         sd_conc_large_tail(false), 
+        aerosol_independent_of_rhod(false), 
         sd_const_multi(0),
         sd_const_multi_dry_sizes(0),
         dt(0),   

include/libcloudph++/lgrngn/particles.hpp

@@ -38,6 +38,8 @@ namespace libcloudphxx
         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>(),
+        const arrinfo_t<real_t> RH = arrinfo_t<real_t>(),
+        const arrinfo_t<real_t> T = 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); 
@@ -102,6 +104,8 @@ namespace libcloudphxx
         const arrinfo_t<real_t> courant_x,
         const arrinfo_t<real_t> courant_y,
         const arrinfo_t<real_t> courant_z,
+        const arrinfo_t<real_t> RH,
+        const arrinfo_t<real_t> T,
         std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem
       );
 
@@ -176,6 +180,8 @@ namespace libcloudphxx
         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>(),
+        const arrinfo_t<real_t> RH = arrinfo_t<real_t>(),
+        const arrinfo_t<real_t> T = 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> >()
       );
       void step_async(

src/impl/particles_impl.ipp

@@ -142,6 +142,10 @@ namespace libcloudphxx
       // is it a pure const_multi run, i.e. no sd_conc
       bool pure_const_multi;
 
+      // true if RH/T is passed directly to step_sync
+      bool external_RH;
+      bool external_T;
+
       // timestep counter
       n_t stp_ctr;
 
@@ -247,6 +251,8 @@ namespace libcloudphxx
         zero(0),
         n_part(0),
         sorted(false), 
+        external_RH(false), 
+        external_T(false), 
         u01(tmp_device_real_part),
         n_user_params(opts_init.kernel_parameters.size()),
         un(tmp_device_n_part),

src/impl/particles_impl_hskpng_Tpr.ipp

@@ -79,12 +79,13 @@ namespace libcloudphxx
     void particles_t<real_t, device>::impl::hskpng_Tpr()
     {   
       // T  = common::theta_dry::T<real_t>(th, rhod);
-      thrust::transform(
-        th.begin(), th.end(),      // input - first arg
-        rhod.begin(),              // input - second arg
-        T.begin(),                 // output
-        detail::common__theta_dry__T<real_t>() 
-      );
+      if(!external_T)
+        thrust::transform(
+          th.begin(), th.end(),      // input - first arg
+          rhod.begin(),              // input - second arg
+          T.begin(),                 // output
+          detail::common__theta_dry__T<real_t>() 
+        );
 
       {
         typedef thrust::zip_iterator<
@@ -104,12 +105,13 @@ namespace libcloudphxx
         );
 
         // RH = p_v / p_vs = rhod * rv * R_v * T / p_vs
-        thrust::transform(
-          zip_it_t(thrust::make_tuple(rhod.begin(), rv.begin(), T.begin())),  // input - begin
-          zip_it_t(thrust::make_tuple(rhod.end(),   rv.end(),   T.end()  )),  // input - end
-          RH.begin(),                                                         // output
-          detail::RH<real_t>()
-        );
+        if(!external_RH)
+          thrust::transform(
+            zip_it_t(thrust::make_tuple(rhod.begin(), rv.begin(), T.begin())),  // input - begin
+            zip_it_t(thrust::make_tuple(rhod.end(),   rv.end(),   T.end()  )),  // input - end
+            RH.begin(),                                                         // output
+            detail::RH<real_t>()
+          );
       }
 
       // dynamic viscosity

src/impl/particles_impl_init_hskpng_ncell.ipp

@@ -13,9 +13,9 @@ namespace libcloudphxx
     void particles_t<real_t, device>::impl::init_hskpng_ncell()
     {
       // memory allocation
-      T.resize(n_cell);
+     // T.resize(n_cell);
       p.resize(n_cell);
-      RH.resize(n_cell); 
+     // RH.resize(n_cell); 
       eta.resize(n_cell); 
 
       count_ijk.resize(n_cell);

src/impl/particles_impl_init_n.ipp

@@ -80,15 +80,16 @@ namespace libcloudphxx
         using common::earth::rho_stp;
 
         // correcting STP -> actual ambient conditions
-        thrust::transform(
-          tmp_real.begin(), tmp_real.end(),            // input - 1st arg
-          thrust::make_permutation_iterator( // input - 2nd arg
-            tmp_rhod.begin(), 
-            tmp_ijk.begin()
-          ),
-          tmp_real.begin(),                       // output
-          arg::_1 * arg::_2 / real_t(rho_stp<real_t>() / si::kilograms * si::cubic_metres)
-        );
+        if(!opts_init.aerosol_independent_of_rhod)
+          thrust::transform(
+            tmp_real.begin(), tmp_real.end(),            // input - 1st arg
+            thrust::make_permutation_iterator( // input - 2nd arg
+              tmp_rhod.begin(), 
+              tmp_ijk.begin()
+            ),
+            tmp_real.begin(),                       // output
+            arg::_1 * arg::_2 / real_t(rho_stp<real_t>() / si::kilograms * si::cubic_metres)
+          );
 
 
         // adjust to cell volume

src/impl/particles_impl_init_sync.ipp

@@ -16,6 +16,8 @@ namespace libcloudphxx
       rhod.resize(n_cell);
       th.resize(n_cell);
       rv.resize(n_cell);
+      RH.resize(n_cell); 
+      T.resize(n_cell); 
       for (int i = 0; i < chem_gas_n; ++i)
         ambient_chem[(chem_species_t)i].resize(n_cell);
 

src/particles_multi_gpu_step.ipp

@@ -22,10 +22,12 @@ namespace libcloudphxx
       const arrinfo_t<real_t> courant_1,
       const arrinfo_t<real_t> courant_2,
       const arrinfo_t<real_t> courant_3,
+      const arrinfo_t<real_t> RH,
+      const arrinfo_t<real_t> T,
       std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem
     )
     {
-      pimpl->mcuda_run(&particles_t<real_t, CUDA>::step_sync, opts, th, rv, rhod, courant_1, courant_2, courant_3, ambient_chem);
+      pimpl->mcuda_run(&particles_t<real_t, CUDA>::step_sync, opts, th, rv, rhod, courant_1, courant_2, courant_3, RH, T, ambient_chem);
     }
 
     template <typename real_t>

src/particles_step.ipp

@@ -20,6 +20,8 @@ namespace libcloudphxx
       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)
+      const arrinfo_t<real_t> RH,        // optionally, directly feed RH (e.g. KiD-A tests)
+      const arrinfo_t<real_t> T,         // optionally, directly feed T (e.g. KiD-A tests)
       std::map<enum chem_species_t, arrinfo_t<real_t> > ambient_chem
     )
     {
@@ -67,24 +69,36 @@ 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);
       }
 
+      if (!RH.is_null()) pimpl->init_e2l(RH, &pimpl->RH);
+      if (!T.is_null()) pimpl->init_e2l(T, &pimpl->T);
+
       // syncing in Eulerian fields (if not null)
       pimpl->sync(th,             pimpl->th);
       pimpl->sync(rv,             pimpl->rv);
       pimpl->sync(courant_x,      pimpl->courant_x);
       pimpl->sync(courant_y,      pimpl->courant_y);
       pimpl->sync(courant_z,      pimpl->courant_z);
       pimpl->sync(rhod,           pimpl->rhod);
+      pimpl->sync(RH,             pimpl->RH);
+      pimpl->sync(T,              pimpl->T);
+
+      pimpl->external_RH = RH.is_null() ? false : true;
+      pimpl->external_T = T.is_null() ? false : true;
 
       nancheck(pimpl->th, " th after sync-in");
       nancheck(pimpl->rv, " rv after sync-in");
       nancheck(pimpl->courant_x, " courant_x after sync-in");
       nancheck(pimpl->courant_y, " courant_y after sync-in");
       nancheck(pimpl->courant_z, " courant_z after sync-in");
       nancheck(pimpl->rhod, " rhod after sync-in");
+      nancheck(pimpl->RH, " RH after sync-in");
+      nancheck(pimpl->T, " T after sync-in");
 
       assert(*thrust::min_element(pimpl->rv.begin(), pimpl->rv.end()) >= 0);
       assert(*thrust::min_element(pimpl->th.begin(), pimpl->th.end()) >= 0);
       assert(*thrust::min_element(pimpl->rhod.begin(), pimpl->rhod.end()) >= 0);
+      assert(*thrust::min_element(pimpl->RH.begin(), pimpl->RH.end()) >= 0);
+      assert(*thrust::min_element(pimpl->T.begin(), pimpl->T.end()) >= 0);
 
       // check if courants are greater than 1 since it would break the predictor-corrector (halo of size 1 in the x direction) 
       assert(pimpl->opts_init.adve_scheme != as_t::pred_corr || (courant_x.is_null() || ((*(thrust::min_element(pimpl->courant_x.begin(), pimpl->courant_x.end()))) >= real_t(-1.) )) );