From 8cd4a2980e10bf3f8f657944c81734f770af5203 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Fri, 8 Jan 2016 10:10:35 +0100 Subject: [PATCH 1/2] add the rc adjustent function --- bindings/python/lgrngn.hpp | 12 +++ bindings/python/lib.cpp | 3 + include/libcloudph++/lgrngn/particles.hpp | 11 ++ src/particles.tpp | 1 + src/particles_impl_rc_adjust.ipp | 121 ++++++++++++++++++++++ src/particles_pimpl_ctor.ipp | 1 + src/particles_step.ipp | 10 ++ 7 files changed, 159 insertions(+) create mode 100644 src/particles_impl_rc_adjust.ipp diff --git a/bindings/python/lgrngn.hpp b/bindings/python/lgrngn.hpp index 67ca5fc68..65c524067 100644 --- a/bindings/python/lgrngn.hpp +++ b/bindings/python/lgrngn.hpp @@ -122,6 +122,18 @@ namespace libcloudphxx ); } + // + template + void step_rc_adjust( + lgr::particles_proto_t *arg, + const bp::numeric::array &rc_adjust + ) + { + arg->step_rc_adjust( + np2ai(rc_adjust, sz(*arg)) + ); + } + template const lgr::opts_init_t get_oi( lgr::particles_proto_t *arg diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index cd751df9a..3f60aeb79 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -278,6 +278,9 @@ BOOST_PYTHON_MODULE(libcloudphxx) bp::arg("Cz") = bp::numeric::array(bp::object()), bp::arg("ambient_chem") = bp::dict() )) + .def("step_rc_adjust", &lgrngn::step_rc_adjust, ( + bp::arg("rc_adjust") = bp::numeric::array(bp::object()) + )) .def("step_async", &lgr::particles_proto_t::step_async) .def("diag_sd_conc", &lgr::particles_proto_t::diag_sd_conc) .def("diag_all", &lgr::particles_proto_t::diag_all) diff --git a/include/libcloudph++/lgrngn/particles.hpp b/include/libcloudph++/lgrngn/particles.hpp index 74ce7a153..10aade293 100644 --- a/include/libcloudph++/lgrngn/particles.hpp +++ b/include/libcloudph++/lgrngn/particles.hpp @@ -50,6 +50,13 @@ namespace libcloudphxx return 0; } + // add some cloud water + virtual void step_rc_adjust( + arrinfo_t rc_adj + ) { + assert(false); + } + // method for accessing super-droplet statistics virtual void diag_sd_conc() { assert(false); } virtual void diag_all() { assert(false); } @@ -99,6 +106,10 @@ namespace libcloudphxx const opts_t & ); + void step_rc_adjust( + arrinfo_t + ); + // diagnostic methods void diag_sd_conc(); void diag_dry_rng( diff --git a/src/particles.tpp b/src/particles.tpp index 484792e4b..5c47774e4 100644 --- a/src/particles.tpp +++ b/src/particles.tpp @@ -64,6 +64,7 @@ #include "particles_impl_fill_outbuf.ipp" #include "particles_impl_sync.ipp" +#include "particles_impl_rc_adjust.ipp" #include "particles_impl_adve.ipp" #include "particles_impl_cond.ipp" diff --git a/src/particles_impl_rc_adjust.ipp b/src/particles_impl_rc_adjust.ipp new file mode 100644 index 000000000..a27a86215 --- /dev/null +++ b/src/particles_impl_rc_adjust.ipp @@ -0,0 +1,121 @@ +// vim:filetype=cpp +/** @file + * @copyright University of Warsaw + * @section LICENSE + * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) + */ + +namespace libcloudphxx +{ + namespace lgrngn + { + namespace detail + { + struct adj_rw2 + { + template + BOOST_GPU_ENABLED + bool operator()(real_t rw2, real_t m_added) + { + // turn into added volume + m_added /= libcloudphxx::common::moist_air::rho_w() / si::kilograms * si::cubic_metres; + return pow(pow(rw2, real_t(3./2)) + 3. / 4. / +#if !defined(__NVCC__) + pi() +#else + CUDART_PI +#endif + * m_added, real_t(2./3)); + } + }; + + struct is_greater_than_zero + { + template + BOOST_GPU_ENABLED + bool operator()(real_t x) + { + return x > real_t(0.); + } + }; + struct divide + { + template + BOOST_GPU_ENABLED + bool operator()(real_t x, n_t y) + { + return x / real_t(y); + } + }; + }; + template + void particles_t::impl::rc_adjust( + const arrinfo_t &_rc_adjust + ) + { + // use tmp vector to store adjustment + thrust_device::vector &rc_adjust(tmp_device_real_cell); + + // save input adjustment to the temp vector + sync(_rc_adjust, rc_adjust); + + // alias for filtered SDs with RH>=Sc + thrust_device::vector &n_filtered(tmp_device_real_part); + + // turn n_filtered into activated? = yes/no + thrust::replace_if( + n_filtered.begin(), + n_filtered.end(), + detail::is_greater_than_zero(), + real_t(1) + ); + + hskpng_sort(); + + // count number of activated SDs per cell + thrust::pair< + thrust_device::vector::iterator, + thrust_device::vector::iterator + > n = thrust::reduce_by_key( + sorted_ijk.begin(), sorted_ijk.end(), // input - keys + n_filtered.begin(), // input - values + count_ijk.begin(), // output - keys + count_num.begin() // output - values + ); + count_n = n.first - count_ijk.begin(); + + // divide water adjustment by number of activated SDs + thrust::transform( + thrust::make_permutation_iterator(rc_adjust.begin(), count_ijk.begin()), + thrust::make_permutation_iterator(rc_adjust.begin(), count_ijk.end()), + thrust::make_permutation_iterator(count_num.begin(), count_ijk.begin()), + thrust::make_permutation_iterator(rc_adjust.begin(), count_ijk.begin()), + detail::divide() + ); + + // adjustment is in mass_of_added_water[kg] / mass_of_dry_air[kg] + // multiply it now by mass of dry air in the cell + thrust::transform( + rc_adjust.begin(), rc_adjust.end(), + dv.begin(), + rc_adjust.begin(), + thrust::multiplies() + ); + thrust::transform( + rc_adjust.begin(), rc_adjust.end(), + rhod.begin(), + rc_adjust.begin(), + thrust::multiplies() + ); + + // apply the adjustment - change rw of activated SDs + thrust::transform( + thrust::make_permutation_iterator(rw2.begin(), sorted_ijk.begin()), + thrust::make_permutation_iterator(rw2.begin(), sorted_ijk.end()), + thrust::make_permutation_iterator(rc_adjust.begin(), sorted_ijk.begin()), // mass to be added to each SD + thrust::make_permutation_iterator(rw2.begin(), sorted_ijk.begin()), + detail::adj_rw2() + ); + } + }; +}; diff --git a/src/particles_pimpl_ctor.ipp b/src/particles_pimpl_ctor.ipp index 9761bad4d..09af40f77 100644 --- a/src/particles_pimpl_ctor.ipp +++ b/src/particles_pimpl_ctor.ipp @@ -322,6 +322,7 @@ namespace libcloudphxx void init_sstp(); void init_kernel(); void init_vterm(); + void rc_adjust(const arrinfo_t &); void fill_outbuf(); diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 84b232419..533ff3b8e 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -222,5 +222,15 @@ namespace libcloudphxx return ret; } + + template + void particles_t::step_rc_adjust(arrinfo_t rc_adjust) + { + // select activated SDs, stored in tmp_device_real_part + diag_RH_ge_Sc(); + + // apply the adjustment + pimpl->rc_adjust(rc_adjust); + } }; }; From 8c6b6cb69d74284bc11a98452607945774b793b6 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Fri, 8 Jan 2016 15:27:28 +0100 Subject: [PATCH 2/2] fix rc adjustment, add test for it --- include/libcloudph++/lgrngn/particles.hpp | 4 +- src/particles_impl_rc_adjust.ipp | 52 ++++++----------- src/particles_pimpl_ctor.ipp | 2 +- src/particles_step.ipp | 5 +- tests/python/unit/CMakeLists.txt | 2 +- tests/python/unit/test_rc_adjustment.py | 69 +++++++++++++++++++++++ 6 files changed, 90 insertions(+), 44 deletions(-) create mode 100644 tests/python/unit/test_rc_adjustment.py diff --git a/include/libcloudph++/lgrngn/particles.hpp b/include/libcloudph++/lgrngn/particles.hpp index 10aade293..18687a112 100644 --- a/include/libcloudph++/lgrngn/particles.hpp +++ b/include/libcloudph++/lgrngn/particles.hpp @@ -52,7 +52,7 @@ namespace libcloudphxx // add some cloud water virtual void step_rc_adjust( - arrinfo_t rc_adj + const arrinfo_t rc_adj ) { assert(false); } @@ -107,7 +107,7 @@ namespace libcloudphxx ); void step_rc_adjust( - arrinfo_t + const arrinfo_t ); // diagnostic methods diff --git a/src/particles_impl_rc_adjust.ipp b/src/particles_impl_rc_adjust.ipp index a27a86215..77925840a 100644 --- a/src/particles_impl_rc_adjust.ipp +++ b/src/particles_impl_rc_adjust.ipp @@ -15,34 +15,24 @@ namespace libcloudphxx { template BOOST_GPU_ENABLED - bool operator()(real_t rw2, real_t m_added) + real_t operator()(real_t rw2, real_t m_added) { // turn into added volume - m_added /= libcloudphxx::common::moist_air::rho_w() / si::kilograms * si::cubic_metres; return pow(pow(rw2, real_t(3./2)) + 3. / 4. / #if !defined(__NVCC__) pi() #else CUDART_PI #endif - * m_added, real_t(2./3)); - } - }; - - struct is_greater_than_zero - { - template - BOOST_GPU_ENABLED - bool operator()(real_t x) - { - return x > real_t(0.); + * m_added / (libcloudphxx::common::moist_air::rho_w() / si::kilograms * si::cubic_metres), + real_t(2./3)); } }; struct divide { template BOOST_GPU_ENABLED - bool operator()(real_t x, n_t y) + real_t operator()(real_t x, n_t y) { return x / real_t(y); } @@ -55,36 +45,26 @@ namespace libcloudphxx { // use tmp vector to store adjustment thrust_device::vector &rc_adjust(tmp_device_real_cell); - + if(!l2e.count(&rc_adjust)) + init_e2l(_rc_adjust, &rc_adjust); // save input adjustment to the temp vector sync(_rc_adjust, rc_adjust); - // alias for filtered SDs with RH>=Sc - thrust_device::vector &n_filtered(tmp_device_real_part); - - // turn n_filtered into activated? = yes/no - thrust::replace_if( - n_filtered.begin(), - n_filtered.end(), - detail::is_greater_than_zero(), - real_t(1) - ); - hskpng_sort(); - // count number of activated SDs per cell + // count number of particles per cell thrust::pair< thrust_device::vector::iterator, thrust_device::vector::iterator - > n = thrust::reduce_by_key( + > np = thrust::reduce_by_key( sorted_ijk.begin(), sorted_ijk.end(), // input - keys - n_filtered.begin(), // input - values + n.begin(), // input - values count_ijk.begin(), // output - keys count_num.begin() // output - values ); - count_n = n.first - count_ijk.begin(); + count_n = np.first - count_ijk.begin(); - // divide water adjustment by number of activated SDs + // divide water adjustment by number of drops thrust::transform( thrust::make_permutation_iterator(rc_adjust.begin(), count_ijk.begin()), thrust::make_permutation_iterator(rc_adjust.begin(), count_ijk.end()), @@ -108,12 +88,12 @@ namespace libcloudphxx thrust::multiplies() ); - // apply the adjustment - change rw of activated SDs + // apply the adjustment - change rw of drops thrust::transform( - thrust::make_permutation_iterator(rw2.begin(), sorted_ijk.begin()), - thrust::make_permutation_iterator(rw2.begin(), sorted_ijk.end()), - thrust::make_permutation_iterator(rc_adjust.begin(), sorted_ijk.begin()), // mass to be added to each SD - thrust::make_permutation_iterator(rw2.begin(), sorted_ijk.begin()), + rw2.begin(), + rw2.end(), + thrust::make_permutation_iterator(rc_adjust.begin(), ijk.begin()), // mass to be added to each drop + rw2.begin(), detail::adj_rw2() ); } diff --git a/src/particles_pimpl_ctor.ipp b/src/particles_pimpl_ctor.ipp index 09af40f77..9daea1f0e 100644 --- a/src/particles_pimpl_ctor.ipp +++ b/src/particles_pimpl_ctor.ipp @@ -322,7 +322,7 @@ namespace libcloudphxx void init_sstp(); void init_kernel(); void init_vterm(); - void rc_adjust(const arrinfo_t &); + void rc_adjust(const arrinfo_t&); void fill_outbuf(); diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 533ff3b8e..53ec9cc9a 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -224,11 +224,8 @@ namespace libcloudphxx } template - void particles_t::step_rc_adjust(arrinfo_t rc_adjust) + void particles_t::step_rc_adjust(const arrinfo_t rc_adjust) { - // select activated SDs, stored in tmp_device_real_part - diag_RH_ge_Sc(); - // apply the adjustment pimpl->rc_adjust(rc_adjust); } diff --git a/tests/python/unit/CMakeLists.txt b/tests/python/unit/CMakeLists.txt index ea10879f6..d1418aac1 100644 --- a/tests/python/unit/CMakeLists.txt +++ b/tests/python/unit/CMakeLists.txt @@ -1,5 +1,5 @@ # non-pytest tests -foreach(test api_blk_1m api_blk_2m api_lgrngn api_common segfault_20150216 test_col_kernels test_terminal_velocities test_SD_removal test_uniform_init test_source test_chem_coal) +foreach(test api_blk_1m api_blk_2m api_lgrngn api_common segfault_20150216 test_col_kernels test_terminal_velocities test_SD_removal test_uniform_init test_source test_chem_coal test_rc_adjustment) #TODO: indicate that tests depend on the lib add_test( NAME ${test} diff --git a/tests/python/unit/test_rc_adjustment.py b/tests/python/unit/test_rc_adjustment.py new file mode 100644 index 000000000..2b2c106b3 --- /dev/null +++ b/tests/python/unit/test_rc_adjustment.py @@ -0,0 +1,69 @@ +#test if after initialization we have approximately the same water content in each cell + +import sys +sys.path.insert(0, "../../bindings/python/") + +from libcloudphxx import lgrngn +import numpy as np +import time +from math import pi + +# initial exponential distribution in droplet volume +# as a function of ln(r) +def expvolumelnr(lnr): + r=np.exp(lnr) + return n_zero * 3.*np.power(r,3)/np.power(r_zero,3)*np.exp(- np.power((r/r_zero),3)); + +#initial conditions, ca. 2g / m^3 +r_zero = 15e-6 +n_zero = 1.42e8 + +opts_init = lgrngn.opts_init_t() +opts_init.coal_switch=0 +opts_init.sedi_switch=0 +opts_init.cond_switch=0 +opts_init.adve_switch=0 +opts_init.dt = 1 +opts_init.dx = 1 +opts_init.dz = 3 +opts_init.dy = 2 +opts_init.nx = 2 +opts_init.nz = 3 +opts_init.ny = 4 +opts_init.x1 = opts_init.dx * opts_init.nx +opts_init.z1 = opts_init.dz * opts_init.nz +opts_init.y1 = opts_init.dy * opts_init.ny +opts_init.rng_seed = int(time.time()) + +th = 300 * np.ones((opts_init.nx, opts_init.ny, opts_init.nz)) +rv = 0.01 * np.ones((opts_init.nx, opts_init.ny, opts_init.nz)) +rc_adj = 0.001 * np.ones((opts_init.nx, opts_init.ny, opts_init.nz)) +rhod = 1. * np.ones((opts_init.nx, opts_init.ny, opts_init.nz)) + .1 * np.mgrid[1:1+opts_init.nx, 1:1+opts_init.ny, 1:1+opts_init.nz][1] # different densities, hence different water content + +kappa = 1e-6 + +opts_init.dry_distros = {kappa:expvolumelnr} + +opts_init.sd_conc = 64 +opts_init.n_sd_max = opts_init.sd_conc * opts_init.nx * opts_init.ny * opts_init.nz + +try: + prtcls = lgrngn.factory(lgrngn.backend_t.OpenMP, opts_init) +except: + prtcls = lgrngn.factory(lgrngn.backend_t.serial, opts_init) + +prtcls.init(th, rv, rhod) + +prtcls.diag_all() +prtcls.diag_wet_mom(3) # gives specific moment (divided by rhod) +water_pre =1e3* 4/3*pi* np.frombuffer(prtcls.outbuf())#.mean() # dropping a constant +print water_pre + +prtcls.step_rc_adjust(rc_adj) + +prtcls.diag_wet_mom(3) # gives specific moment (divided by rhod) +water_post = 1e3*4/3*pi*np.frombuffer(prtcls.outbuf())#.mean() # dropping a constant +water_post -= rc_adj[0][0][0] * np.ones((opts_init.nx * opts_init.ny * opts_init.nz)) # +print water_post +if(not np.allclose(water_pre, water_post,rtol=0,atol=1e-8)): + raise Exception("Error in adding rc adjustment")