diff --git a/bindings/python/lgrngn.hpp b/bindings/python/lgrngn.hpp index bb7ce21bc..655dfdc1f 100644 --- a/bindings/python/lgrngn.hpp +++ b/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(Cx, sz(*arg)), np2ai(Cy, sz(*arg)), np2ai(Cz, sz(*arg)), + np2ai(RH, sz(*arg)), + np2ai(T, sz(*arg)), map ); } diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index 1e010f6a2..5cd230482 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -266,6 +266,7 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("adve_scheme", &lgr::opts_init_t::adve_scheme) .def_readwrite("sd_conc", &lgr::opts_init_t::sd_conc) .def_readwrite("sd_conc_large_tail", &lgr::opts_init_t::sd_conc_large_tail) + .def_readwrite("aerosol_independent_of_rhod", &lgr::opts_init_t::aerosol_independent_of_rhod) .def_readwrite("sd_const_multi", &lgr::opts_init_t::sd_const_multi) .def_readwrite("sd_const_multi_dry_sizes", &lgr::opts_init_t::sd_const_multi_dry_sizes) .def_readwrite("src_sd_conc", &lgr::opts_init_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::step_async) diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index 73d00f9ec..1d1caa937 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/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> // n(ln(rd)) @ STP + std::shared_ptr> // 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), diff --git a/include/libcloudph++/lgrngn/particles.hpp b/include/libcloudph++/lgrngn/particles.hpp index 66bfda42c..c70968eee 100644 --- a/include/libcloudph++/lgrngn/particles.hpp +++ b/include/libcloudph++/lgrngn/particles.hpp @@ -38,6 +38,8 @@ namespace libcloudphxx const arrinfo_t courant_x = arrinfo_t(), const arrinfo_t courant_y = arrinfo_t(), const arrinfo_t courant_z = arrinfo_t(), + const arrinfo_t RH = arrinfo_t(), + const arrinfo_t T = arrinfo_t(), std::map > ambient_chem = std::map >() ) { assert(false); @@ -102,6 +104,8 @@ namespace libcloudphxx const arrinfo_t courant_x, const arrinfo_t courant_y, const arrinfo_t courant_z, + const arrinfo_t RH, + const arrinfo_t T, std::map > ambient_chem ); @@ -176,6 +180,8 @@ namespace libcloudphxx const arrinfo_t courant_x = arrinfo_t(), const arrinfo_t courant_y = arrinfo_t(), const arrinfo_t courant_z = arrinfo_t(), + const arrinfo_t RH = arrinfo_t(), + const arrinfo_t T = arrinfo_t(), std::map > ambient_chem = std::map >() ); void step_async( diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index b290c70ae..f5c3725da 100644 --- a/src/impl/particles_impl.ipp +++ b/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), diff --git a/src/impl/particles_impl_hskpng_Tpr.ipp b/src/impl/particles_impl_hskpng_Tpr.ipp index 476e51f1b..52dac6e42 100644 --- a/src/impl/particles_impl_hskpng_Tpr.ipp +++ b/src/impl/particles_impl_hskpng_Tpr.ipp @@ -79,12 +79,13 @@ namespace libcloudphxx void particles_t::impl::hskpng_Tpr() { // T = common::theta_dry::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() - ); + 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() + ); { 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() - ); + 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() + ); } // dynamic viscosity diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index 2af945aba..99d0ad662 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -13,9 +13,9 @@ namespace libcloudphxx void particles_t::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); diff --git a/src/impl/particles_impl_init_n.ipp b/src/impl/particles_impl_init_n.ipp index d24d8bb5a..b96751dcb 100644 --- a/src/impl/particles_impl_init_n.ipp +++ b/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() / 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() / si::kilograms * si::cubic_metres) + ); // adjust to cell volume diff --git a/src/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index fe51f5b89..d3a7f8481 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/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); diff --git a/src/particles_multi_gpu_step.ipp b/src/particles_multi_gpu_step.ipp index 63f9e4391..b7d5d4e00 100644 --- a/src/particles_multi_gpu_step.ipp +++ b/src/particles_multi_gpu_step.ipp @@ -22,10 +22,12 @@ namespace libcloudphxx const arrinfo_t courant_1, const arrinfo_t courant_2, const arrinfo_t courant_3, + const arrinfo_t RH, + const arrinfo_t T, std::map > ambient_chem ) { - pimpl->mcuda_run(&particles_t::step_sync, opts, th, rv, rhod, courant_1, courant_2, courant_3, ambient_chem); + pimpl->mcuda_run(&particles_t::step_sync, opts, th, rv, rhod, courant_1, courant_2, courant_3, RH, T, ambient_chem); } template diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 132b8cc8d..d28d72b0f 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -20,6 +20,8 @@ namespace libcloudphxx const arrinfo_t courant_x, // defaults to NULL-NULL pair (e.g. kinematic model) const arrinfo_t courant_y, // defaults to NULL-NULL pair (e.g. kinematic model) const arrinfo_t courant_z, // defaults to NULL-NULL pair (e.g. kinematic model) + const arrinfo_t RH, // optionally, directly feed RH (e.g. KiD-A tests) + const arrinfo_t T, // optionally, directly feed T (e.g. KiD-A tests) std::map > ambient_chem ) { @@ -67,6 +69,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); } + 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); @@ -74,6 +79,11 @@ namespace libcloudphxx 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"); @@ -81,10 +91,14 @@ namespace libcloudphxx 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.) )) );