diff --git a/CMakeLists.txt b/CMakeLists.txt index 2bba307c6..250e14c59 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -262,6 +262,8 @@ if (EXISTS "${CMAKE_SOURCE_DIR}/.git") ) endif() +include_directories(${libcloudphxx_INCLUDE_DIRS}) + add_subdirectory(src) enable_testing() add_subdirectory(tests) diff --git a/bindings/python/lgrngn.hpp b/bindings/python/lgrngn.hpp index 2975c9f88..2106b641b 100644 --- a/bindings/python/lgrngn.hpp +++ b/bindings/python/lgrngn.hpp @@ -113,6 +113,7 @@ namespace libcloudphxx const bp_array &Cx, const bp_array &Cy, const bp_array &Cz, + const bp_array &diss_rate, bp::dict &ambient_chem ) { @@ -136,6 +137,7 @@ namespace libcloudphxx np2ai(Cx, sz(*arg)), np2ai(Cy, sz(*arg)), np2ai(Cz, sz(*arg)), + np2ai(diss_rate, sz(*arg)), map ); } @@ -150,6 +152,7 @@ namespace libcloudphxx const bp_array &Cx, const bp_array &Cy, const bp_array &Cz, + const bp_array &diss_rate, bp::dict &ambient_chem ) { @@ -172,6 +175,7 @@ namespace libcloudphxx np2ai(Cx, sz(*arg)), np2ai(Cy, sz(*arg)), np2ai(Cz, sz(*arg)), + np2ai(diss_rate, sz(*arg)), map ); } diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index 8e956c375..bbb3b7510 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -244,6 +244,9 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("chem_dsc", &lgr::opts_t::chem_dsc) .def_readwrite("chem_rct", &lgr::opts_t::chem_rct) .def_readwrite("RH_max", &lgr::opts_t::RH_max) + .def_readwrite("turb_adve", &lgr::opts_t::turb_adve) + .def_readwrite("turb_cond", &lgr::opts_t::turb_cond) + .def_readwrite("turb_coal", &lgr::opts_t::turb_coal) ; bp::class_>("opts_init_t") .add_property("dry_distros", &lgrngn::get_dd, &lgrngn::set_dd) @@ -269,6 +272,9 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("coal_switch", &lgr::opts_init_t::coal_switch) .def_readwrite("sedi_switch", &lgr::opts_init_t::sedi_switch) .def_readwrite("src_switch", &lgr::opts_init_t::src_switch) + .def_readwrite("turb_adve_switch", &lgr::opts_init_t::turb_adve_switch) + .def_readwrite("turb_cond_switch", &lgr::opts_init_t::turb_cond_switch) + .def_readwrite("turb_coal_switch", &lgr::opts_init_t::turb_coal_switch) .def_readwrite("exact_sstp_cond", &lgr::opts_init_t::exact_sstp_cond) .def_readwrite("sstp_cond", &lgr::opts_init_t::sstp_cond) .def_readwrite("sstp_coal", &lgr::opts_init_t::sstp_coal) @@ -308,6 +314,7 @@ 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("diss_rate") = BP_ARR_FROM_BP_OBJ, bp::arg("ambient_chem") = bp::dict() )) .def("sync_in", &lgrngn::sync_in, ( @@ -317,6 +324,7 @@ 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("diss_rate") = BP_ARR_FROM_BP_OBJ, bp::arg("ambient_chem") = bp::dict() )) .def("step_cond", &lgrngn::step_cond, ( diff --git a/include/libcloudph++/common/turbulence.hpp b/include/libcloudph++/common/turbulence.hpp new file mode 100644 index 000000000..3ac1f7748 --- /dev/null +++ b/include/libcloudph++/common/turbulence.hpp @@ -0,0 +1,136 @@ +#pragma once + +namespace libcloudphxx +{ + namespace common + { + // Grabowski and Abade 2017 + namespace turbulence + { + libcloudphxx_const_derived(si::dimensionless, C_E, real_t(0.845)); + libcloudphxx_const_derived(si::dimensionless, C_tau, real_t(1.5)); + libcloudphxx_const_derived(si::dimensionless, cube_root_of_two_pi, pow(real_t(2) * +#if !defined(__NVCC__) + pi() +#else + CUDART_PI +#endif + , real_t(1./3.)) + ); + +#if !defined(__NVCC__) + using std::pow; + using std::sqrt; +#endif + + typedef divide_typeof_helper< + si::dimensionless, + si::length + >::type one_over_length; + libcloudphxx_const(one_over_length, a_1, real_t(3e-4), real_t(1) / si::meters); + + typedef divide_typeof_helper< + si::area, + si::time + >::type area_over_time; + libcloudphxx_const(area_over_time, a_2, real_t(2.8e-4), si::square_meters / si::seconds); + + typedef divide_typeof_helper< + si::area, + multiply_typeof_helper::type + >::type area_over_time_squared; + + typedef divide_typeof_helper< + area_over_time_squared, + si::time + >::type area_over_time_cubed; + + typedef divide_typeof_helper< + si::dimensionless, + si::time + >::type one_over_time; + + typedef divide_typeof_helper< + si::dimensionless, + si::area + >::type one_over_area; + + template + quantity length_scale( + quantity dx + ) + { return quantity(dx);} + + template + quantity length_scale( + quantity dx, + quantity dz + ) + //{ return quantity(sqrt(dx*dz));} + { return quantity(dz);} + + template + quantity length_scale( + quantity dx, + quantity dy, + quantity dz + ) + //{ return quantity(pow((dx/si::metres)*(dy/si::metres)*(dz/si::metres), real_t(1./3.)) * si::metres);} + { return quantity(dz);} + + template + BOOST_GPU_ENABLED + quantity tke( + const quantity &diss_rate, // dissipation rate (epsilon) + const quantity &L // characteristic length-scale + ) + { + return quantity(pow((L * diss_rate) / si::cubic_metres * si::seconds * si::seconds * si::seconds / C_E(), real_t(2./3.)) * si::metres * si::metres / si::seconds / si::seconds); + }; + + template + BOOST_GPU_ENABLED + quantity tau( + const quantity &tke, + const quantity &L // characteristic length-scale + ) + { + return quantity(L / cube_root_of_two_pi() * sqrt(C_tau() / tke)); + }; + + template + BOOST_GPU_ENABLED + quantity update_turb_vel( + const quantity &wp, + const quantity &tau, + const quantity &dt, + const quantity &tke, + const quantity &r_normal + ) + { + quantity exp_m_dt_ov_tau(exp(-dt / tau)); + return quantity(wp * exp_m_dt_ov_tau + sqrt((real_t(1) - exp_m_dt_ov_tau * exp_m_dt_ov_tau) * real_t(2./3.) * tke ) * r_normal); + }; + + template + BOOST_GPU_ENABLED + quantity tau_relax( + const quantity &wet_mom_1_over_vol + ) + { + return quantity(real_t(1) / ( a_2() * wet_mom_1_over_vol ) ); + }; + + template + BOOST_GPU_ENABLED + quantity dot_turb_ss( + const quantity &ssp, + const quantity &wp, + const quantity &tau_rlx + ) + { + return quantity(a_1() * wp - ssp / tau_rlx); + }; + }; + }; +}; diff --git a/include/libcloudph++/lgrngn/opts.hpp b/include/libcloudph++/lgrngn/opts.hpp index 58849bed1..04d67e4a4 100644 --- a/include/libcloudph++/lgrngn/opts.hpp +++ b/include/libcloudph++/lgrngn/opts.hpp @@ -19,7 +19,7 @@ namespace libcloudphxx struct opts_t { // process toggling - bool adve, sedi, cond, coal, src, rcyc; + bool adve, sedi, cond, coal, src, rcyc, turb_adve, turb_cond, turb_coal; // RH limit for drop growth real_t RH_max; @@ -32,6 +32,7 @@ namespace libcloudphxx opts_t() : adve(true), sedi(true), cond(true), coal(true), src(false), rcyc(false), chem_dsl(false), chem_dsc(false), chem_rct(false), + turb_adve(false), turb_cond(false), turb_coal(false), RH_max(44) // :) (anything greater than 1.1 would be enough { } diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index e8da5eee8..b390c2312 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -102,7 +102,10 @@ namespace libcloudphxx coal_switch, // if false no coalescence throughout the whole simulation sedi_switch, // if false no sedimentation throughout the whole simulation src_switch, // if false no source throughout the whole simulation - exact_sstp_cond; // if true, use per-particle sstp_cond logic, if false, use per-cell + exact_sstp_cond, // if true, use per-particle sstp_cond logic, if false, use per-cell + turb_adve_switch, // if true, turbulent motion of SDs is modeled + turb_cond_switch, // if true, turbulent condensation of SDs is modeled + turb_coal_switch; // if true, turbulent coalescence kernels can be used int sstp_chem; real_t chem_rho; @@ -122,6 +125,8 @@ namespace libcloudphxx // large-scale horizontal wind divergence [1/s], used to calculate subsidence rate as -div_LS*z real_t div_LS; + real_t rd_min; // minimal dry radius of droplets (works only for init from spectrum) + // ctor with defaults (C++03 compliant) ... opts_init_t() : nx(0), ny(0), nz(0), @@ -140,6 +145,9 @@ namespace libcloudphxx coal_switch(true), // coalescence turned on by default src_switch(false), // source turned off by default exact_sstp_cond(false), + turb_cond_switch(false), + turb_adve_switch(false), + turb_coal_switch(false), RH_max(.95), // value seggested in Lebo and Seinfeld 2011 chem_rho(0), // dry particle density //TODO add checking if the user gave a different value (np w init) (was 1.8e-3) rng_seed(44), @@ -152,7 +160,8 @@ namespace libcloudphxx n_sd_max(0), src_sd_conc(0), src_z1(0), - div_LS(0.) + div_LS(0.), + rd_min(0.) {} // dtor (just to silence -Winline warnings) diff --git a/include/libcloudph++/lgrngn/particles.hpp b/include/libcloudph++/lgrngn/particles.hpp index 1317b73e3..e4a083d31 100644 --- a/include/libcloudph++/lgrngn/particles.hpp +++ b/include/libcloudph++/lgrngn/particles.hpp @@ -39,6 +39,7 @@ 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 diss_rate = arrinfo_t(), // TKE dissipation rate (epsilon) std::map > ambient_chem = std::map >() ) { assert(false); @@ -51,6 +52,7 @@ 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 diss_rate = arrinfo_t(), std::map > ambient_chem = std::map >() ) { assert(false); @@ -127,6 +129,7 @@ namespace libcloudphxx const arrinfo_t courant_x, const arrinfo_t courant_y, const arrinfo_t courant_z, + const arrinfo_t diss_rate, std::map > ambient_chem ); @@ -137,6 +140,7 @@ namespace libcloudphxx const arrinfo_t courant_x, const arrinfo_t courant_y, const arrinfo_t courant_z, + const arrinfo_t diss_rate, std::map > ambient_chem ); @@ -221,6 +225,7 @@ namespace libcloudphxx const arrinfo_t courant_x, const arrinfo_t courant_y, const arrinfo_t courant_z, + const arrinfo_t diss_rate, std::map > ambient_chem ); @@ -231,6 +236,7 @@ namespace libcloudphxx const arrinfo_t courant_x, const arrinfo_t courant_y, const arrinfo_t courant_z, + const arrinfo_t diss_rate, std::map > ambient_chem ); diff --git a/models/kinematic_2D/src/kin_cloud_2d_lgrngn_chem.hpp b/models/kinematic_2D/src/kin_cloud_2d_lgrngn_chem.hpp index 765173f7a..4004b5cd2 100644 --- a/models/kinematic_2D/src/kin_cloud_2d_lgrngn_chem.hpp +++ b/models/kinematic_2D/src/kin_cloud_2d_lgrngn_chem.hpp @@ -259,6 +259,7 @@ class kin_cloud_2d_lgrngn_chem : public kin_cloud_2d_lgrngn libcloudphxx::lgrngn::arrinfo_t(), libcloudphxx::lgrngn::arrinfo_t(), libcloudphxx::lgrngn::arrinfo_t(), + libcloudphxx::lgrngn::arrinfo_t(), ambient_chem ); diff --git a/src/detail/kernel_onishi_nograv.hpp b/src/detail/kernel_onishi_nograv.hpp index 4de64e9a9..fb8093ccd 100644 --- a/src/detail/kernel_onishi_nograv.hpp +++ b/src/detail/kernel_onishi_nograv.hpp @@ -29,6 +29,8 @@ namespace libcloudphxx BOOST_GPU_ENABLED real_t kernel_onishi_nograv(const real_t &r1, const real_t &r2, const real_t &Re_l, const real_t &eps, real_t dnu, real_t ratio_den) { + if(eps < 1e-10) return 0.; + #if !defined(__NVCC__) using std::max; using std::pow; diff --git a/src/detail/kernels.hpp b/src/detail/kernels.hpp index 553694dea..b231b2f9a 100644 --- a/src/detail/kernels.hpp +++ b/src/detail/kernels.hpp @@ -202,14 +202,14 @@ namespace libcloudphxx //ctor BOOST_GPU_ENABLED - kernel_onishi(thrust_device::pointer k_params, real_t r_max) : kernel_geometric(k_params, 2, r_max) {} + kernel_onishi(thrust_device::pointer k_params, real_t r_max) : kernel_geometric(k_params, 1, r_max) {} kernel_onishi() = default; BOOST_GPU_ENABLED virtual real_t calc(const tpl_calc_wrap &tpl_wrap) const { enum { n_a_ix, n_b_ix, rw2_a_ix, rw2_b_ix, vt_a_ix, vt_b_ix, rd3_a_ix, rd3_b_ix }; - enum { rhod_ix, eta_ix }; + enum { rhod_ix, eta_ix, diss_rate_ix }; #if !defined(__NVCC__) using std::sqrt; @@ -218,7 +218,7 @@ namespace libcloudphxx real_t rwb = sqrt( thrust::get(tpl_wrap.get_rw())); real_t onishi_nograv = detail::kernel_onishi_nograv // value of the turbulent onishi kernel that neglects gravitational settling ( - rwa, rwb, kernel_base::k_params[1], kernel_base::k_params[0], // k_params[0] - epsilon, k_params[1] - Re_lambda + rwa, rwb, kernel_base::k_params[0], thrust::get(tpl_wrap.get_ro_calc()), // k_params[0] - Re_lambda thrust::get(tpl_wrap.get_ro_calc()) / thrust::get(tpl_wrap.get_ro_calc()), // kinetic viscosity common::moist_air::rho_w() / si::kilograms * si::cubic_metres / thrust::get(tpl_wrap.get_ro_calc()) // ratio of water to air density ); diff --git a/src/detail/tpl_calc_wrapper.hpp b/src/detail/tpl_calc_wrapper.hpp index da8d6a305..840e69401 100644 --- a/src/detail/tpl_calc_wrapper.hpp +++ b/src/detail/tpl_calc_wrapper.hpp @@ -1,6 +1,8 @@ #pragma once #include +// TODO: enums used to get from these tuples are defined separately in at least two other places, define them once here + namespace libcloudphxx { namespace lgrngn @@ -21,7 +23,8 @@ namespace libcloudphxx typedef thrust::tuple< real_t, // rhod - real_t // eta + real_t, // eta + real_t // tke dissipation rate > tpl_ro_calc_t; tpl_rw_t tpl_rw; diff --git a/src/detail/urand.hpp b/src/detail/urand.hpp index 7ca11b42a..7fe52b05a 100644 --- a/src/detail/urand.hpp +++ b/src/detail/urand.hpp @@ -23,9 +23,11 @@ namespace libcloudphxx // serial version using C++11's using engine_t = std::mt19937; using dist_u01_t = std::uniform_real_distribution; + using dist_normal01_t = std::normal_distribution; using dist_un_t = std::uniform_int_distribution; engine_t engine; dist_u01_t dist_u01; + dist_normal01_t dist_normal01; dist_un_t dist_un; struct fnctr_u01 @@ -35,6 +37,13 @@ namespace libcloudphxx real_t operator()() { return dist_u01(engine); } }; + struct fnctr_normal01 + { + engine_t &engine; + dist_normal01_t &dist_normal01; + real_t operator()() { return dist_normal01(engine); } + }; + struct fnctr_un { engine_t &engine; @@ -45,7 +54,7 @@ namespace libcloudphxx public: // ctor - rng(int seed) : engine(seed), dist_u01(0,1), dist_un(0, std::numeric_limits::max()) {} + rng(int seed) : engine(seed), dist_u01(0,1), dist_normal01(0,1), dist_un(0, std::numeric_limits::max()) {} void generate_n( thrust_device::vector &u01, @@ -55,6 +64,14 @@ namespace libcloudphxx std::generate_n(u01.begin(), n, fnctr_u01({.engine = engine, .dist_u01 = dist_u01})); } + void generate_normal_n( + thrust_device::vector &normal01, + const thrust_size_t n + ) { + // note: generate_n copies the third argument!!! + std::generate_n(normal01.begin(), n, fnctr_normal01({.engine = engine, .dist_normal01 = dist_normal01})); + } + void generate_n( thrust_device::vector &un, const thrust_size_t n @@ -117,6 +134,26 @@ namespace libcloudphxx _unused(status); } + void generate_normal_n( + thrust_device::vector &v, + const thrust_size_t n + ) + { + int status = curandGenerateNormal(gen, thrust::raw_pointer_cast(v.data()), n, float(0), float(1)); + assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/); + _unused(status); + } + + void generate_normal_n( + thrust_device::vector &v, + const thrust_size_t n + ) + { + int status = curandGenerateNormalDouble(gen, thrust::raw_pointer_cast(v.data()), n, double(0), double(1)); + assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/); + _unused(status); + } + void generate_n( thrust_device::vector &v, const thrust_size_t n diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 65684f064..3cc513736 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -14,6 +14,8 @@ #include #include +#include + #include namespace libcloudphxx @@ -77,11 +79,19 @@ namespace libcloudphxx x, // x spatial coordinate (for 1D, 2D and 3D) y, // y spatial coordinate (for 3D) z, // z spatial coordinate (for 2D and 3D) + up, // turbulent perturbation of velocity + vp, // turbulent perturbation of velocity + wp, // turbulent perturbation of velocity + ssp, // turbulent perturbation of supersaturation + dot_ssp, // time derivative of the turbulent perturbation of supersaturation sstp_tmp_rv, // either rv_old or advection-caused change in water vapour mixing ratio sstp_tmp_th, // ditto for theta sstp_tmp_rh, // ditto for rho sstp_tmp_p; // ditto for pressure + const int no_of_n_vctrs_copied = 1; + const int no_of_real_vctrs_copied = 15; + // dry radii distribution characteristics real_t log_rd_min, // logarithm of the lower bound of the distr log_rd_max, // logarithm of the upper bound of the distr @@ -137,7 +147,10 @@ namespace libcloudphxx T, // temperature [K] p, // pressure [Pa] RH, // relative humisity - eta;// dynamic viscosity + eta,// dynamic viscosity + diss_rate; // turbulent kinetic energy dissipation rate + + real_t L; // extent of a cell, needed in tubulence, TODO: should be a ncell array, cause some cells are smaller // sorting needed only for diagnostics and coalescence bool sorted; @@ -277,6 +290,11 @@ namespace libcloudphxx n_dims == 2 ? halo_size * (opts_init.nz + 1): // 2D halo_size * (opts_init.nz + 1) * opts_init.ny // 3D ), + L( + n_dims == 1 ? common::turbulence::length_scale(opts_init.dx * si::metres) / si::metres: // 1D + n_dims == 2 ? common::turbulence::length_scale(opts_init.dx * si::metres, opts_init.dz * si::metres) / si::metres: // 2D + common::turbulence::length_scale(opts_init.dx * si::metres, opts_init.dy * si::metres, opts_init.dz * si::metres)/ si::metres // 3D + ), adve_scheme(opts_init.adve_scheme), pure_const_multi (((opts_init.sd_conc) == 0) && (opts_init.sd_const_multi > 0 || opts_init.dry_sizes.size() > 0)) // coal prob can be greater than one only in sd_conc simulations { @@ -398,6 +416,9 @@ namespace libcloudphxx void hskpng_vterm_all(); void hskpng_vterm_invalid(); + void hskpng_tke(); + void hskpng_turb_vel(const bool only_vertical = false); + void hskpng_turb_dot_ss(); void hskpng_remove_n0(); void hskpng_resize_npart(); @@ -435,18 +456,21 @@ namespace libcloudphxx ); void adve(); + void turb_adve(); template void adve_calc(bool, thrust_size_t = 0); void sedi(); void cond_dm3_helper(); - void cond(const real_t &dt, const real_t &RH_max); - void cond_sstp(const real_t &dt, const real_t &RH_max); + void cond(const real_t &dt, const real_t &RH_max, const bool turb_cond); + void cond_sstp(const real_t &dt, const real_t &RH_max, const bool turb_cond); + template + void cond_sstp_hlpr(const real_t &dt, const real_t &RH_max, const thrust_device::vector &Tp, const pres_iter &pi, const RH_iter &rhi); void update_th_rv(thrust_device::vector &); void update_state(thrust_device::vector &, thrust_device::vector &); void update_pstate(thrust_device::vector &, thrust_device::vector &); - void coal(const real_t &dt); + void coal(const real_t &dt, const bool &turb_coal); void chem_vol_ante(); void chem_flag_ante(); @@ -462,6 +486,7 @@ namespace libcloudphxx void sstp_step(const int &step); void sstp_step_exact(const int &step); + void sstp_step_ssp(const real_t &dt); void sstp_save(); void sstp_step_chem(const int &step); void sstp_save_chem(); diff --git a/src/impl/particles_impl_bcnd.ipp b/src/impl/particles_impl_bcnd.ipp index 62990cad1..eca8dd9c8 100644 --- a/src/impl/particles_impl_bcnd.ipp +++ b/src/impl/particles_impl_bcnd.ipp @@ -113,8 +113,20 @@ namespace libcloudphxx arg::_1 >= opts_init.x1 ) - rgt_id.begin(); - if(lft_count > in_n_bfr.size() || rgt_count > in_n_bfr.size()) - throw std::runtime_error(detail::formatter() << "Overflow of the in int buffer, bfr size: " << in_n_bfr.size() << " to be copied left: " << lft_count << " right: " << rgt_count); // TODO: resize buffers? + if(lft_count*no_of_n_vctrs_copied > in_n_bfr.size() || rgt_count*no_of_n_vctrs_copied > in_n_bfr.size()) + { + n_t new_size = lft_count > rgt_count ? + 1.1 * lft_count : + 1.1 * rgt_count; + + std::cerr << "Overflow of the buffer, bfr size: " << in_n_bfr.size() << " to be copied left: " << lft_count << " right: " << rgt_count << "; resizing to: " << new_size << std::endl; + + in_n_bfr.resize(no_of_n_vctrs_copied * new_size); + out_n_bfr.resize(no_of_n_vctrs_copied * new_size); + + in_real_bfr.resize(no_of_real_vctrs_copied * new_size); + out_real_bfr.resize(no_of_real_vctrs_copied * new_size); + } } // hardcoded periodic boundary in y! (TODO - as an option) diff --git a/src/impl/particles_impl_coal.ipp b/src/impl/particles_impl_coal.ipp index 84a800f5a..6e39ddddf 100644 --- a/src/impl/particles_impl_coal.ipp +++ b/src/impl/particles_impl_coal.ipp @@ -135,9 +135,10 @@ namespace libcloudphxx // read-only parameters passed to the calc function typedef thrust::tuple< real_t, // rhod (dry air density) - real_t // eta (dynamic viscosity) + real_t, // eta (dynamic viscosity) + real_t // tke dissipation rate > tpl_ro_calc_t; - enum { rhod_ix, eta_ix }; + enum { rhod_ix, eta_ix, diss_rate_ix }; const real_t dt; const kernel_base *p_kernel; @@ -236,7 +237,7 @@ namespace libcloudphxx }; template - void particles_t::impl::coal(const real_t &dt) + void particles_t::impl::coal(const real_t &dt, const bool &turb_coal) { // prerequisites hskpng_shuffle_and_sort(); // to get random neighbours by default @@ -331,13 +332,6 @@ namespace libcloudphxx > > zip_rw_t; - typedef thrust::zip_iterator< - thrust::tuple< - pi_real_t, // rhod - pi_real_t // eta - > - > zip_ro_calc_t; //read-only parameters passed to the calc() function, later also epsilon and Re_lambda - zip_ro_t zip_ro_it( thrust::make_tuple( // u01 @@ -355,14 +349,31 @@ namespace libcloudphxx ) ); - zip_ro_calc_t zip_ro_calc_it( - thrust::make_tuple( - // rhod - thrust::make_permutation_iterator(rhod.begin(), sorted_ijk.begin()), - // eta - thrust::make_permutation_iterator(eta.begin(), sorted_ijk.begin()) + auto zip_ro_calc_it = + thrust::make_zip_iterator( + thrust::make_tuple( + // rhod + thrust::make_permutation_iterator(rhod.begin(), sorted_ijk.begin()), + // eta + thrust::make_permutation_iterator(eta.begin(), sorted_ijk.begin()), + // tke dissipation rate + thrust::make_permutation_iterator(thrust::make_constant_iterator(0), sorted_ijk.begin()) + ) ) - ); + ; + + auto zip_ro_calc_turb_it = + thrust::make_zip_iterator( + thrust::make_tuple( + // rhod + thrust::make_permutation_iterator(rhod.begin(), sorted_ijk.begin()), + // eta + thrust::make_permutation_iterator(eta.begin(), sorted_ijk.begin()), + // tke dissipation rate + thrust::make_permutation_iterator(diss_rate.begin(), sorted_ijk.begin()) + ) + ) + ; zip_rw_t zip_rw_it( thrust::make_tuple( @@ -384,11 +395,20 @@ namespace libcloudphxx ) ); - thrust::for_each( - thrust::make_zip_iterator(thrust::make_tuple(zip_ro_it, zip_rw_it, zip_ro_calc_it)), - thrust::make_zip_iterator(thrust::make_tuple(zip_ro_it, zip_rw_it, zip_ro_calc_it)) + n_part - 1, - detail::collider(dt, p_kernel, pure_const_multi, increase_sstp_coal) - ); + + if(turb_coal) + thrust::for_each( + thrust::make_zip_iterator(thrust::make_tuple(zip_ro_it, zip_rw_it, zip_ro_calc_turb_it)), + thrust::make_zip_iterator(thrust::make_tuple(zip_ro_it, zip_rw_it, zip_ro_calc_turb_it)) + n_part - 1, + detail::collider(dt, p_kernel, pure_const_multi, increase_sstp_coal) + ); + else + thrust::for_each( + thrust::make_zip_iterator(thrust::make_tuple(zip_ro_it, zip_rw_it, zip_ro_calc_it)), + thrust::make_zip_iterator(thrust::make_tuple(zip_ro_it, zip_rw_it, zip_ro_calc_it)) + n_part - 1, + detail::collider(dt, p_kernel, pure_const_multi, increase_sstp_coal) + ); + // nancheck(n, "n - post coalescence"); nancheck(rw2, "rw2 - post coalescence"); nancheck(rd3, "rd3 - post coalescence"); diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index 2a83e6569..99b322e74 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -12,9 +12,12 @@ namespace libcloudphxx template void particles_t::impl::cond( const real_t &dt, - const real_t &RH_max + const real_t &RH_max, + const bool turb_cond ) { + namespace arg = thrust::placeholders; + // --- calc liquid water content before cond --- hskpng_sort(); thrust_device::vector &drv(tmp_device_real_cell); @@ -33,25 +36,54 @@ namespace libcloudphxx thrust::negate() ); + auto hlpr_zip_iter = thrust::make_zip_iterator(thrust::make_tuple( + thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), + thrust::make_permutation_iterator(rv.begin(), ijk.begin()), + thrust::make_permutation_iterator(T.begin(), ijk.begin()), + thrust::make_permutation_iterator(eta.begin(), ijk.begin()), + rd3.begin(), + kpa.begin(), + vt.begin() + )); + // calculating drop growth in a timestep using backward Euler - thrust::transform( - rw2.begin(), rw2.end(), // input - 1st arg (zip not as 1st arg not to write zip.end() - thrust::make_zip_iterator( // input - 2nd arg - thrust::make_tuple( - thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), - thrust::make_permutation_iterator(rv.begin(), ijk.begin()), - thrust::make_permutation_iterator(T.begin(), ijk.begin()), - thrust::make_permutation_iterator(p.begin(), ijk.begin()), - thrust::make_permutation_iterator(RH.begin(), ijk.begin()), - thrust::make_permutation_iterator(eta.begin(), ijk.begin()), - rd3.begin(), - kpa.begin(), - vt.begin() - ) - ), - rw2.begin(), // output - detail::advance_rw2(dt, RH_max) - ); + // TODO: both calls almost identical, use std::bind or sth? + if(turb_cond) + { + thrust_device::vector &RH_plus_ssp(tmp_device_real_part2); + thrust::transform( + ssp.begin(), ssp.end(), + thrust::make_permutation_iterator(RH.begin(), ijk.begin()), + RH_plus_ssp.begin(), + arg::_1 + arg::_2 + ); + + thrust::transform( + rw2.begin(), rw2.end(), // input - 1st arg (zip not as 1st arg not to write zip.end() + thrust::make_zip_iterator( // input - 2nd arg + thrust::make_tuple( + hlpr_zip_iter, + thrust::make_permutation_iterator(p.begin(), ijk.begin()), + RH_plus_ssp.begin() + ) + ), + rw2.begin(), // output + detail::advance_rw2(dt, RH_max) + ); + } + else + thrust::transform( + rw2.begin(), rw2.end(), // input - 1st arg (zip not as 1st arg not to write zip.end() + thrust::make_zip_iterator( // input - 2nd arg + thrust::make_tuple( + hlpr_zip_iter, + thrust::make_permutation_iterator(p.begin(), ijk.begin()), + thrust::make_permutation_iterator(RH.begin(), ijk.begin()) + ) + ), + rw2.begin(), // output + detail::advance_rw2(dt, RH_max) + ); nancheck(rw2, "rw2 after condensation (no sub-steps"); // calculating the 3rd wet moment after condensation diff --git a/src/impl/particles_impl_cond_common.ipp b/src/impl/particles_impl_cond_common.ipp index e795a9e09..01ea1f1b3 100644 --- a/src/impl/particles_impl_cond_common.ipp +++ b/src/impl/particles_impl_cond_common.ipp @@ -69,20 +69,20 @@ namespace libcloudphxx advance_rw2_minfun( const real_t &dt, const real_t &rw2, - const thrust::tuple &tpl, + const thrust::tuple, real_t, real_t> &tpl, const real_t &RH_max ) : dt(dt * si::seconds), rw2_old(rw2 * si::square_metres), - rhod( thrust::get<0>(tpl) * si::kilograms / si::cubic_metres), - rv( thrust::get<1>(tpl)), - T( thrust::get<2>(tpl) * si::kelvins), - p( thrust::get<3>(tpl) * si::pascals), - RH( thrust::get<4>(tpl)), - eta( thrust::get<5>(tpl) * si::pascals * si::seconds), - rd3( thrust::get<6>(tpl) * si::cubic_metres), - kpa( thrust::get<7>(tpl)), - vt( thrust::get<8>(tpl) * si::metres_per_second), + rhod( thrust::get<0>(thrust::get<0>(tpl)) * si::kilograms / si::cubic_metres), + rv( thrust::get<1>(thrust::get<0>(tpl))), + T( thrust::get<2>(thrust::get<0>(tpl)) * si::kelvins), + eta( thrust::get<3>(thrust::get<0>(tpl)) * si::pascals * si::seconds), + rd3( thrust::get<4>(thrust::get<0>(tpl)) * si::cubic_metres), + kpa( thrust::get<5>(thrust::get<0>(tpl))), + vt( thrust::get<6>(thrust::get<0>(tpl)) * si::metres_per_second), + p( thrust::get<1>(tpl) * si::pascals), + RH( thrust::get<2>(tpl)), RH_max(RH_max) {} @@ -100,7 +100,9 @@ namespace libcloudphxx using common::mean_free_path::lambda_K; using common::ventil::Sh; using common::ventil::Nu; +#if !defined(__NVCC__) using std::sqrt; +#endif const quantity rw = sqrt(real_t(rw2 / si::square_metres)) * si::metres; const quantity rw3 = rw * rw * rw;; @@ -152,7 +154,7 @@ namespace libcloudphxx BOOST_GPU_ENABLED real_t operator()( const real_t &rw2_old, - const thrust::tuple &tpl + const thrust::tuple, real_t, real_t> &tpl ) const { #if !defined(__NVCC__) using std::min; @@ -161,6 +163,7 @@ namespace libcloudphxx using std::abs; #endif + auto& tpl_in = thrust::get<0>(tpl); const advance_rw2_minfun f(dt, rw2_old, tpl, RH_max); const real_t drw2 = dt * f.drw2_dt(rw2_old * si::square_metres) * si::seconds / si::square_metres; @@ -171,22 +174,22 @@ namespace libcloudphxx printf("rw2_old: %g\n",rw2_old); printf("dt: %g\n",dt); printf("RH_max: %g\n",RH_max); - printf("rhod: %g\n",thrust::get<0>(tpl)); - printf("rv: %g\n",thrust::get<1>(tpl)); - printf("T: %g\n",thrust::get<2>(tpl)); - printf("p: %g\n",thrust::get<3>(tpl)); - printf("RH: %g\n",thrust::get<4>(tpl)); - printf("eta: %g\n",thrust::get<5>(tpl)); - printf("rd3: %g\n",thrust::get<6>(tpl)); - printf("kpa: %g\n",thrust::get<7>(tpl)); - printf("vt: %g\n",thrust::get<8>(tpl)); + printf("rhod: %g\n",thrust::get<0>(tpl_in)); + printf("rv: %g\n",thrust::get<1>(tpl_in)); + printf("T: %g\n",thrust::get<2>(tpl_in)); + printf("eta: %g\n",thrust::get<3>(tpl_in)); + printf("rd3: %g\n",thrust::get<4>(tpl_in)); + printf("kpa: %g\n",thrust::get<5>(tpl_in)); + printf("vt: %g\n",thrust::get<6>(tpl_in)); + printf("p: %g\n",thrust::get<1>(tpl)); + printf("RH: %g\n",thrust::get<2>(tpl)); assert(0); } #endif if (drw2 == 0) return rw2_old; - const real_t rd2 = pow(thrust::get<6>(tpl), real_t(2./3)); + const real_t rd2 = pow(thrust::get<4>(tpl_in), real_t(2./3)); const real_t a = max(rd2, rw2_old + min(real_t(0), config.cond_mlt * drw2)), diff --git a/src/impl/particles_impl_cond_sstp.ipp b/src/impl/particles_impl_cond_sstp.ipp index 7227b70be..4d9a3a53f 100644 --- a/src/impl/particles_impl_cond_sstp.ipp +++ b/src/impl/particles_impl_cond_sstp.ipp @@ -9,11 +9,72 @@ namespace libcloudphxx { namespace lgrngn { + namespace detail + { + template + struct RH_sgs : thrust::unary_function&, real_t> + { + RH resolved_RH; + + RH_sgs(RH_formula_t::RH_formula_t RH_formula): + resolved_RH(RH_formula) + {} + + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple &tpl) + { + return resolved_RH(thrust::make_tuple(thrust::get<0>(tpl), thrust::get<1>(tpl), thrust::get<2>(tpl))) + thrust::get<3>(tpl); + } + }; + }; + + template + template + void particles_t::impl::cond_sstp_hlpr( + const real_t &dt, + const real_t &RH_max, + const thrust_device::vector &Tp, + const pres_iter &pi, + const RH_iter &rhi + ) { + auto hlpr_zip_iter = thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_rh.begin(), + sstp_tmp_rv.begin(), + Tp.begin(), + // particle-specific eta + thrust::make_transform_iterator( + Tp.begin(), + detail::common__vterm__visc() + ), + rd3.begin(), + kpa.begin(), + vt.begin() + )); + + // calculating drop growth in a timestep using backward Euler + thrust::transform( + rw2.begin(), rw2.end(), // input - 1st arg (zip not as 1st arg not to write zip.end() + thrust::make_zip_iterator(thrust::make_tuple( // input - 2nd arg + hlpr_zip_iter, + // particle-specific p + pi, + // particle-specific RH + rhi + )), + rw2.begin(), // output + detail::advance_rw2(dt, RH_max) + ); + } + + template void particles_t::impl::cond_sstp( const real_t &dt, - const real_t &RH_max - ) { + const real_t &RH_max, + const bool turb_cond + ) { + namespace arg = thrust::placeholders; + // prerequisite hskpng_sort(); // particle's local change in rv @@ -55,91 +116,80 @@ namespace libcloudphxx // calculating drop growth in a timestep using backward Euler - // TODO: these two function calls only differ in the way pressure is calculated, - // find a way to reuse it (e.g. std::bind?) - // TODO2: in const_p we don't substep pressure if(!const_p) { // particle-specific pressure iterator, used twice // TODO: store the value somewhere? - auto pressure_iter = - thrust::make_transform_iterator( - thrust::make_zip_iterator( - thrust::make_tuple( - sstp_tmp_rh.begin(), - sstp_tmp_rv.begin(), - Tp.begin() - )), - detail::common__theta_dry__p() - ); + auto pressure_iter = thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_rh.begin(), + sstp_tmp_rv.begin(), + Tp.begin() + )), + detail::common__theta_dry__p() + ); - thrust::transform( - rw2.begin(), rw2.end(), // input - 1st arg (zip not as 1st arg not to write zip.end() - thrust::make_zip_iterator( // input - 2nd arg - thrust::make_tuple( - sstp_tmp_rh.begin(), - sstp_tmp_rv.begin(), - Tp.begin(), - // particle-specific p - pressure_iter, - // particle-specific RH - thrust::make_transform_iterator( - thrust::make_zip_iterator( - thrust::make_tuple( - pressure_iter, - sstp_tmp_rv.begin(), - Tp.begin() - )), - detail::RH(opts_init.RH_formula) - ), - // particle-specific eta - thrust::make_transform_iterator( + if(turb_cond) + cond_sstp_hlpr(dt, RH_max, Tp, + // particle-specific p + pressure_iter, + // particle-specific RH, resolved + SGS + thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple( + pressure_iter, + sstp_tmp_rv.begin(), Tp.begin(), - detail::common__vterm__visc() - ), - rd3.begin(), - kpa.begin(), - vt.begin() - ) - ), - rw2.begin(), // output - detail::advance_rw2(dt, RH_max) - ); + ssp.begin() + )), + detail::RH_sgs(opts_init.RH_formula) + ) + ); + else // no RH SGS + cond_sstp_hlpr(dt, RH_max, Tp, + // particle-specific p + pressure_iter, + // particle-specific RH, resolved + thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple( + pressure_iter, + sstp_tmp_rv.begin(), + Tp.begin() + )), + detail::RH(opts_init.RH_formula) + ) + ); } - else + else // const_p { - thrust::transform( - rw2.begin(), rw2.end(), // input - 1st arg (zip not as 1st arg not to write zip.end() - thrust::make_zip_iterator( // input - 2nd arg - thrust::make_tuple( - sstp_tmp_rh.begin(), - sstp_tmp_rv.begin(), - Tp.begin(), - // particle-specific p - sstp_tmp_p.begin(), - // particle-specific RH - thrust::make_transform_iterator( - thrust::make_zip_iterator( - thrust::make_tuple( - sstp_tmp_p.begin(), - sstp_tmp_rv.begin(), - Tp.begin() - )), - detail::RH(opts_init.RH_formula) - ), - // particle-specific eta - thrust::make_transform_iterator( + if(turb_cond) + cond_sstp_hlpr(dt, RH_max, Tp, + // particle-specific p + sstp_tmp_p.begin(), + // particle-specific RH, resolved + SGS + thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_p.begin(), + sstp_tmp_rv.begin(), Tp.begin(), - detail::common__vterm__visc() - ), - rd3.begin(), - kpa.begin(), - vt.begin() - ) - ), - rw2.begin(), // output - detail::advance_rw2(dt, RH_max) - ); + ssp.begin() + )), + detail::RH_sgs(opts_init.RH_formula) + ) + ); + else // no RH SGS + cond_sstp_hlpr(dt, RH_max, Tp, + // particle-specific p + sstp_tmp_p.begin(), + // particle-specific RH, resolved + thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_p.begin(), + sstp_tmp_rv.begin(), + Tp.begin() + )), + detail::RH(opts_init.RH_formula) + ) + ); } // calc rw3_new - rw3_old diff --git a/src/impl/particles_impl_dist_analysis.ipp b/src/impl/particles_impl_dist_analysis.ipp index 23797f6a8..74cafaafd 100644 --- a/src/impl/particles_impl_dist_analysis.ipp +++ b/src/impl/particles_impl_dist_analysis.ipp @@ -55,6 +55,10 @@ namespace libcloudphxx else if (n_max == 0) rd_max /= 1.01; else found_optimal_range = true; } +#if !defined(__NVCC__) + using std::max; +#endif + log_rd_min = max(log_rd_min, real_t(log(opts_init.rd_min))); // user-defined lower limit for rd }; template @@ -89,6 +93,10 @@ namespace libcloudphxx detail::eval_and_add(n_of_lnrd_stp, -init_dist_bound_value)(real_t(log(config.rd_max_init))), config.eps_tolerance, n_iter ); +#if !defined(__NVCC__) + using std::max; +#endif + log_rd_min = max(log_rd_min, real_t(log(opts_init.rd_min))); // user-defined lower limit for rd } }; }; diff --git a/src/impl/particles_impl_hskpng_remove.ipp b/src/impl/particles_impl_hskpng_remove.ipp index c6e2a7b79..522f137f3 100644 --- a/src/impl/particles_impl_hskpng_remove.ipp +++ b/src/impl/particles_impl_hskpng_remove.ipp @@ -6,6 +6,7 @@ */ #include +#include namespace libcloudphxx { @@ -20,25 +21,39 @@ namespace libcloudphxx { // TODO: init these vctrs once per run thrust_device::vector * real_t_vctrs_a[] = {&rd3, &rw2, &kpa, &vt}; - std::vector*> real_t_vctrs(&real_t_vctrs_a[0], &real_t_vctrs_a[0]+4); - std::vector*> n_t_vctrs; - n_t_vctrs.push_back(&ijk); + std::set*> real_t_vctrs(&real_t_vctrs_a[0], &real_t_vctrs_a[0]+4); + std::set*> n_t_vctrs; + n_t_vctrs.insert(&ijk); - if (opts_init.nx != 0) n_t_vctrs.push_back(&i); - if (opts_init.ny != 0) n_t_vctrs.push_back(&j); - if (opts_init.nz != 0) n_t_vctrs.push_back(&k); + if (opts_init.nx != 0) n_t_vctrs.insert(&i); + if (opts_init.ny != 0) n_t_vctrs.insert(&j); + if (opts_init.nz != 0) n_t_vctrs.insert(&k); - if (opts_init.nx != 0) real_t_vctrs.push_back(&x); - if (opts_init.ny != 0) real_t_vctrs.push_back(&y); - if (opts_init.nz != 0) real_t_vctrs.push_back(&z); + if (opts_init.nx != 0) real_t_vctrs.insert(&x); + if (opts_init.ny != 0) real_t_vctrs.insert(&y); + if (opts_init.nz != 0) real_t_vctrs.insert(&z); + + if(opts_init.turb_adve_switch) + { + if (opts_init.nx != 0) real_t_vctrs.insert(&up); + if (opts_init.ny != 0) real_t_vctrs.insert(&vp); + if (opts_init.nz != 0) real_t_vctrs.insert(&wp); + } + + if(opts_init.turb_cond_switch) + { + real_t_vctrs.insert(&wp); + real_t_vctrs.insert(&ssp); + real_t_vctrs.insert(&dot_ssp); + } if(opts_init.sstp_cond>1 && opts_init.exact_sstp_cond) { - real_t_vctrs.push_back(&sstp_tmp_th); - real_t_vctrs.push_back(&sstp_tmp_rv); - real_t_vctrs.push_back(&sstp_tmp_rh); + real_t_vctrs.insert(&sstp_tmp_th); + real_t_vctrs.insert(&sstp_tmp_rv); + real_t_vctrs.insert(&sstp_tmp_rh); if(const_p) - real_t_vctrs.push_back(&sstp_tmp_p); + real_t_vctrs.insert(&sstp_tmp_p); } namespace arg = thrust::placeholders; diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index d12e85eed..0a74ac83d 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -33,11 +33,25 @@ namespace libcloudphxx if (opts_init.ny != 0) y.resize(n_part); if (opts_init.nz != 0) z.resize(n_part); + if(opts_init.turb_adve_switch) + { + if (opts_init.nx != 0) up.resize(n_part); + if (opts_init.ny != 0) vp.resize(n_part); + if (opts_init.nz != 0) wp.resize(n_part); + } + + if(opts_init.turb_cond_switch) + { + wp.resize(n_part); + ssp.resize(n_part); + dot_ssp.resize(n_part); + } + if(opts_init.chem_switch || opts_init.sstp_cond > 1 || n_dims >= 2) { tmp_device_real_part1.resize(n_part); } - if((opts_init.sstp_cond>1 && opts_init.exact_sstp_cond) || n_dims==3) + if((opts_init.sstp_cond>1 && opts_init.exact_sstp_cond) || n_dims==3 || opts_init.turb_cond_switch) { tmp_device_real_part2.resize(n_part); } diff --git a/src/impl/particles_impl_hskpng_tke.ipp b/src/impl/particles_impl_hskpng_tke.ipp new file mode 100644 index 000000000..c5ead9db1 --- /dev/null +++ b/src/impl/particles_impl_hskpng_tke.ipp @@ -0,0 +1,39 @@ +// vim:filetype=cpp +/** @file + * @copyright University of Warsaw + * @section LICENSE + * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) + */ + +#include + +namespace libcloudphxx +{ + namespace lgrngn + { + namespace detail + { + template + struct common__turbulence__tke + { + const quantity L; + common__turbulence__tke(const real_t &L): + L(L * si::metres){} + + BOOST_GPU_ENABLED + real_t operator()(const real_t &diss_rate) + { + return common::turbulence::tke( + diss_rate * si::metres * si::metres / si::seconds / si::seconds / si::seconds, + L) / si::metres / si::metres * si::seconds * si::seconds; + } + }; + }; + // calc the SGS TKE, in place of dissipation rate + template + void particles_t::impl::hskpng_tke() + { + thrust::transform(diss_rate.begin(), diss_rate.end(), diss_rate.begin(), detail::common__turbulence__tke(L)); + } + }; +}; diff --git a/src/impl/particles_impl_hskpng_turb_ss.ipp b/src/impl/particles_impl_hskpng_turb_ss.ipp new file mode 100644 index 000000000..429402088 --- /dev/null +++ b/src/impl/particles_impl_hskpng_turb_ss.ipp @@ -0,0 +1,94 @@ +// vim:filetype=cpp +/** @file + * @copyright University of Warsaw + * @section LICENSE + * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) + */ + +#include + +namespace libcloudphxx +{ + namespace lgrngn + { + namespace detail + { + struct common__turbulence__tau_relax + { + template + BOOST_GPU_ENABLED + real_t operator()(const real_t &wet_mom_1, const real_t &dv) + { + return common::turbulence::tau_relax(wet_mom_1 / dv / si::square_metres) / si::seconds; + } + }; + + template + struct common__turbulence__turb_dot_ss + { + template + BOOST_GPU_ENABLED + real_t operator()(tpl_t tpl) + { + assert(thrust::get<2>(tpl) > real_t(0)); + return common::turbulence::dot_turb_ss( + quantity(thrust::get<0>(tpl)), + thrust::get<1>(tpl) * si::metres / si::seconds, + thrust::get<2>(tpl) * si::seconds + ) * si::seconds; + } + }; + }; + + // calc the SGS turbulent supersaturation + template + void particles_t::impl::hskpng_turb_dot_ss() + { + thrust_device::vector &tau_rlx(tmp_device_real_cell); // tau_rlx needs to have length n_cell +#if !defined(NDEBUG) + // fill with a dummy value for debugging + thrust::fill(tau_rlx.begin(), tau_rlx.end(), -44); +#endif + + // calc relaxation times (only in cells that contain any SDs) stored in count_mom + moms_all(); + moms_calc(rw2.begin(), real_t(1./2), false); + thrust::transform( + count_mom.begin(), + count_mom.begin() + count_n, + thrust::make_permutation_iterator( + dv.begin(), + count_ijk.begin() + ), + count_mom.begin(), + detail::common__turbulence__tau_relax() + ); + + // copy the calculated relaxation times to the array of length n_cell + thrust::copy( + count_mom.begin(), + count_mom.begin() + count_n, + thrust::make_permutation_iterator( + tau_rlx.begin(), + count_ijk.begin() + ) + ); + + // calc ss perturb for each SD + thrust::transform( + thrust::make_zip_iterator(thrust::make_tuple( + ssp.begin(), + wp.begin(), + thrust::make_permutation_iterator(tau_rlx.begin(), ijk.begin()) + )), + thrust::make_zip_iterator(thrust::make_tuple( + ssp.begin(), + wp.begin(), + thrust::make_permutation_iterator(tau_rlx.begin(), ijk.begin()) + )) + n_part, + dot_ssp.begin(), + detail::common__turbulence__turb_dot_ss() + ); + } + }; +}; diff --git a/src/impl/particles_impl_hskpng_turb_vel.ipp b/src/impl/particles_impl_hskpng_turb_vel.ipp new file mode 100644 index 000000000..8ec33ec47 --- /dev/null +++ b/src/impl/particles_impl_hskpng_turb_vel.ipp @@ -0,0 +1,85 @@ +// vim:filetype=cpp +/** @file + * @copyright University of Warsaw + * @section LICENSE + * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) + */ + +#include + +namespace libcloudphxx +{ + namespace lgrngn + { + namespace detail + { + template + struct common__turbulence__tau + { + const quantity L; + common__turbulence__tau(const real_t &L): + L(L * si::metres){} + + BOOST_GPU_ENABLED + real_t operator()(const real_t &tke) + { + return common::turbulence::tau( + tke * si::metres * si::metres / si::seconds / si::seconds, + L) / si::seconds; + } + }; + + template + struct common__turbulence__update_turb_vel + { + const quantity dt; + common__turbulence__update_turb_vel(const real_t &dt): + dt(dt * si::seconds){} + + template + BOOST_GPU_ENABLED + real_t operator()(tpl_t tpl) + { + return common::turbulence::update_turb_vel( + thrust::get<0>(tpl) * si::metres / si::seconds, + thrust::get<1>(tpl) * si::seconds, + dt, + thrust::get<2>(tpl) * si::metres * si::metres / si::seconds / si::seconds, + quantity(thrust::get<3>(tpl)) + ) / si::metres * si::seconds; + } + }; + }; + // calc the SGS turbulent velocity component + template + void particles_t::impl::hskpng_turb_vel(const bool only_vertical) + { + thrust_device::vector &tau(tmp_device_real_cell); + thrust_device::vector &tke(diss_rate); // should be called after hskpng_tke, which replaces diss_rate with tke + thrust::transform(tke.begin(), tke.end(), tau.begin(), detail::common__turbulence__tau(L)); + + thrust_device::vector &r_normal(tmp_device_real_part); + thrust_device::vector * vel_turbs_vctrs_a[] = {&up, &wp, &vp}; + for(int i = (only_vertical ? 1 : 0); i < (only_vertical ? 2 : n_dims); ++i) + { + rng.generate_normal_n(r_normal, n_part); // generate a random number for wach particle with a normal distribution with mean 0 and std dev 1 + thrust::transform( + thrust::make_zip_iterator(thrust::make_tuple( + vel_turbs_vctrs_a[i]->begin(), + thrust::make_permutation_iterator(tau.begin(), ijk.begin()), + thrust::make_permutation_iterator(tke.begin(), ijk.begin()), + r_normal.begin() + )), + thrust::make_zip_iterator(thrust::make_tuple( + vel_turbs_vctrs_a[i]->begin(), + thrust::make_permutation_iterator(tau.begin(), ijk.begin()), + thrust::make_permutation_iterator(tke.begin(), ijk.begin()), + r_normal.begin() + )) + n_part, + vel_turbs_vctrs_a[i]->begin(), + detail::common__turbulence__update_turb_vel(opts_init.dt) + ); + } + } + }; +}; diff --git a/src/impl/particles_impl_init_hskpng_npart.ipp b/src/impl/particles_impl_init_hskpng_npart.ipp index d06a74d87..5b858ab49 100644 --- a/src/impl/particles_impl_init_hskpng_npart.ipp +++ b/src/impl/particles_impl_init_hskpng_npart.ipp @@ -23,6 +23,21 @@ namespace libcloudphxx if (opts_init.ny != 0) y.reserve(opts_init.n_sd_max); if (opts_init.nz != 0) z.reserve(opts_init.n_sd_max); + // using resize instead of reserve is ok, because hskpng_resize is called right after this init + if(opts_init.turb_adve_switch) + { + if (opts_init.nx != 0) up.resize(opts_init.n_sd_max, 0.); // init with no perturbation + if (opts_init.ny != 0) vp.resize(opts_init.n_sd_max, 0.); + if (opts_init.nz != 0) wp.resize(opts_init.n_sd_max, 0.); + } + + if(opts_init.turb_cond_switch) + { + wp.resize(opts_init.n_sd_max, 0.); + ssp.resize(opts_init.n_sd_max, 0.); + dot_ssp.resize(opts_init.n_sd_max, 0.); + } + vt.resize(opts_init.n_sd_max, 0.); // so that it may be safely used in condensation before first update sorted_id.reserve(opts_init.n_sd_max); @@ -41,7 +56,7 @@ namespace libcloudphxx { tmp_device_real_part1.reserve(opts_init.n_sd_max); } - if((opts_init.sstp_cond>1 && opts_init.exact_sstp_cond) || n_dims==3) + if((opts_init.sstp_cond>1 && opts_init.exact_sstp_cond) || n_dims==3 || opts_init.turb_cond_switch) { tmp_device_real_part2.reserve(opts_init.n_sd_max); } @@ -61,11 +76,11 @@ namespace libcloudphxx // reserve memory for in/out buffers if(opts_init.dev_count > 1) { - in_n_bfr.resize(opts_init.n_sd_max / opts_init.nx / config.bfr_fraction); // for n - out_n_bfr.resize(opts_init.n_sd_max / opts_init.nx / config.bfr_fraction); + in_n_bfr.resize(no_of_n_vctrs_copied * opts_init.n_sd_max / opts_init.nx / config.bfr_fraction); // for n + out_n_bfr.resize(no_of_n_vctrs_copied * opts_init.n_sd_max / opts_init.nx / config.bfr_fraction); - in_real_bfr.resize(11 * opts_init.n_sd_max / opts_init.nx / config.bfr_fraction); // for rd3 rw2 kpa vt x y z sstp_tmp_th/rv/rh/p - out_real_bfr.resize(11 * opts_init.n_sd_max / opts_init.nx / config.bfr_fraction); + in_real_bfr.resize(no_of_real_vctrs_copied * opts_init.n_sd_max / opts_init.nx / config.bfr_fraction); // for rd3 rw2 kpa vt x y z sstp_tmp_th/rv/rh/p + out_real_bfr.resize(no_of_real_vctrs_copied * opts_init.n_sd_max / opts_init.nx / config.bfr_fraction); } } }; diff --git a/src/impl/particles_impl_init_kernel.ipp b/src/impl/particles_impl_init_kernel.ipp index 7e7efdfb2..9b4ab8eee 100644 --- a/src/impl/particles_impl_init_kernel.ipp +++ b/src/impl/particles_impl_init_kernel.ipp @@ -182,10 +182,12 @@ namespace libcloudphxx //Onishi turbulent kernel (Onishi 2015 JAS) with Hall, Davis and Jones (no van der Waals) efficiencies case(kernel_t::onishi_hall_davis_no_waals): - if(n_user_params != 2) + if(n_user_params != 1) { - throw std::runtime_error("Please supply two kernel parameters: rate of dissipation epsilon [m^2/s^3] and Taylor microscale Reynolds number."); + throw std::runtime_error("Please supply one kernel parameter: Taylor microscale Reynolds number."); } + if(!opts_init.turb_coal_switch) + throw std::runtime_error("To use the turbulent Onishis kernel, set turb_coal_switch=True"); //read in kernel efficiencies to a temporary container detail::hall_davis_no_waals_efficiencies (tmp_kernel_eff); @@ -205,10 +207,12 @@ namespace libcloudphxx //Onishi turbulent kernel (Onishi 2015 JAS) with Hall efficiencies case(kernel_t::onishi_hall): - if(n_user_params != 2) + if(n_user_params != 1) { - throw std::runtime_error("Please supply two kernel parameters: rate of dissipation epsilon [m^2/s^3] and Taylor microscale Reynolds number."); + throw std::runtime_error("Please supply one kernel parameter: Taylor microscale Reynolds number."); } + if(!opts_init.turb_coal_switch) + throw std::runtime_error("To use the turbulent Onishis kernel, set turb_coal_switch=True"); //read in kernel efficiencies to a temporary container detail::hall_efficiencies (tmp_kernel_eff); diff --git a/src/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index ca3e0cc8d..616792eab 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/src/impl/particles_impl_init_sync.ipp @@ -17,11 +17,14 @@ namespace libcloudphxx p.resize(n_cell); th.resize(n_cell); rv.resize(n_cell); - for (int i = 0; i < chem_gas_n; ++i) - ambient_chem[(chem_species_t)i].resize(n_cell); + if(opts_init.chem_switch) + for (int i = 0; i < chem_gas_n; ++i) + ambient_chem[(chem_species_t)i].resize(n_cell); + if(opts_init.turb_cond_switch || opts_init.turb_adve_switch || opts_init.turb_coal_switch) + diss_rate.resize(n_cell); // memory allocation for vector fields (Arakawa-C grid) - // 1-cell halo in the x dimensions (which could be distmem boundary) + // halo in the x dimensions (which could be distmem boundary) switch (n_dims) { case 3: diff --git a/src/impl/particles_impl_rcyc.ipp b/src/impl/particles_impl_rcyc.ipp index f77a9732a..b68704d35 100644 --- a/src/impl/particles_impl_rcyc.ipp +++ b/src/impl/particles_impl_rcyc.ipp @@ -118,6 +118,20 @@ namespace libcloudphxx for (int i = 0; i < chem_all; ++i) detail::copy_prop(chem_bgn[i], sorted_id, n_flagged); } + + if(opts_init.turb_adve_switch) + { + if (opts_init.nx > 0) detail::copy_prop(up.begin(), sorted_id, n_flagged); + if (opts_init.ny > 0) detail::copy_prop(vp.begin(), sorted_id, n_flagged); + if (opts_init.nz > 0) detail::copy_prop(wp.begin(), sorted_id, n_flagged); + } + + if(opts_init.turb_cond_switch) + { + if(!(opts_init.turb_adve_switch && opts_init.nz > 0)) detail::copy_prop(wp.begin(), sorted_id, n_flagged); + detail::copy_prop(ssp.begin(), sorted_id, n_flagged); + detail::copy_prop(dot_ssp.begin(), sorted_id, n_flagged); + } { namespace arg = thrust::placeholders; diff --git a/src/impl/particles_impl_sstp.ipp b/src/impl/particles_impl_sstp.ipp index e69d66044..26464c5ae 100644 --- a/src/impl/particles_impl_sstp.ipp +++ b/src/impl/particles_impl_sstp.ipp @@ -140,5 +140,21 @@ namespace libcloudphxx } } } + + + template + void particles_t::impl::sstp_step_ssp( + const real_t &dt + ) + { + namespace arg = thrust::placeholders; + + thrust::transform( + ssp.begin(), ssp.end(), + dot_ssp.begin(), + ssp.begin(), + arg::_1 + arg::_2 * dt + ); + } }; }; diff --git a/src/impl/particles_impl_turb_adve.ipp b/src/impl/particles_impl_turb_adve.ipp new file mode 100644 index 000000000..8782c129f --- /dev/null +++ b/src/impl/particles_impl_turb_adve.ipp @@ -0,0 +1,36 @@ +// 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 + { + // calc the SGS turbulent velocity component + template + void particles_t::impl::turb_adve() + { + thrust_device::vector * vel_pos_a[] = {&x, &y, &z}; + std::vector*> vel_pos(&vel_pos_a[0], &vel_pos_a[0]+n_dims); + + thrust_device::vector * vel_turbs_vctrs_a[] = {&up, &wp, &vp}; + std::vector*> vel_turbs_vctrs(&vel_turbs_vctrs_a[0], &vel_turbs_vctrs_a[0]+n_dims); + + namespace arg = thrust::placeholders; + + for(int i=0; ibegin(), + vel_pos[i]->end(), + vel_turbs_vctrs[i]->begin(), + vel_pos[i]->begin(), + arg::_1 + arg::_2 * opts_init.dt + ); + } + } + }; +}; diff --git a/src/impl_multi_gpu/particles_multi_gpu_impl_step_async_and_copy.ipp b/src/impl_multi_gpu/particles_multi_gpu_impl_step_async_and_copy.ipp index cab421ead..593fc0cee 100644 --- a/src/impl_multi_gpu/particles_multi_gpu_impl_step_async_and_copy.ipp +++ b/src/impl_multi_gpu/particles_multi_gpu_impl_step_async_and_copy.ipp @@ -66,6 +66,10 @@ namespace libcloudphxx thrust_device::vector &x(particles[dev_id]->pimpl->x); thrust_device::vector &y(particles[dev_id]->pimpl->y); thrust_device::vector &z(particles[dev_id]->pimpl->z); + thrust_device::vector &up(particles[dev_id]->pimpl->up); + thrust_device::vector &vp(particles[dev_id]->pimpl->vp); + thrust_device::vector &wp(particles[dev_id]->pimpl->wp); + thrust_device::vector &ssp(particles[dev_id]->pimpl->ssp); thrust_device::vector &rd3(particles[dev_id]->pimpl->rd3); thrust_device::vector &rw2(particles[dev_id]->pimpl->rw2); thrust_device::vector &kpa(particles[dev_id]->pimpl->kpa); @@ -135,6 +139,18 @@ namespace libcloudphxx if(particles[dev_id]->pimpl->const_p) real_t_vctrs.push_back(&sstp_tmp_p); } + if(glob_opts_init.turb_adve_switch) + { + if(glob_opts_init.nx > 0) real_t_vctrs.push_back(&up); + if(glob_opts_init.ny > 0) real_t_vctrs.push_back(&vp); + if(glob_opts_init.nz > 0) real_t_vctrs.push_back(&wp); + } + else if(glob_opts_init.turb_cond_switch) + if(glob_opts_init.nz > 0) real_t_vctrs.push_back(&wp); + + if(glob_opts_init.turb_cond_switch) + if(glob_opts_init.nz > 0) real_t_vctrs.push_back(&ssp); + const int real_vctrs_count = real_t_vctrs.size(); assert(out_real_bfr.size() >= lft_count * real_vctrs_count); assert(in_real_bfr.size() >= lft_count * real_vctrs_count); diff --git a/src/particles.tpp b/src/particles.tpp index 4ddd40cb8..df6c0435a 100644 --- a/src/particles.tpp +++ b/src/particles.tpp @@ -69,6 +69,9 @@ #include "impl/particles_impl_hskpng_ijk.ipp" #include "impl/particles_impl_hskpng_Tpr.ipp" #include "impl/particles_impl_hskpng_vterm.ipp" +#include "impl/particles_impl_hskpng_turb_vel.ipp" +#include "impl/particles_impl_hskpng_turb_ss.ipp" +#include "impl/particles_impl_hskpng_tke.ipp" #include "impl/particles_impl_hskpng_sort.ipp" #include "impl/particles_impl_hskpng_count.ipp" #include "impl/particles_impl_hskpng_remove.ipp" @@ -79,6 +82,7 @@ #include "impl/particles_impl_sync.ipp" #include "impl/particles_impl_bcnd.ipp" // bcnd has to be b4 adve for periodic struct; move it to separate file in detail... #include "impl/particles_impl_adve.ipp" +#include "impl/particles_impl_turb_adve.ipp" #include "impl/particles_impl_cond_common.ipp" #include "impl/particles_impl_cond.ipp" #include "impl/particles_impl_cond_sstp.ipp" diff --git a/src/particles_multi_gpu_step.ipp b/src/particles_multi_gpu_step.ipp index b7e879122..d4f68b9f3 100644 --- a/src/particles_multi_gpu_step.ipp +++ b/src/particles_multi_gpu_step.ipp @@ -22,10 +22,11 @@ namespace libcloudphxx const arrinfo_t courant_1, const arrinfo_t courant_2, const arrinfo_t courant_3, + const arrinfo_t diss_rate, 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, diss_rate, ambient_chem); } template @@ -36,10 +37,11 @@ namespace libcloudphxx const arrinfo_t courant_1, const arrinfo_t courant_2, const arrinfo_t courant_3, + const arrinfo_t diss_rate, std::map > ambient_chem ) { - pimpl->mcuda_run(&particles_t::sync_in, th, rv, rhod, courant_1, courant_2, courant_3, ambient_chem); + pimpl->mcuda_run(&particles_t::sync_in, th, rv, rhod, courant_1, courant_2, courant_3, diss_rate, ambient_chem); } template diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 6919082e3..3d6790e48 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -20,10 +20,11 @@ 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 diss_rate, // defaults to NULL-NULL pair std::map > ambient_chem ) { - sync_in(th, rv, rhod, courant_x, courant_y, courant_z, ambient_chem); + sync_in(th, rv, rhod, courant_x, courant_y, courant_z, diss_rate, ambient_chem); step_cond(opts, th, rv, ambient_chem); } @@ -35,6 +36,7 @@ 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 diss_rate, // defaults to NULL-NULL pair std::map > ambient_chem ) { @@ -68,6 +70,12 @@ namespace libcloudphxx if (!pimpl->opts_init.chem_switch && ambient_chem.size() != 0) throw std::runtime_error("chemistry was switched off and ambient_chem is not empty"); + + if ( (pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_switch || pimpl->opts_init.turb_coal_switch) && diss_rate.is_null()) + throw std::runtime_error("turbulent advection, coalescence and condesation are not switched off and diss_rate is empty"); + + if ( !(pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_switch || pimpl->opts_init.turb_coal_switch) && !diss_rate.is_null()) + throw std::runtime_error("turbulent advection, coalescence and condesation are switched off and diss_rate is not empty"); // if (pimpl->l2e[&pimpl->courant_x].size() == 0) // TODO: y, z,... @@ -82,6 +90,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 (pimpl->l2e[&pimpl->diss_rate].size() == 0) + if (!diss_rate.is_null()) pimpl->init_e2l(diss_rate, &pimpl->diss_rate); + // did rhod change pimpl->var_rho = !rhod.is_null(); @@ -91,6 +102,7 @@ namespace libcloudphxx pimpl->sync(courant_x, pimpl->courant_x); pimpl->sync(courant_y, pimpl->courant_y); pimpl->sync(courant_z, pimpl->courant_z); + pimpl->sync(diss_rate, pimpl->diss_rate); pimpl->sync(rhod, pimpl->rhod); nancheck(pimpl->th, " th after sync-in"); @@ -98,11 +110,16 @@ namespace libcloudphxx 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->diss_rate, " diss_rate after sync-in"); nancheck(pimpl->rhod, " rhod after sync-in"); + if(pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_switch || pimpl->opts_init.turb_coal_switch) + {nancheck(pimpl->diss_rate, " diss_rate 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); + if(pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_switch || pimpl->opts_init.turb_coal_switch) + {assert(*thrust::min_element(pimpl->diss_rate.begin(), pimpl->diss_rate.end()) >= 0);} // check if courants are greater than 2 since it would break the predictor-corrector (halo of size 2 in the x direction) if(!(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(-2.) )) )) @@ -144,6 +161,9 @@ namespace libcloudphxx if (!pimpl->should_now_run_cond) throw std::runtime_error("please call sync_in() before calling step_cond()"); + if(opts.turb_cond && !pimpl->opts_init.turb_cond_switch) + throw std::runtime_error("turb_cond_swtich=False, but turb_cond==True"); + pimpl->should_now_run_cond = false; // condensation/evaporation @@ -155,7 +175,9 @@ namespace libcloudphxx for (int step = 0; step < pimpl->opts_init.sstp_cond; ++step) { pimpl->sstp_step_exact(step); - pimpl->cond_sstp(pimpl->opts_init.dt / pimpl->opts_init.sstp_cond, opts.RH_max); + if(opts.turb_cond) + pimpl->sstp_step_ssp(pimpl->opts_init.dt / pimpl->opts_init.sstp_cond); + pimpl->cond_sstp(pimpl->opts_init.dt / pimpl->opts_init.sstp_cond, opts.RH_max, opts.turb_cond); } // copy sstp_tmp_rv and th to rv and th pimpl->update_state(pimpl->rv, pimpl->sstp_tmp_rv); @@ -167,8 +189,10 @@ namespace libcloudphxx for (int step = 0; step < pimpl->opts_init.sstp_cond; ++step) { pimpl->sstp_step(step); + if(opts.turb_cond) + pimpl->sstp_step_ssp(pimpl->opts_init.dt / pimpl->opts_init.sstp_cond); pimpl->hskpng_Tpr(); - pimpl->cond(pimpl->opts_init.dt / pimpl->opts_init.sstp_cond, opts.RH_max); + pimpl->cond(pimpl->opts_init.dt / pimpl->opts_init.sstp_cond, opts.RH_max, opts.turb_cond); } } @@ -289,6 +313,12 @@ namespace libcloudphxx if(opts.sedi && !pimpl->opts_init.sedi_switch) throw std::runtime_error("all sedimentation was switched off in opts_init"); + if(opts.turb_adve && !pimpl->opts_init.turb_adve_switch) + throw std::runtime_error("turb_adve_switch=False, but turb_adve==True"); + + if(opts.turb_adve && pimpl->n_dims==0) + throw std::runtime_error("turbulent advection does not work in 0D"); + if (opts.chem_dsl) { // saving rv to be used as rv_old @@ -309,7 +339,7 @@ namespace libcloudphxx for (int step = 0; step < pimpl->opts_init.sstp_coal; ++step) { // collide - pimpl->coal(pimpl->opts_init.dt / pimpl->opts_init.sstp_coal); + pimpl->coal(pimpl->opts_init.dt / pimpl->opts_init.sstp_coal, opts.turb_coal); // update invalid vterm if (step + 1 != pimpl->opts_init.sstp_coal) @@ -325,15 +355,40 @@ namespace libcloudphxx } } + if (opts.turb_adve || opts.turb_cond) + { + // calc tke (diss_rate now holds TKE, not dissipation rate! Hence this must be done after coal, which requires diss rate) + pimpl->hskpng_tke(); + } + if (opts.turb_adve) + { + // calc turbulent perturbation of velocity + pimpl->hskpng_turb_vel(); + } + else if (opts.turb_cond) + { + // calc turbulent perturbation only of vertical velocity + pimpl->hskpng_turb_vel(true); + } + + if(opts.turb_cond) + { + // calculate the time derivatie of the turbulent supersaturation perturbation; applied in the next step during condensation substepping - is the delay a problem? + pimpl->hskpng_turb_dot_ss(); + } + // advection, it invalidates i,j,k and ijk! if (opts.adve) pimpl->adve(); // revert to the desired adve scheme (in case we used eulerian this timestep for halo reasons) pimpl->adve_scheme = pimpl->opts_init.adve_scheme; + // apply turbulent perturbation of velocity, TODO: add it to advection velocity (turb_vel_calc would need to be called couple times in the pred-corr advection + diss_rate would need a halo) + if (opts.turb_adve) pimpl->turb_adve(); + // sedimentation has to be done after advection, so that negative z doesnt crash hskpng_ijk in adve if (opts.sedi) { - // advection with terminal velocity + // advection with terminal velocity, TODO: add it to the advection velocity (makes a difference for predictor-corrector) pimpl->sedi(); } diff --git a/tests/python/physics/coalescence_onishi_hall.py b/tests/python/physics/coalescence_onishi_hall.py index 80cce74c2..eb97f1955 100755 --- a/tests/python/physics/coalescence_onishi_hall.py +++ b/tests/python/physics/coalescence_onishi_hall.py @@ -43,6 +43,7 @@ def expvolumelnr(lnr): th = 300 * np.ones((1,)) rv = 1. * np.ones((1,)) rhod = 1.22419 * np.ones((1,)) +diss_rate = epsilon * np.ones((1,)) opts_init.dry_distros = {0.:expvolumelnr} @@ -58,7 +59,7 @@ def expvolumelnr(lnr): opts_init.kernel = lgrngn.kernel_t.hall else: opts_init.kernel = lgrngn.kernel_t.onishi_hall - opts_init.kernel_parameters = np.array([epsilon, re_lambda]) + opts_init.kernel_parameters = np.array([re_lambda]) for zz in range(n_runs): opts_init.rng_seed = int(time.time()) @@ -78,7 +79,10 @@ def expvolumelnr(lnr): t=0 t10=0 while(t10 == 0): - prtcls.step_sync(Opts,th,rv,rhod) + if(kernel_t == 0): + prtcls.step_sync(Opts,th,rv,rhod) + else: + prtcls.step_sync(Opts,th,rv,rhod, diss_rate = diss_rate) prtcls.step_async(Opts) rain_volume = diag_rain_volume() t+=opts_init.dt diff --git a/tests/python/unit/api_lgrngn.py b/tests/python/unit/api_lgrngn.py index ba36ed047..96d7b0f86 100644 --- a/tests/python/unit/api_lgrngn.py +++ b/tests/python/unit/api_lgrngn.py @@ -347,6 +347,28 @@ def lognormal(lnr): opts_init.dry_sizes = dict() +# ---------- +# 0D turb +print "0D turb" +eps = arr_t([ 1e-4]) +opts_init.turb_cond_switch=True +opts.turb_cond=True +prtcls = lgrngn.factory(backend, opts_init) +prtcls.init(th, rv, rhod) +prtcls.step_sync(opts, th, rv, diss_rate=eps) +prtcls.step_async(opts) + +prtcls.diag_all() +prtcls.diag_sd_conc() +assert len(frombuffer(prtcls.outbuf())) == 1 +print frombuffer(prtcls.outbuf()) +assert (frombuffer(prtcls.outbuf()) > 0).all() +assert sum(frombuffer(prtcls.outbuf())) == 64 + +opts_init.turb_cond_switch=False +opts.turb_cond=False + + # ---------- # 1D (periodic horizontal domain) @@ -386,7 +408,28 @@ def lognormal(lnr): assert (frombuffer(prtcls.outbuf()) > 0).all() assert sum(frombuffer(prtcls.outbuf())) == opts_init.nx * opts_init.sd_conc +print "1D turb" +eps = arr_t([ 1e-4, 1e-4, 1e-4]) +opts_init.turb_adve_switch=True +opts.turb_adve=True +opts_init.turb_cond_switch=True +opts.turb_cond=True +prtcls = lgrngn.factory(backend, opts_init) +prtcls.init(th, rv, rhod) +prtcls.step_sync(opts, th, rv, diss_rate=eps) +prtcls.step_async(opts) + +prtcls.diag_all() +prtcls.diag_sd_conc() +assert len(frombuffer(prtcls.outbuf())) == opts_init.nx +print frombuffer(prtcls.outbuf()) +assert (frombuffer(prtcls.outbuf()) > 0).all() +assert sum(frombuffer(prtcls.outbuf())) == opts_init.nx * opts_init.sd_conc +opts_init.turb_adve_switch=False +opts.turb_adve=False +opts_init.turb_cond_switch=False +opts.turb_cond=False # ---------- # 2D (periodic horizontal domain) @@ -436,6 +479,29 @@ def lognormal(lnr): #TODO: test profile vs. 2D array +print "2D turb" +eps = arr_t([[ 1e-4, 1e-4], [ 1e-4, 1e-4]]) +opts_init.turb_adve_switch=True +opts.turb_adve=True +opts_init.turb_cond_switch=True +opts.turb_cond=True +prtcls = lgrngn.factory(backend, opts_init) +prtcls.init(th, rv, rhod) +prtcls.step_sync(opts, th, rv, diss_rate=eps) +prtcls.step_async(opts) + +prtcls.diag_all() +prtcls.diag_sd_conc() +assert len(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx +print frombuffer(prtcls.outbuf()) +assert (frombuffer(prtcls.outbuf()) > 0).all() +assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.sd_conc +assert opts_init.nx == prtcls.opts_init.nx + +opts_init.turb_adve_switch=False +opts.turb_adve=False +opts_init.turb_cond_switch=False +opts.turb_cond=False # ---------- @@ -473,6 +539,29 @@ def lognormal(lnr): assert (frombuffer(prtcls.outbuf()) > 0).all() assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny * opts_init.sd_conc +print "3D turb" +eps = arr_t([eps, eps]) +opts_init.turb_adve_switch=True +opts.turb_adve=False +opts_init.turb_cond_switch=True +opts.turb_cond=True +prtcls = lgrngn.factory(backend, opts_init) +prtcls.init(th, rv, rhod) +prtcls.step_sync(opts, th, rv, diss_rate=eps) +prtcls.step_async(opts) + +prtcls.diag_all() +prtcls.diag_sd_conc() +assert len(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny +print frombuffer(prtcls.outbuf()) +assert (frombuffer(prtcls.outbuf()) > 0).all() +assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny * opts_init.sd_conc + +opts_init.turb_adve_switch=False +opts.turb_adve=False +opts_init.turb_cond_switch=False +opts.turb_cond=False + # 3D with large_tail option print "3D large tail" opts_init.sd_conc_large_tail = 1 diff --git a/tests/python/unit/col_kernels.py b/tests/python/unit/col_kernels.py index 41863e592..f76b305a1 100644 --- a/tests/python/unit/col_kernels.py +++ b/tests/python/unit/col_kernels.py @@ -8,6 +8,7 @@ rhod = 1. * np.ones((1,)) th = 300. * np.ones((1,)) rv = 0.01 * np.ones((1,)) +diss_rate = 0.04 * np.ones((1,)) def lognormal(lnr): mean_r = .04e-6 / 2 @@ -31,7 +32,8 @@ def lognormal(lnr): opts_init.kernel = kernel opts_init.kernel_parameters = np.array([]) if(kernel == lgrngn.kernel_t.onishi_hall_davis_no_waals or kernel == lgrngn.kernel_t.onishi_hall): - opts_init.kernel_parameters = np.array([0.04, 100]); + opts_init.turb_coal_switch = True + opts_init.kernel_parameters = np.array([100.]); if(kernel == lgrngn.kernel_t.golovin): opts_init.kernel_parameters = np.array([1.]) if(kernel == lgrngn.kernel_t.geometric): @@ -47,6 +49,8 @@ def lognormal(lnr): except: prtcls = lgrngn.factory(lgrngn.backend_t.serial, opts_init) + print opts_init.turb_coal_switch + prtcls.init(th, rv, rhod) Opts = lgrngn.opts_t() @@ -58,5 +62,10 @@ def lognormal(lnr): Opts.chem_dsc = False Opts.chem_rct = False - prtcls.step_sync(Opts,th,rv,rhod) + if(kernel == lgrngn.kernel_t.onishi_hall_davis_no_waals or kernel == lgrngn.kernel_t.onishi_hall): + Opts.turb_coal = True + prtcls.step_sync(Opts,th,rv,rhod, diss_rate = diss_rate) + else: + prtcls.step_sync(Opts,th,rv,rhod) + prtcls.step_async(Opts) diff --git a/tests/python/unit/lgrngn_subsidence.py b/tests/python/unit/lgrngn_subsidence.py index c2e3d9860..f2988aa95 100644 --- a/tests/python/unit/lgrngn_subsidence.py +++ b/tests/python/unit/lgrngn_subsidence.py @@ -80,7 +80,7 @@ def lognormal(lnr): assert(tab_out[0][3] == 0.) tolerance = 1. / sqrt(Opts_init.sd_conc) print "relative tolerance: ", tolerance -assert np.isclose(res_2, tab_out[0][2], atol=0., rtol=4*tolerance) +assert np.isclose(res_2, tab_out[0][2], atol=0., rtol=5*tolerance) assert np.isclose(res_01, tab_out[0][1], atol=0., rtol=tolerance) assert np.isclose(res_01, tab_out[0][0], atol=0., rtol=tolerance) diff --git a/tests/python/unit/lgrngn_turb_adve.py b/tests/python/unit/lgrngn_turb_adve.py new file mode 100644 index 000000000..975abc869 --- /dev/null +++ b/tests/python/unit/lgrngn_turb_adve.py @@ -0,0 +1,71 @@ +import sys +sys.path.insert(0, "../../bindings/python/") + +from libcloudphxx import lgrngn + +import numpy as np +from math import exp, log, sqrt, pi +from time import time + +def lognormal(lnr): + mean_r = .04e-6 / 2 + stdev = 1.4 + n_tot = 60e6 + return n_tot * exp( + -pow((lnr - log(mean_r)), 2) / 2 / pow(log(stdev),2) + ) / log(stdev) / sqrt(2*pi); + +Opts_init = lgrngn.opts_init_t() +kappa = .61 +Opts_init.dry_distros = {kappa:lognormal} +Opts_init.coal_switch = False +Opts_init.sedi_switch = False +Opts_init.turb_adve_switch = True + +Opts_init.dt = 1 + +Opts_init.nz = 6 +Opts_init.nx = 6 +Opts_init.dz = 1 +Opts_init.dx = 1 +Opts_init.z1 = Opts_init.nz * Opts_init.dz +Opts_init.x1 = Opts_init.nx * Opts_init.dx + +Opts_init.rng_seed = int(time()) +Opts_init.sd_conc = 100 +Opts_init.n_sd_max = Opts_init.sd_conc * (Opts_init.nx * Opts_init.nz) + +Backend = lgrngn.backend_t.serial + +Opts = lgrngn.opts_t() +Opts.adve = False +Opts.turb_adve = True +Opts.sedi = False +Opts.cond = False +Opts.coal = False +Opts.chem = False +Opts.rcyc = False + +Rhod = 1. * np.ones((Opts_init.nx, Opts_init.nz)) +Th = 300. * np.ones((Opts_init.nx, Opts_init.nz)) +Rv = 0.01 * np.ones((Opts_init.nx, Opts_init.nz)) +diss_rate = 1e-4 * np.ones((Opts_init.nx, Opts_init.nz)) + +prtcls = lgrngn.factory(Backend, Opts_init) +prtcls.init(Th, Rv, Rhod) + +prtcls.diag_all() +prtcls.diag_sd_conc() +tab_in = np.copy(np.frombuffer(prtcls.outbuf()).reshape(Opts_init.nx, Opts_init.nz)) +print "at init \n", tab_in + +for it in range(100): + prtcls.step_sync(Opts, Th, Rv, Rhod, diss_rate = diss_rate) + prtcls.step_async(Opts) + +prtcls.diag_all() +prtcls.diag_sd_conc() +tab_out = np.copy(np.frombuffer(prtcls.outbuf()).reshape(Opts_init.nx, Opts_init.nz)) +print "after 100s \n", tab_out + +assert(np.array_equal(tab_in,tab_out) == False) # turbulence should have moved SDs around diff --git a/tests/python/unit/sstp_cond.py b/tests/python/unit/sstp_cond.py index 5f10f4f51..777f060c6 100644 --- a/tests/python/unit/sstp_cond.py +++ b/tests/python/unit/sstp_cond.py @@ -3,7 +3,7 @@ from libcloudphxx import lgrngn -from numpy import array as arr_t, frombuffer, repeat, zeros, float64, ones, roll, pi, allclose +from numpy import array as arr_t, frombuffer, repeat, zeros, float64, ones, roll, pi, allclose, copy from math import exp, log, sqrt, pi @@ -15,89 +15,108 @@ def lognormal(lnr): -pow((lnr - log(mean_r)), 2) / 2 / pow(log(stdev),2) ) / log(stdev) / sqrt(2*pi); -opts_init = lgrngn.opts_init_t() -kappa = .61 -opts_init.dry_distros = {kappa:lognormal} -opts_init.coal_switch=0 -opts_init.sedi_switch=0 -opts_init.dt = 1 -opts_init.sd_conc = 64 -opts_init.n_sd_max = 512 -opts_init.rng_seed = 396 -opts_init.exact_sstp_cond = True # test would fail with per-cell sstp logic -spinup = 20 - -backend = lgrngn.backend_t.serial - -opts = lgrngn.opts_t() -opts.sedi=0 -opts.coal=0 -opts.cond=1 - -# 1D (periodic horizontal domain) -rhod = arr_t([ 1., 1.]) -C = arr_t([ 1., 1., 1.]) - -opts_init.nx = 2 -opts_init.dx = 1 -opts_init.x1 = opts_init.nx * opts_init.dx - -for sstp_cond in [1,2,5]: - print 'sstp_cond = ' + str(sstp_cond) - opts.adve=0 - opts_init.sstp_cond = sstp_cond - prtcls = lgrngn.factory(backend, opts_init) - th = arr_t([300., 300.]) - rv = arr_t([ .0025, .0095]) # first cell subsaturated, second cell supersaturated - prtcls.init(th, rv, rhod, Cx=C) - - #equilibrium wet moment post spinup - prtcls.diag_all() - prtcls.diag_wet_mom(3); - wet_post_init = frombuffer(prtcls.outbuf()).copy() - water_post_init = 1000. * 4./3. * pi * wet_post_init + rv +def test(turb_cond): + print 'turb_cond = ', turb_cond + opts_init = lgrngn.opts_init_t() + kappa = .61 + opts_init.dry_distros = {kappa:lognormal} + opts_init.coal_switch=0 + opts_init.sedi_switch=0 + opts_init.dt = 1 + opts_init.sd_conc = 64 + opts_init.n_sd_max = 512 + opts_init.rng_seed = 396 + opts_init.exact_sstp_cond = True # test would fail with per-cell sstp logic + opts_init.turb_cond_switch = turb_cond + spinup = 20 - #spinup to get equilibrium - prtcls.step_sync(opts, th, rv) - for it in range(spinup): - prtcls.step_async(opts) - prtcls.step_sync(opts, th, rv) + backend = lgrngn.backend_t.serial - #equilibrium wet moment post spinup - prtcls.diag_all() - prtcls.diag_wet_mom(3); - wet_post_spin = frombuffer(prtcls.outbuf()).copy() - water_post_spin = 1000. * 4./3. * pi * wet_post_spin + rv - assert allclose(water_post_spin, water_post_init, atol=0, rtol=1e-10) #some discrepancy due to water density + opts = lgrngn.opts_t() + opts.sedi=0 + opts.coal=0 + opts.cond=1 + opts.turb_cond=turb_cond - #advect SDs - opts.adve=1 - prtcls.step_async(opts) + opts_init.nx = 2 + opts_init.dx = 1 + opts_init.x1 = opts_init.nx * opts_init.dx - #equilibrium wet moment post adve - prtcls.diag_all() - prtcls.diag_wet_mom(3); - wet_post_adve = frombuffer(prtcls.outbuf()).copy() - wet_post_adve_roll = roll(wet_post_adve,1).copy() - assert all(wet_post_adve_roll == wet_post_spin) + opts_init.nz = 1 + opts_init.dz = 1 + opts_init.z1 = opts_init.nz * opts_init.dz - #advect rv - tmp = rv.copy() - rv[0] = tmp[1] - rv[1] = tmp[0] - #advect th - tmp = th.copy() - th[0] = tmp[1] - th[1] = tmp[0] - - #condensation with advected SDs and rv - prtcls.step_sync(opts, th, rv) + rhod = 1. * ones((opts_init.nx, opts_init.nz)) + Cx = 1. * ones((opts_init.nx+1, opts_init.nz)) + Cz = 0. * ones((opts_init.nx, opts_init.nz+1)) + eps = 1e-4 * ones((opts_init.nx, opts_init.nz)) + + for sstp_cond in [1,2,5]: + print 'sstp_cond = ' + str(sstp_cond) + opts.adve=0 + opts_init.sstp_cond = sstp_cond + prtcls = lgrngn.factory(backend, opts_init) + th = 300 * ones((opts_init.nx, opts_init.nz)) + rv = .0025 * ones((opts_init.nx, opts_init.nz)) + rv[1,0] = 0.0095 + prtcls.init(th, rv, rhod, Cx=Cx, Cz=Cz) + + #equilibrium wet moment post spinup + prtcls.diag_all() + prtcls.diag_wet_mom(3); + wet_post_init = copy(frombuffer(prtcls.outbuf()).reshape(opts_init.nx, opts_init.nz)) + water_post_init = 1000. * 4./3. * pi * wet_post_init + rv + + #spinup to get equilibrium + if(turb_cond): + prtcls.step_sync(opts, th, rv, diss_rate=eps) + else: + prtcls.step_sync(opts, th, rv) + for it in range(spinup): + prtcls.step_async(opts) + if(turb_cond): + prtcls.step_sync(opts, th, rv, diss_rate=eps) + else: + prtcls.step_sync(opts, th, rv) + + #equilibrium wet moment post spinup + prtcls.diag_all() + prtcls.diag_wet_mom(3); + wet_post_spin = copy(frombuffer(prtcls.outbuf()).reshape(opts_init.nx, opts_init.nz)) + water_post_spin = 1000. * 4./3. * pi * wet_post_spin + rv + assert allclose(water_post_spin, water_post_init, atol=0, rtol=1e-10) #some discrepancy due to water density + + #advect SDs + opts.adve=1 + prtcls.step_async(opts) + + #equilibrium wet moment post adve + prtcls.diag_all() + prtcls.diag_wet_mom(3); + wet_post_adve = copy(frombuffer(prtcls.outbuf()).reshape(opts_init.nx, opts_init.nz)) + wet_post_adve_roll = roll(wet_post_adve,1).copy() + assert all(wet_post_adve_roll == wet_post_spin) + + #advect rv + tmp = rv.copy() + rv[0,0] = tmp[1,0] + rv[1,0] = tmp[0,0] + #advect th + tmp = th.copy() + th[0,0] = tmp[1,0] + th[1,0] = tmp[0,0] - #wet mom post adve and cond - prtcls.diag_all() - prtcls.diag_wet_mom(3); - wet_post_adve_cond = frombuffer(prtcls.outbuf()).copy() - print wet_post_adve - print wet_post_adve_cond - assert allclose(wet_post_adve, wet_post_adve_cond, atol=0, rtol=3e-2) + #condensation with advected SDs and rv + if(turb_cond): + prtcls.step_sync(opts, th, rv, diss_rate=eps) + else: + prtcls.step_sync(opts, th, rv) + + #wet mom post adve and cond + prtcls.diag_all() + prtcls.diag_wet_mom(3); + wet_post_adve_cond = copy(frombuffer(prtcls.outbuf()).reshape(opts_init.nx, opts_init.nz)) + assert allclose(wet_post_adve, wet_post_adve_cond, atol=0, rtol=3e-2) +test(False) +test(True)