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..18687a112 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( + const 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( + const 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..77925840a --- /dev/null +++ b/src/particles_impl_rc_adjust.ipp @@ -0,0 +1,101 @@ +// 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 + real_t operator()(real_t rw2, real_t m_added) + { + // turn into added volume + return pow(pow(rw2, real_t(3./2)) + 3. / 4. / +#if !defined(__NVCC__) + pi() +#else + CUDART_PI +#endif + * m_added / (libcloudphxx::common::moist_air::rho_w() / si::kilograms * si::cubic_metres), + real_t(2./3)); + } + }; + struct divide + { + template + BOOST_GPU_ENABLED + real_t 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); + if(!l2e.count(&rc_adjust)) + init_e2l(_rc_adjust, &rc_adjust); + // save input adjustment to the temp vector + sync(_rc_adjust, rc_adjust); + + hskpng_sort(); + + // count number of particles per cell + thrust::pair< + thrust_device::vector::iterator, + thrust_device::vector::iterator + > np = thrust::reduce_by_key( + sorted_ijk.begin(), sorted_ijk.end(), // input - keys + n.begin(), // input - values + count_ijk.begin(), // output - keys + count_num.begin() // output - values + ); + count_n = np.first - count_ijk.begin(); + + // 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()), + 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 drops + thrust::transform( + 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 9761bad4d..9daea1f0e 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..53ec9cc9a 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -222,5 +222,12 @@ namespace libcloudphxx return ret; } + + template + void particles_t::step_rc_adjust(const arrinfo_t rc_adjust) + { + // 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")