From ec2f6e5982eefd1709b93206b362fb0670289fbe Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 27 Sep 2018 11:28:05 +0200 Subject: [PATCH 01/44] add turb switch, turb vel arrays, diss_rate array --- include/libcloudph++/lgrngn/opts_init.hpp | 4 ++- include/libcloudph++/lgrngn/particles.hpp | 6 +++++ src/impl/particles_impl.ipp | 13 +++++++++- src/impl/particles_impl_hskpng_remove.ipp | 7 ++++++ src/impl/particles_impl_hskpng_resize.ipp | 7 ++++++ src/impl/particles_impl_init_hskpng_npart.ipp | 7 ++++++ src/impl/particles_impl_init_sync.ipp | 9 ++++--- src/particles_step.ipp | 25 +++++++++++++++++-- 8 files changed, 71 insertions(+), 7 deletions(-) diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index c36d075a1..6132119dc 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -103,7 +103,8 @@ 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_switch; // if true, turbulent motion of SDs is modeled int sstp_chem; real_t chem_rho; @@ -142,6 +143,7 @@ namespace libcloudphxx coal_switch(true), // coalescence turned on by default src_switch(false), // source turned off by default exact_sstp_cond(false), + turb_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), 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/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 6760a78ac..a942bc6f0 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -76,6 +76,9 @@ 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 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 @@ -136,7 +139,10 @@ namespace libcloudphxx T, // temperature [K] p, // pressure [Pa] RH, // relative humisity - eta;// dynamic viscosity + eta,// dynamic viscosity + TKE;// turbulent kinetic energy + + 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; @@ -276,6 +282,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 ? opts_init.dx: // 1D + n_dims == 2 ? common::turbulence::length_scale(opts_init.dx, opts_init.dz): // 2D + common::turbulence::length_scale(opts_init.dx, opts_init.dy, opts_init.dz) // 3D + ), pure_const_multi (((opts_init.sd_conc) == 0) && (opts_init.sd_const_multi > 0 || opts_init.sd_const_multi_dry_sizes > 0)) // coal prob can be greater than one only in sd_conc simulations { // note: there could be less tmp data spaces if _cell vectors diff --git a/src/impl/particles_impl_hskpng_remove.ipp b/src/impl/particles_impl_hskpng_remove.ipp index c6e2a7b79..f4292816c 100644 --- a/src/impl/particles_impl_hskpng_remove.ipp +++ b/src/impl/particles_impl_hskpng_remove.ipp @@ -32,6 +32,13 @@ namespace libcloudphxx 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.turb_switch) + { + if (opts_init.nx != 0) real_t_vctrs.push_back(&up); + if (opts_init.ny != 0) real_t_vctrs.push_back(&vp); + if (opts_init.nz != 0) real_t_vctrs.push_back(&wp); + } + if(opts_init.sstp_cond>1 && opts_init.exact_sstp_cond) { real_t_vctrs.push_back(&sstp_tmp_th); diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index d12e85eed..bf57851c0 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -33,6 +33,13 @@ namespace libcloudphxx if (opts_init.ny != 0) y.resize(n_part); if (opts_init.nz != 0) z.resize(n_part); + if(opts_init.turb_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.chem_switch || opts_init.sstp_cond > 1 || n_dims >= 2) { tmp_device_real_part1.resize(n_part); diff --git a/src/impl/particles_impl_init_hskpng_npart.ipp b/src/impl/particles_impl_init_hskpng_npart.ipp index d06a74d87..f92b4c4a4 100644 --- a/src/impl/particles_impl_init_hskpng_npart.ipp +++ b/src/impl/particles_impl_init_hskpng_npart.ipp @@ -23,6 +23,13 @@ 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); + if(opts_init.turb_switch) + { + if (opts_init.nx != 0) up.reserve(opts_init.n_sd_max); + if (opts_init.ny != 0) vp.reserve(opts_init.n_sd_max); + if (opts_init.nz != 0) wp.reserve(opts_init.n_sd_max); + } + 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); diff --git a/src/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index ca3e0cc8d..191270c82 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_swtich) + for (int i = 0; i < chem_gas_n; ++i) + ambient_chem[(chem_species_t)i].resize(n_cell); + if(opts_init.turb_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/particles_step.ipp b/src/particles_step.ipp index 96c5ab8af..87c67257d 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_switch && diss_rate.is_null()) + throw std::runtime_error("turbulence was not switched off and diss_rate is empty"); + + if (!pimpl->opts_init.turb_switch && !diss_rate.is_null()) + throw std::runtime_error("turbulence was 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) + 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,14 @@ 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"); 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(opts_init.turb_switch) + assert(*thrust::min_element(pimpl->rhod.begin(), pimpl->rhod.end()) >= 0); // check if courants are greater than 2 since it would break the predictor-corrector (halo of size 2 in the x direction) assert(pimpl->opts_init.adve_scheme != as_t::pred_corr || (courant_x.is_null() || ((*(thrust::min_element(pimpl->courant_x.begin(), pimpl->courant_x.end()))) >= real_t(-2.) )) ); @@ -313,13 +328,19 @@ namespace libcloudphxx } } + // calc turbulent perturbation of velocity + if (opts.adve && opts.turb_switch) pimpl->hskpng_turb_vel_calc(); + // advection, it invalidates i,j,k and ijk! if (opts.adve) pimpl->adve(); + // 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.adve && opts.turb_switch) 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(); } From c1d2eccafd8205677691eb8baf631bbc0147e63f Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 27 Sep 2018 11:42:50 +0200 Subject: [PATCH 02/44] up,vp,wp copying during rcyc and distmem mcuda --- include/libcloudph++/lgrngn/opts.hpp | 3 ++- src/impl/particles_impl_init_hskpng_npart.ipp | 6 +++--- src/impl/particles_impl_rcyc.ipp | 7 +++++++ .../particles_multi_gpu_impl_step_async_and_copy.ipp | 9 +++++++++ src/particles_step.ipp | 7 +++++-- 5 files changed, 26 insertions(+), 6 deletions(-) diff --git a/include/libcloudph++/lgrngn/opts.hpp b/include/libcloudph++/lgrngn/opts.hpp index 58849bed1..fd26b012e 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; // 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), RH_max(44) // :) (anything greater than 1.1 would be enough { } diff --git a/src/impl/particles_impl_init_hskpng_npart.ipp b/src/impl/particles_impl_init_hskpng_npart.ipp index f92b4c4a4..f76b37937 100644 --- a/src/impl/particles_impl_init_hskpng_npart.ipp +++ b/src/impl/particles_impl_init_hskpng_npart.ipp @@ -25,9 +25,9 @@ namespace libcloudphxx if(opts_init.turb_switch) { - if (opts_init.nx != 0) up.reserve(opts_init.n_sd_max); - if (opts_init.ny != 0) vp.reserve(opts_init.n_sd_max); - if (opts_init.nz != 0) wp.reserve(opts_init.n_sd_max); + 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.); } vt.resize(opts_init.n_sd_max, 0.); // so that it may be safely used in condensation before first update diff --git a/src/impl/particles_impl_rcyc.ipp b/src/impl/particles_impl_rcyc.ipp index dcc945294..6ec24cb17 100644 --- a/src/impl/particles_impl_rcyc.ipp +++ b/src/impl/particles_impl_rcyc.ipp @@ -112,6 +112,13 @@ namespace libcloudphxx for (int i = 0; i < chem_all; ++i) detail::copy_prop(chem_bgn[i], sorted_id, n_flagged); } + + if(opts_init.turb_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); + } { namespace arg = thrust::placeholders; 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..2d19a0861 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,9 @@ 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 &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 +138,12 @@ namespace libcloudphxx if(particles[dev_id]->pimpl->const_p) real_t_vctrs.push_back(&sstp_tmp_p); } + if(glob_opts_init.turb_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); + } 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_step.ipp b/src/particles_step.ipp index 87c67257d..48a53bf3e 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -292,6 +292,9 @@ 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_switch) + throw std::runtime_error("all turbulence sgs was switched off in opts_init, but turb_adve==True"); + if (opts.chem_dsl) { // saving rv to be used as rv_old @@ -329,13 +332,13 @@ namespace libcloudphxx } // calc turbulent perturbation of velocity - if (opts.adve && opts.turb_switch) pimpl->hskpng_turb_vel_calc(); + if (opts.turb_adve) pimpl->hskpng_turb_vel_calc(); // advection, it invalidates i,j,k and ijk! if (opts.adve) pimpl->adve(); // 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.adve && opts.turb_switch) pimpl->turb_adve(); + 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) From 8dc2bdfc4929956c3eec7742d54707cdb6810092 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 27 Sep 2018 14:13:56 +0200 Subject: [PATCH 03/44] working tk ecalc --- src/impl/particles_impl.ipp | 12 ++++++++---- src/impl/particles_impl_hskpng_resize.ipp | 2 +- src/impl/particles_impl_init_sync.ipp | 2 +- src/particles_step.ipp | 13 +++++++++---- 4 files changed, 19 insertions(+), 10 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index a942bc6f0..8bbd13569 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -14,6 +14,8 @@ #include #include +#include + #include namespace libcloudphxx @@ -140,7 +142,7 @@ namespace libcloudphxx p, // pressure [Pa] RH, // relative humisity eta,// dynamic viscosity - TKE;// turbulent kinetic energy + 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 @@ -283,9 +285,9 @@ namespace libcloudphxx halo_size * (opts_init.nz + 1) * opts_init.ny // 3D ), L( - n_dims == 1 ? opts_init.dx: // 1D - n_dims == 2 ? common::turbulence::length_scale(opts_init.dx, opts_init.dz): // 2D - common::turbulence::length_scale(opts_init.dx, opts_init.dy, opts_init.dz) // 3D + 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 ), pure_const_multi (((opts_init.sd_conc) == 0) && (opts_init.sd_const_multi > 0 || opts_init.sd_const_multi_dry_sizes > 0)) // coal prob can be greater than one only in sd_conc simulations { @@ -404,6 +406,8 @@ namespace libcloudphxx void hskpng_vterm_all(); void hskpng_vterm_invalid(); + void hskpng_tke(); + void hskpng_turb_vel(); void hskpng_remove_n0(); void hskpng_resize_npart(); diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index bf57851c0..6c09d7bde 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -33,7 +33,7 @@ namespace libcloudphxx if (opts_init.ny != 0) y.resize(n_part); if (opts_init.nz != 0) z.resize(n_part); - if(opts_init.turb_switch) + if(opts_init.turb_switch) // TODO: they are not needed by turb_coal, but diss_rate is - make more switches { if (opts_init.nx != 0) up.resize(n_part); if (opts_init.ny != 0) vp.resize(n_part); diff --git a/src/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index 191270c82..a19c6d6fc 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/src/impl/particles_impl_init_sync.ipp @@ -17,7 +17,7 @@ namespace libcloudphxx p.resize(n_cell); th.resize(n_cell); rv.resize(n_cell); - if(opts_init.chem_swtich) + 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_switch) diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 48a53bf3e..b4d48918a 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -116,7 +116,7 @@ namespace libcloudphxx 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(opts_init.turb_switch) + if(pimpl->opts_init.turb_switch) assert(*thrust::min_element(pimpl->rhod.begin(), pimpl->rhod.end()) >= 0); // check if courants are greater than 2 since it would break the predictor-corrector (halo of size 2 in the x direction) @@ -331,14 +331,19 @@ namespace libcloudphxx } } - // calc turbulent perturbation of velocity - if (opts.turb_adve) pimpl->hskpng_turb_vel_calc(); + if (opts.turb_adve) + { + // calc tke (diss_rate now holds TKE, not dissipation rate!) + pimpl->hskpng_tke(); + // calc turbulent perturbation of velocity +// pimpl->hskpng_turb_vel(); + } // advection, it invalidates i,j,k and ijk! if (opts.adve) pimpl->adve(); // 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(); +// 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) From 33da0ff30340df9750508d92ead9ec7ea07155bb Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 27 Sep 2018 15:57:55 +0200 Subject: [PATCH 04/44] working? turb adve --- src/detail/urand.hpp | 39 ++++++++- src/impl/particles_impl.ipp | 1 + src/impl/particles_impl_hskpng_tke.ipp | 39 +++++++++ src/impl/particles_impl_hskpng_turb_vel.ipp | 87 +++++++++++++++++++++ src/impl/particles_impl_turb_adve.ipp | 36 +++++++++ src/particles.tpp | 2 + src/particles_step.ipp | 4 +- 7 files changed, 205 insertions(+), 3 deletions(-) create mode 100644 src/impl/particles_impl_hskpng_tke.ipp create mode 100644 src/impl/particles_impl_hskpng_turb_vel.ipp create mode 100644 src/impl/particles_impl_turb_adve.ipp diff --git a/src/detail/urand.hpp b/src/detail/urand.hpp index 7ca11b42a..42719cd8d 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_normal01_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_normal01_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 8bbd13569..28c4984ea 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -445,6 +445,7 @@ namespace libcloudphxx ); void adve(); + void turb_adve(); template void adve_calc(bool, thrust_size_t = 0); void sedi(); diff --git a/src/impl/particles_impl_hskpng_tke.ipp b/src/impl/particles_impl_hskpng_tke.ipp new file mode 100644 index 000000000..57812a7df --- /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(), common__turbulence__tke(L)); + } + }; +}; 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..59e46d604 --- /dev/null +++ b/src/impl/particles_impl_hskpng_turb_vel.ipp @@ -0,0 +1,87 @@ +// 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() + { + 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}; + std::vector*> vel_turbs_vctrs(&vel_turbs_vctrs_a[0], &vel_turbs_vctrs_a[0]+n_dims); + for(auto vctr_ptr : vel_turbs_vctrs) + { + 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( + vctr_ptr->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( + up.begin(), + thrust::make_permutation_iterator(tau.begin(), ijk.begin()), + thrust::make_permutation_iterator(tke.begin(), ijk.begin()), + r_normal.begin() + )) + n_part, + vctr_ptr->begin(), + detail::common__turbulence__update_turb_vel(opts_init.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/particles.tpp b/src/particles.tpp index 4ddd40cb8..508fd8b6c 100644 --- a/src/particles.tpp +++ b/src/particles.tpp @@ -69,6 +69,7 @@ #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_sort.ipp" #include "impl/particles_impl_hskpng_count.ipp" #include "impl/particles_impl_hskpng_remove.ipp" @@ -79,6 +80,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_step.ipp b/src/particles_step.ipp index b4d48918a..716da4b58 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -336,14 +336,14 @@ namespace libcloudphxx // calc tke (diss_rate now holds TKE, not dissipation rate!) pimpl->hskpng_tke(); // calc turbulent perturbation of velocity -// pimpl->hskpng_turb_vel(); + pimpl->hskpng_turb_vel(); } // advection, it invalidates i,j,k and ijk! if (opts.adve) pimpl->adve(); // 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(); + 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) From 3d6dad566765e2c21ad37a8bcbc9e3410c190b83 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 27 Sep 2018 16:05:43 +0200 Subject: [PATCH 05/44] track turbulence.hpp --- include/libcloudph++/common/turbulence.hpp | 99 ++++++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 include/libcloudph++/common/turbulence.hpp diff --git a/include/libcloudph++/common/turbulence.hpp b/include/libcloudph++/common/turbulence.hpp new file mode 100644 index 000000000..bb5ac6e46 --- /dev/null +++ b/include/libcloudph++/common/turbulence.hpp @@ -0,0 +1,99 @@ +#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.)) + ); + + typedef divide_typeof_helper< + si::dimensionless, + si::length + >::type one_over_length; + libcloudphxx_const(one_over_length, a_1, real_t(3e-4), 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; + + template + quantity length_scale( + quantity dx + ) + { return dx;} + + template + quantity length_scale( + quantity dx, + quantity dz + ) + { return quantity(sqrt(dx*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);} + + template + BOOST_GPU_ENABLED + quantity tke( + const quantity &diss_rate, // dissipation rate (epsilon) + const quantity &L // characteristic length-scale + ) + { + return pow(L * diss_rate / C_E(), real_t(2./3.)); + }; + + 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 &vt, + const quantity &tau, + const quantity &dt, + const quantity &tke, + const quantity &r_normal + ) + { + quantity exp_m_dt_ov_tau(exp(-dt / tau)); + return quantity(vt * 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); + }; + }; + }; +}; From a82f47745bde3ad6813d733662db015819cd89dd Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 27 Sep 2018 16:38:53 +0200 Subject: [PATCH 06/44] turb fix: a02 fixes --- src/detail/urand.hpp | 4 ++-- src/impl/particles_impl_cond_common.ipp | 2 ++ src/particles_multi_gpu_step.ipp | 6 ++++-- 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/detail/urand.hpp b/src/detail/urand.hpp index 42719cd8d..7fe52b05a 100644 --- a/src/detail/urand.hpp +++ b/src/detail/urand.hpp @@ -134,7 +134,7 @@ namespace libcloudphxx _unused(status); } - void generate_normal01_n( + void generate_normal_n( thrust_device::vector &v, const thrust_size_t n ) @@ -144,7 +144,7 @@ namespace libcloudphxx _unused(status); } - void generate_normal01_n( + void generate_normal_n( thrust_device::vector &v, const thrust_size_t n ) diff --git a/src/impl/particles_impl_cond_common.ipp b/src/impl/particles_impl_cond_common.ipp index e795a9e09..7ea1677e2 100644 --- a/src/impl/particles_impl_cond_common.ipp +++ b/src/impl/particles_impl_cond_common.ipp @@ -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;; 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 From 986014446594c14651dbc4803740c97690a61874 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 28 Sep 2018 11:55:57 +0200 Subject: [PATCH 07/44] python bindings for turb adve --- bindings/python/lgrngn.hpp | 4 ++++ bindings/python/lib.cpp | 3 +++ 2 files changed, 7 insertions(+) diff --git a/bindings/python/lgrngn.hpp b/bindings/python/lgrngn.hpp index 839ef44e4..729688f27 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 2336a5bd0..271e93b4f 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -235,6 +235,7 @@ BOOST_PYTHON_MODULE(libcloudphxx) // classes bp::class_>("opts_t") .def_readwrite("adve", &lgr::opts_t::adve) + .def_readwrite("turb_adve", &lgr::opts_t::turb_adve) .def_readwrite("sedi", &lgr::opts_t::sedi) .def_readwrite("cond", &lgr::opts_t::cond) .def_readwrite("coal", &lgr::opts_t::coal) @@ -269,6 +270,7 @@ 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_switch", &lgr::opts_init_t::turb_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) @@ -318,6 +320,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, ( From 0da6ca06e72bbb544852bbad67fd96f1c94cc458 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 28 Sep 2018 12:20:14 +0200 Subject: [PATCH 08/44] hskpng_tke fixes: --- include/libcloudph++/common/turbulence.hpp | 2 +- src/impl/particles_impl_hskpng_tke.ipp | 6 +++--- src/particles.tpp | 1 + 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/include/libcloudph++/common/turbulence.hpp b/include/libcloudph++/common/turbulence.hpp index bb5ac6e46..86386ac27 100644 --- a/include/libcloudph++/common/turbulence.hpp +++ b/include/libcloudph++/common/turbulence.hpp @@ -68,7 +68,7 @@ namespace libcloudphxx const quantity &L // characteristic length-scale ) { - return pow(L * diss_rate / C_E(), real_t(2./3.)); + return 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 diff --git a/src/impl/particles_impl_hskpng_tke.ipp b/src/impl/particles_impl_hskpng_tke.ipp index 57812a7df..c5ead9db1 100644 --- a/src/impl/particles_impl_hskpng_tke.ipp +++ b/src/impl/particles_impl_hskpng_tke.ipp @@ -27,13 +27,13 @@ namespace libcloudphxx 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(), common__turbulence__tke(L)); + thrust::transform(diss_rate.begin(), diss_rate.end(), diss_rate.begin(), detail::common__turbulence__tke(L)); } }; }; diff --git a/src/particles.tpp b/src/particles.tpp index 508fd8b6c..90ccf4ef4 100644 --- a/src/particles.tpp +++ b/src/particles.tpp @@ -70,6 +70,7 @@ #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_tke.ipp" #include "impl/particles_impl_hskpng_sort.ipp" #include "impl/particles_impl_hskpng_count.ipp" #include "impl/particles_impl_hskpng_remove.ipp" From 7fa7154fbcbabf31eb881efdacbdce4d1d4126a9 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 28 Sep 2018 12:35:48 +0200 Subject: [PATCH 09/44] fix bindings for turb --- bindings/python/lib.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index 271e93b4f..b7e32d024 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -311,6 +311,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, ( From 63d87f94cb5aa12ffcc841d0a143f7b81ec82b4e Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 28 Sep 2018 14:20:01 +0200 Subject: [PATCH 10/44] turb_vel calc: fis the loop --- src/impl/particles_impl_hskpng_turb_vel.ipp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/impl/particles_impl_hskpng_turb_vel.ipp b/src/impl/particles_impl_hskpng_turb_vel.ipp index 59e46d604..9861f54d4 100644 --- a/src/impl/particles_impl_hskpng_turb_vel.ipp +++ b/src/impl/particles_impl_hskpng_turb_vel.ipp @@ -59,26 +59,24 @@ namespace libcloudphxx 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}; - std::vector*> vel_turbs_vctrs(&vel_turbs_vctrs_a[0], &vel_turbs_vctrs_a[0]+n_dims); - for(auto vctr_ptr : vel_turbs_vctrs) + for(int i=0; ibegin(), + 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( - up.begin(), + 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, - vctr_ptr->begin(), + vel_turbs_vctrs_a[i]->begin(), detail::common__turbulence__update_turb_vel(opts_init.dt) ); } From 59a3fb875d009e5d6f54ce202791c00cdbb3e678 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 28 Sep 2018 14:27:44 +0200 Subject: [PATCH 11/44] test for turb adve --- tests/python/unit/lgrngn_turb_adve.py | 71 +++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) create mode 100644 tests/python/unit/lgrngn_turb_adve.py diff --git a/tests/python/unit/lgrngn_turb_adve.py b/tests/python/unit/lgrngn_turb_adve.py new file mode 100644 index 000000000..9b74a377b --- /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_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 From 93ad91d73b1f52e7565a5fe58061b4c30c4d6c3d Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 28 Sep 2018 15:29:12 +0200 Subject: [PATCH 12/44] construct quantities to be returned in turbulence.hpp --- include/libcloudph++/common/turbulence.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/include/libcloudph++/common/turbulence.hpp b/include/libcloudph++/common/turbulence.hpp index 86386ac27..27fa56f60 100644 --- a/include/libcloudph++/common/turbulence.hpp +++ b/include/libcloudph++/common/turbulence.hpp @@ -18,6 +18,11 @@ namespace libcloudphxx , real_t(1./3.)) ); +#if !defined(__NVCC__) + using std::pow; + using std::sqrt; +#endif + typedef divide_typeof_helper< si::dimensionless, si::length @@ -44,7 +49,7 @@ namespace libcloudphxx quantity length_scale( quantity dx ) - { return dx;} + { return quantity(dx);} template quantity length_scale( @@ -68,7 +73,7 @@ namespace libcloudphxx const quantity &L // characteristic length-scale ) { - return 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; + 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 From 953f062143faece53f0a2538686c08e68a477208 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 28 Sep 2018 15:55:03 +0200 Subject: [PATCH 13/44] api_lgrngn test: test turb api --- tests/python/unit/api_lgrngn.py | 55 +++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/tests/python/unit/api_lgrngn.py b/tests/python/unit/api_lgrngn.py index 5cca42ce6..8fbc230a3 100644 --- a/tests/python/unit/api_lgrngn.py +++ b/tests/python/unit/api_lgrngn.py @@ -366,7 +366,24 @@ 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_switch=True +opts.turb_adve=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_switch=False +opts.turb_adve=False # ---------- # 2D (periodic horizontal domain) @@ -416,6 +433,25 @@ 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_switch=True +opts.turb_adve=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_switch=False +opts.turb_adve=False # ---------- @@ -453,6 +489,25 @@ 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_switch=True +opts.turb_adve=False +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_switch=False +opts.turb_adve=False + # 3D with large_tail option print "3D large tail" opts_init.sd_conc_large_tail = 1 From 9c1d39eee2cb7a68822e029685d3a3583135145e Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 28 Sep 2018 15:55:26 +0200 Subject: [PATCH 14/44] prtlc_step: init_e2l of diss rate only if it is not null --- src/particles_step.ipp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 716da4b58..d950ccd9b 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -91,7 +91,7 @@ namespace libcloudphxx } if (pimpl->l2e[&pimpl->diss_rate].size() == 0) - pimpl->init_e2l(diss_rate, &pimpl->diss_rate); + if (!diss_rate.is_null()) pimpl->init_e2l(diss_rate, &pimpl->diss_rate); // did rhod change pimpl->var_rho = !rhod.is_null(); From 7363f00d932df039110ab2dfc4f66ddfe505906b Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 1 Oct 2018 12:16:19 +0200 Subject: [PATCH 15/44] ptrcls_step: throw error if 0D with turb_adve --- src/particles_step.ipp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/particles_step.ipp b/src/particles_step.ipp index d950ccd9b..f3a1e2545 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -295,6 +295,9 @@ namespace libcloudphxx if(opts.turb_adve && !pimpl->opts_init.turb_switch) throw std::runtime_error("all turbulence sgs was switched off in opts_init, 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 From 3aa4fcb49191bea79c9ddf775177421c65defa4a Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 1 Oct 2018 12:54:28 +0200 Subject: [PATCH 16/44] separate turb_adve/cond switches, add the ssp array --- include/libcloudph++/lgrngn/opts.hpp | 4 +- include/libcloudph++/lgrngn/opts_init.hpp | 6 ++- src/impl/particles_impl.ipp | 1 + src/impl/particles_impl_hskpng_remove.ipp | 7 ++- src/impl/particles_impl_hskpng_resize.ipp | 7 ++- src/impl/particles_impl_init_hskpng_npart.ipp | 8 +++- src/impl/particles_impl_init_sync.ipp | 2 +- src/impl/particles_impl_rcyc.ipp | 7 ++- ...les_multi_gpu_impl_step_async_and_copy.ipp | 8 +++- src/particles_step.ipp | 44 +++++++++++++++---- 10 files changed, 75 insertions(+), 19 deletions(-) diff --git a/include/libcloudph++/lgrngn/opts.hpp b/include/libcloudph++/lgrngn/opts.hpp index fd26b012e..b0ede176e 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, turb_adve; + bool adve, sedi, cond, coal, src, rcyc, turb_adve, turb_cond; // RH limit for drop growth real_t RH_max; @@ -32,7 +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_adve(false), turb_cond(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 6132119dc..52be15b76 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -104,7 +104,8 @@ namespace libcloudphxx 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 - turb_switch; // if true, turbulent motion of SDs is modeled + turb_adve_switch, // if true, turbulent motion of SDs is modeled + turb_cond_switch; // if true, turbulent condensation of SDs is modeled int sstp_chem; real_t chem_rho; @@ -143,7 +144,8 @@ namespace libcloudphxx coal_switch(true), // coalescence turned on by default src_switch(false), // source turned off by default exact_sstp_cond(false), - turb_switch(false), + turb_cond_switch(false), + turb_adve_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), diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 28c4984ea..ee3884de9 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -81,6 +81,7 @@ namespace libcloudphxx up, // turbulent perturbation of velocity vp, // turbulent perturbation of velocity wp, // turbulent perturbation of velocity + ssp, // 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 diff --git a/src/impl/particles_impl_hskpng_remove.ipp b/src/impl/particles_impl_hskpng_remove.ipp index f4292816c..15e7b52d4 100644 --- a/src/impl/particles_impl_hskpng_remove.ipp +++ b/src/impl/particles_impl_hskpng_remove.ipp @@ -32,12 +32,17 @@ namespace libcloudphxx 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.turb_switch) + if(opts_init.turb_adve_switch) { if (opts_init.nx != 0) real_t_vctrs.push_back(&up); if (opts_init.ny != 0) real_t_vctrs.push_back(&vp); if (opts_init.nz != 0) real_t_vctrs.push_back(&wp); } + else if(opts_init.turb_cond_switch) + if (opts_init.nz != 0) real_t_vctrs.push_back(&wp); + + if(opts_init.turb_cond_switch) + if (opts_init.nz != 0) real_t_vctrs.push_back(&ssp); if(opts_init.sstp_cond>1 && opts_init.exact_sstp_cond) { diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index 6c09d7bde..4980a1d8c 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -33,12 +33,17 @@ namespace libcloudphxx if (opts_init.ny != 0) y.resize(n_part); if (opts_init.nz != 0) z.resize(n_part); - if(opts_init.turb_switch) // TODO: they are not needed by turb_coal, but diss_rate is - make more switches + 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); } + else if(opts_init.turb_cond_switch) + if (opts_init.nz != 0) wp.resize(n_part); + + if(opts_init.turb_cond_switch) + if (opts_init.nz != 0) ssp.resize(n_part); if(opts_init.chem_switch || opts_init.sstp_cond > 1 || n_dims >= 2) { diff --git a/src/impl/particles_impl_init_hskpng_npart.ipp b/src/impl/particles_impl_init_hskpng_npart.ipp index f76b37937..e6e4e191d 100644 --- a/src/impl/particles_impl_init_hskpng_npart.ipp +++ b/src/impl/particles_impl_init_hskpng_npart.ipp @@ -23,12 +23,18 @@ 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); - if(opts_init.turb_switch) + // 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.); } + else if(opts_init.turb_cond_switch) + if (opts_init.nz != 0) wp.resize(opts_init.n_sd_max, 0.); + + if(opts_init.turb_cond_switch) + if (opts_init.nz != 0) 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 diff --git a/src/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index a19c6d6fc..2c4743b51 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/src/impl/particles_impl_init_sync.ipp @@ -20,7 +20,7 @@ namespace libcloudphxx 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_switch) + if(opts_init.turb_cond_switch || opts_init.turb_adve_switch) diss_rate.resize(n_cell); // memory allocation for vector fields (Arakawa-C grid) diff --git a/src/impl/particles_impl_rcyc.ipp b/src/impl/particles_impl_rcyc.ipp index 6ec24cb17..812f66412 100644 --- a/src/impl/particles_impl_rcyc.ipp +++ b/src/impl/particles_impl_rcyc.ipp @@ -113,12 +113,17 @@ namespace libcloudphxx detail::copy_prop(chem_bgn[i], sorted_id, n_flagged); } - if(opts_init.turb_switch) + 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); } + else if(opts_init.turb_cond_switch) + if (opts_init.nz > 0) detail::copy_prop(wp.begin(), sorted_id, n_flagged); + + if(opts_init.turb_cond_switch) + if (opts_init.nz > 0) detail::copy_prop(ssp.begin(), sorted_id, n_flagged); { namespace arg = thrust::placeholders; 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 2d19a0861..1332db391 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 @@ -138,12 +138,18 @@ namespace libcloudphxx if(particles[dev_id]->pimpl->const_p) real_t_vctrs.push_back(&sstp_tmp_p); } - if(glob_opts_init.turb_switch) + 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_step.ipp b/src/particles_step.ipp index f3a1e2545..48756fbc2 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -71,11 +71,17 @@ 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_switch && diss_rate.is_null()) - throw std::runtime_error("turbulence was not switched off and diss_rate is empty"); + if (pimpl->opts_init.turb_adve_switch && diss_rate.is_null()) + throw std::runtime_error("turbulent advection was not switched off and diss_rate is empty"); - if (!pimpl->opts_init.turb_switch && !diss_rate.is_null()) - throw std::runtime_error("turbulence was switched off and diss_rate is not empty"); + if (!pimpl->opts_init.turb_adve_switch && !diss_rate.is_null()) + throw std::runtime_error("turbulent advection was switched off and diss_rate is not empty"); + + if (pimpl->opts_init.turb_cond_switch && diss_rate.is_null()) + throw std::runtime_error("turbulent condensation was not switched off and diss_rate is empty"); + + if (!pimpl->opts_init.turb_cond_switch && !diss_rate.is_null()) + throw std::runtime_error("turbulent condensation was switched off and diss_rate is not empty"); // if (pimpl->l2e[&pimpl->courant_x].size() == 0) // TODO: y, z,... @@ -112,12 +118,12 @@ namespace libcloudphxx 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"); + 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_switch) - assert(*thrust::min_element(pimpl->rhod.begin(), pimpl->rhod.end()) >= 0); + 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) assert(pimpl->opts_init.adve_scheme != as_t::pred_corr || (courant_x.is_null() || ((*(thrust::min_element(pimpl->courant_x.begin(), pimpl->courant_x.end()))) >= real_t(-2.) )) ); @@ -292,12 +298,18 @@ 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_switch) - throw std::runtime_error("all turbulence sgs was switched off in opts_init, but turb_adve==True"); + 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_cond && !pimpl->opts_init.turb_cond_switch) + throw std::runtime_error("turb_cond_swtich=False, but turb_cond==True"); if(opts.turb_adve && pimpl->n_dims==0) throw std::runtime_error("turbulent advection does not work in 0D"); + if(opts.turb_cond && pimpl->n_dims<2) + throw std::runtime_error("turbulent condensation works only in 2D and 3D"); + if (opts.chem_dsl) { // saving rv to be used as rv_old @@ -334,13 +346,27 @@ namespace libcloudphxx } } - if (opts.turb_adve) + if (opts.turb_adve || opts.turb_cond) { // calc tke (diss_rate now holds TKE, not dissipation 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_vert_vel(); + } + + if(opts.turb_cond) + { + // calculate the turbulent supersaturation; used in the next step during condensation in step_cond - is it a problem? + //pimpl->hskpng_turb_ss(); + } // advection, it invalidates i,j,k and ijk! if (opts.adve) pimpl->adve(); From 572e85b10ba09dbad32b701182cce8cb24500a9f Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 1 Oct 2018 14:57:59 +0200 Subject: [PATCH 17/44] working calc of ssp --- include/libcloudph++/common/turbulence.hpp | 36 +++++++++++++++++++-- src/impl/particles_impl.ipp | 3 +- src/impl/particles_impl_hskpng_turb_vel.ipp | 4 +-- src/particles.tpp | 1 + src/particles_step.ipp | 4 +-- 5 files changed, 40 insertions(+), 8 deletions(-) diff --git a/include/libcloudph++/common/turbulence.hpp b/include/libcloudph++/common/turbulence.hpp index 27fa56f60..f7bd26118 100644 --- a/include/libcloudph++/common/turbulence.hpp +++ b/include/libcloudph++/common/turbulence.hpp @@ -27,7 +27,7 @@ namespace libcloudphxx si::dimensionless, si::length >::type one_over_length; - libcloudphxx_const(one_over_length, a_1, real_t(3e-4), 1. / si::meters); + libcloudphxx_const(one_over_length, a_1, real_t(3e-4), real_t(1) / si::meters); typedef divide_typeof_helper< si::area, @@ -45,6 +45,16 @@ namespace libcloudphxx 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 @@ -89,7 +99,7 @@ namespace libcloudphxx template BOOST_GPU_ENABLED quantity update_turb_vel( - const quantity &vt, + const quantity &wp, const quantity &tau, const quantity &dt, const quantity &tke, @@ -97,7 +107,27 @@ namespace libcloudphxx ) { quantity exp_m_dt_ov_tau(exp(-dt / tau)); - return quantity(vt * 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); + 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/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index ee3884de9..6d8041888 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -408,7 +408,8 @@ namespace libcloudphxx void hskpng_vterm_all(); void hskpng_vterm_invalid(); void hskpng_tke(); - void hskpng_turb_vel(); + void hskpng_turb_vel(const bool only_vertical = false); + void hskpng_turb_ss(); void hskpng_remove_n0(); void hskpng_resize_npart(); diff --git a/src/impl/particles_impl_hskpng_turb_vel.ipp b/src/impl/particles_impl_hskpng_turb_vel.ipp index 9861f54d4..84ad04483 100644 --- a/src/impl/particles_impl_hskpng_turb_vel.ipp +++ b/src/impl/particles_impl_hskpng_turb_vel.ipp @@ -52,7 +52,7 @@ namespace libcloudphxx }; // calc the SGS turbulent velocity component template - void particles_t::impl::hskpng_turb_vel() + 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 @@ -60,7 +60,7 @@ namespace libcloudphxx thrust_device::vector &r_normal(tmp_device_real_part); thrust_device::vector * vel_turbs_vctrs_a[] = {&up, &wp, &vp}; - for(int i=0; ihskpng_turb_vert_vel(); + pimpl->hskpng_turb_vel(true); } if(opts.turb_cond) { // calculate the turbulent supersaturation; used in the next step during condensation in step_cond - is it a problem? - //pimpl->hskpng_turb_ss(); + pimpl->hskpng_turb_ss(); } // advection, it invalidates i,j,k and ijk! From e63d79e3007987c63f054b06cefcebf37320eed4 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 1 Oct 2018 15:16:45 +0200 Subject: [PATCH 18/44] track hskpng_turb_ss --- src/impl/particles_impl_hskpng_turb_ss.ipp | 73 ++++++++++++++++++++++ 1 file changed, 73 insertions(+) create mode 100644 src/impl/particles_impl_hskpng_turb_ss.ipp 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..becba313f --- /dev/null +++ b/src/impl/particles_impl_hskpng_turb_ss.ipp @@ -0,0 +1,73 @@ +// 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; + } + }; + + struct common__turbulence__turb_dot_ss + { + template + BOOST_GPU_ENABLED + real_t operator()(tpl_t tpl) + { + 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() + { + // calc tau_relax + moms_calc(rw2.begin(), real_t(1./2), false); + thrust_device::vector &tau_rlx(count_mom); + thrust::transform( + tau_rlx.begin(), + dv.begin(), + tau_rlx.end(), + tau_rlx.begin(), + detail::common__turbulence__tau_relax() + ); + + // 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() + ); + } + }; +}; From f5dc095680db13e79ce9817138f7bbacf7437175 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 1 Oct 2018 15:17:07 +0200 Subject: [PATCH 19/44] add the dot_ssp array --- src/impl/particles_impl.ipp | 1 + src/impl/particles_impl_hskpng_remove.ipp | 8 +++++++- src/impl/particles_impl_hskpng_resize.ipp | 8 +++++++- src/impl/particles_impl_init_hskpng_npart.ipp | 8 +++++++- src/impl/particles_impl_rcyc.ipp | 8 +++++++- src/particles_step.ipp | 4 ++-- 6 files changed, 31 insertions(+), 6 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 6d8041888..3a2fb7fb6 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -82,6 +82,7 @@ namespace libcloudphxx 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 diff --git a/src/impl/particles_impl_hskpng_remove.ipp b/src/impl/particles_impl_hskpng_remove.ipp index 15e7b52d4..08e20c25e 100644 --- a/src/impl/particles_impl_hskpng_remove.ipp +++ b/src/impl/particles_impl_hskpng_remove.ipp @@ -42,7 +42,13 @@ namespace libcloudphxx if (opts_init.nz != 0) real_t_vctrs.push_back(&wp); if(opts_init.turb_cond_switch) - if (opts_init.nz != 0) real_t_vctrs.push_back(&ssp); + { + if (opts_init.nz != 0) + { + real_t_vctrs.push_back(&ssp); + real_t_vctrs.push_back(&dot_ssp); + } + } if(opts_init.sstp_cond>1 && opts_init.exact_sstp_cond) { diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index 4980a1d8c..2ba2f676a 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -43,7 +43,13 @@ namespace libcloudphxx if (opts_init.nz != 0) wp.resize(n_part); if(opts_init.turb_cond_switch) - if (opts_init.nz != 0) ssp.resize(n_part); + { + if (opts_init.nz != 0) + { + ssp.resize(n_part); + dot_ssp.resize(n_part); + } + } if(opts_init.chem_switch || opts_init.sstp_cond > 1 || n_dims >= 2) { diff --git a/src/impl/particles_impl_init_hskpng_npart.ipp b/src/impl/particles_impl_init_hskpng_npart.ipp index e6e4e191d..454f16d6a 100644 --- a/src/impl/particles_impl_init_hskpng_npart.ipp +++ b/src/impl/particles_impl_init_hskpng_npart.ipp @@ -34,7 +34,13 @@ namespace libcloudphxx if (opts_init.nz != 0) wp.resize(opts_init.n_sd_max, 0.); if(opts_init.turb_cond_switch) - if (opts_init.nz != 0) ssp.resize(opts_init.n_sd_max, 0.); + { + if (opts_init.nz != 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 diff --git a/src/impl/particles_impl_rcyc.ipp b/src/impl/particles_impl_rcyc.ipp index 812f66412..f340d72ae 100644 --- a/src/impl/particles_impl_rcyc.ipp +++ b/src/impl/particles_impl_rcyc.ipp @@ -123,7 +123,13 @@ namespace libcloudphxx if (opts_init.nz > 0) detail::copy_prop(wp.begin(), sorted_id, n_flagged); if(opts_init.turb_cond_switch) - if (opts_init.nz > 0) detail::copy_prop(ssp.begin(), sorted_id, n_flagged); + { + if (opts_init.nz > 0) + { + 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/particles_step.ipp b/src/particles_step.ipp index bd089e3ea..40a2929a3 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -364,8 +364,8 @@ namespace libcloudphxx if(opts.turb_cond) { - // calculate the turbulent supersaturation; used in the next step during condensation in step_cond - is it a problem? - pimpl->hskpng_turb_ss(); + // 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! From e2eddc5db30a686974db7d51762997d6d1b3a4dd Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 1 Oct 2018 16:09:01 +0200 Subject: [PATCH 20/44] per-cell turb_cond --- src/impl/particles_impl.ipp | 1 + src/impl/particles_impl_cond.ipp | 60 +++++++++++++++------- src/impl/particles_impl_hskpng_turb_ss.ipp | 2 +- src/impl/particles_impl_sstp.ipp | 16 ++++++ src/particles_step.ipp | 2 + 5 files changed, 61 insertions(+), 20 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 3a2fb7fb6..f93ef56e8 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -476,6 +476,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_cond.ipp b/src/impl/particles_impl_cond.ipp index 2a83e6569..cc9d7ae62 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -12,7 +12,8 @@ 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 ) { // --- calc liquid water content before cond --- @@ -34,24 +35,45 @@ namespace libcloudphxx ); // 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::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_transform_iterator(thrust::make_permutation_iterator(RH.begin(), ijk.begin()), ssp.begin(), thrust::add()), + thrust::make_permutation_iterator(eta.begin(), ijk.begin()), + rd3.begin(), + kpa.begin(), + vt.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( + 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) + ); nancheck(rw2, "rw2 after condensation (no sub-steps"); // calculating the 3rd wet moment after condensation diff --git a/src/impl/particles_impl_hskpng_turb_ss.ipp b/src/impl/particles_impl_hskpng_turb_ss.ipp index becba313f..61d0b8111 100644 --- a/src/impl/particles_impl_hskpng_turb_ss.ipp +++ b/src/impl/particles_impl_hskpng_turb_ss.ipp @@ -47,8 +47,8 @@ namespace libcloudphxx thrust_device::vector &tau_rlx(count_mom); thrust::transform( tau_rlx.begin(), - dv.begin(), tau_rlx.end(), + dv.begin(), tau_rlx.begin(), detail::common__turbulence__tau_relax() ); 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/particles_step.ipp b/src/particles_step.ipp index 40a2929a3..2fb3f23b9 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -176,6 +176,8 @@ 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); } From 2cd8b75ed441d69f0ac653a000aac5e44281fe0b Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 1 Oct 2018 17:44:04 +0200 Subject: [PATCH 21/44] per-prtcl turb_cond --- bindings/python/lib.cpp | 6 +- src/impl/particles_impl.ipp | 6 +- src/impl/particles_impl_cond.ipp | 16 ++++- src/impl/particles_impl_cond_sstp.ipp | 62 ++++++++++++++----- src/impl/particles_impl_hskpng_resize.ipp | 2 +- src/impl/particles_impl_hskpng_turb_ss.ipp | 3 +- src/impl/particles_impl_init_hskpng_npart.ipp | 2 +- src/particles_step.ipp | 6 +- tests/python/unit/api_lgrngn.py | 12 ++-- tests/python/unit/lgrngn_turb_adve.py | 2 +- 10 files changed, 84 insertions(+), 33 deletions(-) diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index b7e32d024..36eb74c12 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -235,7 +235,6 @@ BOOST_PYTHON_MODULE(libcloudphxx) // classes bp::class_>("opts_t") .def_readwrite("adve", &lgr::opts_t::adve) - .def_readwrite("turb_adve", &lgr::opts_t::turb_adve) .def_readwrite("sedi", &lgr::opts_t::sedi) .def_readwrite("cond", &lgr::opts_t::cond) .def_readwrite("coal", &lgr::opts_t::coal) @@ -245,6 +244,8 @@ 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) ; bp::class_>("opts_init_t") .add_property("dry_distros", &lgrngn::get_dd, &lgrngn::set_dd) @@ -270,7 +271,8 @@ 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_switch", &lgr::opts_init_t::turb_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("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) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index f93ef56e8..7f9f86e6d 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -410,7 +410,7 @@ namespace libcloudphxx void hskpng_vterm_invalid(); void hskpng_tke(); void hskpng_turb_vel(const bool only_vertical = false); - void hskpng_turb_ss(); + void hskpng_turb_dot_ss(); void hskpng_remove_n0(); void hskpng_resize_npart(); @@ -454,8 +454,8 @@ namespace libcloudphxx 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); 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 &); diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index cc9d7ae62..744448f87 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -13,9 +13,11 @@ namespace libcloudphxx void particles_t::impl::cond( const real_t &dt, const real_t &RH_max, - const bool &turb_cond + const bool turb_cond ) { + namespace arg = thrust::placeholders; + // --- calc liquid water content before cond --- hskpng_sort(); thrust_device::vector &drv(tmp_device_real_cell); @@ -37,6 +39,15 @@ namespace libcloudphxx // calculating drop growth in a timestep using backward Euler // 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 @@ -45,7 +56,7 @@ namespace libcloudphxx 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_transform_iterator(thrust::make_permutation_iterator(RH.begin(), ijk.begin()), ssp.begin(), thrust::add()), + RH_plus_ssp.begin(), thrust::make_permutation_iterator(eta.begin(), ijk.begin()), rd3.begin(), kpa.begin(), @@ -55,6 +66,7 @@ namespace libcloudphxx 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() diff --git a/src/impl/particles_impl_cond_sstp.ipp b/src/impl/particles_impl_cond_sstp.ipp index 7227b70be..0089c4530 100644 --- a/src/impl/particles_impl_cond_sstp.ipp +++ b/src/impl/particles_impl_cond_sstp.ipp @@ -9,11 +9,36 @@ namespace libcloudphxx { namespace lgrngn { + namespace detail + { + template + struct RH_hlpr : thrust::unary_function, real_t>&, real_t> + { + const bool turb_cond; + RH resolved_RH; + + RH_hlpr(RH_formula_t::RH_formula_t RH_formula, const bool turb_cond): + resolved_RH(RH_formula), + turb_cond(turb_cond) + {} + + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple, real_t> &tpl) + { + real_t res = resolved_RH(thrust::get<0>(tpl)); + return turb_cond ? res + thrust::get<1>(tpl) : res; + } + }; + }; + 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 @@ -57,7 +82,6 @@ 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 @@ -82,16 +106,21 @@ namespace libcloudphxx Tp.begin(), // particle-specific p pressure_iter, - // particle-specific RH + // particle-specific RH, resolved + SGS + // TODO: if turb_cond_switch=0, ssp has size=0 - it shouldnt be a problem, but maybe? thrust::make_transform_iterator( thrust::make_zip_iterator( thrust::make_tuple( - pressure_iter, - sstp_tmp_rv.begin(), - Tp.begin() + thrust::make_zip_iterator( + thrust::make_tuple( + pressure_iter, + sstp_tmp_rv.begin(), + Tp.begin() + )), + ssp.begin() )), - detail::RH(opts_init.RH_formula) - ), + detail::RH_hlpr(opts_init.RH_formula, turb_cond) + ), // particle-specific eta thrust::make_transform_iterator( Tp.begin(), @@ -117,16 +146,21 @@ namespace libcloudphxx Tp.begin(), // particle-specific p sstp_tmp_p.begin(), - // particle-specific RH + // particle-specific RH, resolved + SGS + // TODO: if turb_cond_switch=0, ssp has size=0 - it shouldnt be a problem, but maybe? thrust::make_transform_iterator( thrust::make_zip_iterator( thrust::make_tuple( - sstp_tmp_p.begin(), - sstp_tmp_rv.begin(), - Tp.begin() + thrust::make_zip_iterator( + thrust::make_tuple( + sstp_tmp_p.begin(), + sstp_tmp_rv.begin(), + Tp.begin() + )), + ssp.begin() )), - detail::RH(opts_init.RH_formula) - ), + detail::RH_hlpr(opts_init.RH_formula, turb_cond) + ), // particle-specific eta thrust::make_transform_iterator( Tp.begin(), diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index 2ba2f676a..37e71d1a6 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -55,7 +55,7 @@ namespace libcloudphxx { 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_turb_ss.ipp b/src/impl/particles_impl_hskpng_turb_ss.ipp index 61d0b8111..3202b5f4c 100644 --- a/src/impl/particles_impl_hskpng_turb_ss.ipp +++ b/src/impl/particles_impl_hskpng_turb_ss.ipp @@ -23,6 +23,7 @@ namespace libcloudphxx } }; + template struct common__turbulence__turb_dot_ss { template @@ -66,7 +67,7 @@ namespace libcloudphxx thrust::make_permutation_iterator(tau_rlx.begin(), ijk.begin()) )) + n_part, dot_ssp.begin(), - detail::common__turbulence__turb_dot_ss() + detail::common__turbulence__turb_dot_ss() ); } }; diff --git a/src/impl/particles_impl_init_hskpng_npart.ipp b/src/impl/particles_impl_init_hskpng_npart.ipp index 454f16d6a..5c25afab4 100644 --- a/src/impl/particles_impl_init_hskpng_npart.ipp +++ b/src/impl/particles_impl_init_hskpng_npart.ipp @@ -60,7 +60,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); } diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 2fb3f23b9..63009f228 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -164,7 +164,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); @@ -179,7 +181,7 @@ namespace libcloudphxx 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); } } diff --git a/tests/python/unit/api_lgrngn.py b/tests/python/unit/api_lgrngn.py index 8fbc230a3..cdccdd6f7 100644 --- a/tests/python/unit/api_lgrngn.py +++ b/tests/python/unit/api_lgrngn.py @@ -368,7 +368,7 @@ def lognormal(lnr): print "1D turb" eps = arr_t([ 1e-4, 1e-4, 1e-4]) -opts_init.turb_switch=True +opts_init.turb_adve_switch=True opts.turb_adve=True prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) @@ -382,7 +382,7 @@ def lognormal(lnr): assert (frombuffer(prtcls.outbuf()) > 0).all() assert sum(frombuffer(prtcls.outbuf())) == opts_init.nx * opts_init.sd_conc -opts_init.turb_switch=False +opts_init.turb_adve_switch=False opts.turb_adve=False # ---------- @@ -435,7 +435,7 @@ def lognormal(lnr): print "2D turb" eps = arr_t([[ 1e-4, 1e-4], [ 1e-4, 1e-4]]) -opts_init.turb_switch=True +opts_init.turb_adve_switch=True opts.turb_adve=True prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) @@ -450,7 +450,7 @@ def lognormal(lnr): 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_switch=False +opts_init.turb_adve_switch=False opts.turb_adve=False @@ -491,7 +491,7 @@ def lognormal(lnr): print "3D turb" eps = arr_t([eps, eps]) -opts_init.turb_switch=True +opts_init.turb_adve_switch=True opts.turb_adve=False prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) @@ -505,7 +505,7 @@ 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 -opts_init.turb_switch=False +opts_init.turb_adve_switch=False opts.turb_adve=False # 3D with large_tail option diff --git a/tests/python/unit/lgrngn_turb_adve.py b/tests/python/unit/lgrngn_turb_adve.py index 9b74a377b..975abc869 100644 --- a/tests/python/unit/lgrngn_turb_adve.py +++ b/tests/python/unit/lgrngn_turb_adve.py @@ -20,7 +20,7 @@ def lognormal(lnr): Opts_init.dry_distros = {kappa:lognormal} Opts_init.coal_switch = False Opts_init.sedi_switch = False -Opts_init.turb_switch = True +Opts_init.turb_adve_switch = True Opts_init.dt = 1 From 6c8c7f050d27c033def21106655a689f72828648 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 9 Oct 2018 16:48:15 +0200 Subject: [PATCH 22/44] fix per-prtcl cond without turbulence --- src/impl/particles_impl_cond.ipp | 30 ++-- src/impl/particles_impl_cond_common.ipp | 43 ++--- src/impl/particles_impl_cond_sstp.ipp | 182 +++++++++++--------- src/impl/particles_impl_hskpng_turb_ss.ipp | 1 + src/impl/particles_impl_hskpng_turb_vel.ipp | 2 +- src/particles_step.ipp | 32 ++-- tests/python/unit/sstp_cond.py | 48 +++--- 7 files changed, 179 insertions(+), 159 deletions(-) diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index 744448f87..99b322e74 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -36,6 +36,16 @@ 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 // TODO: both calls almost identical, use std::bind or sth? if(turb_cond) @@ -52,15 +62,9 @@ namespace libcloudphxx 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()), + hlpr_zip_iter, thrust::make_permutation_iterator(p.begin(), ijk.begin()), - RH_plus_ssp.begin(), - thrust::make_permutation_iterator(eta.begin(), ijk.begin()), - rd3.begin(), - kpa.begin(), - vt.begin() + RH_plus_ssp.begin() ) ), rw2.begin(), // output @@ -72,15 +76,9 @@ namespace libcloudphxx 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()), + hlpr_zip_iter, 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() + thrust::make_permutation_iterator(RH.begin(), ijk.begin()) ) ), rw2.begin(), // output diff --git a/src/impl/particles_impl_cond_common.ipp b/src/impl/particles_impl_cond_common.ipp index 7ea1677e2..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) {} @@ -154,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; @@ -163,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; @@ -173,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 0089c4530..64dc562a5 100644 --- a/src/impl/particles_impl_cond_sstp.ipp +++ b/src/impl/particles_impl_cond_sstp.ipp @@ -12,21 +12,18 @@ namespace libcloudphxx namespace detail { template - struct RH_hlpr : thrust::unary_function, real_t>&, real_t> + struct RH_sgs : thrust::unary_function&, real_t> { - const bool turb_cond; RH resolved_RH; - RH_hlpr(RH_formula_t::RH_formula_t RH_formula, const bool turb_cond): - resolved_RH(RH_formula), - turb_cond(turb_cond) + RH_sgs(RH_formula_t::RH_formula_t RH_formula): + resolved_RH(RH_formula) {} BOOST_GPU_ENABLED - real_t operator()(const thrust::tuple, real_t> &tpl) + real_t operator()(const thrust::tuple &tpl) { - real_t res = resolved_RH(thrust::get<0>(tpl)); - return turb_cond ? res + thrust::get<1>(tpl) : res; + return resolved_RH(thrust::make_tuple(thrust::get<0>(tpl), thrust::get<1>(tpl), thrust::get<2>(tpl))) + thrust::get<3>(tpl); } }; }; @@ -36,7 +33,7 @@ namespace libcloudphxx const real_t &dt, const real_t &RH_max, const bool turb_cond - ) { + ) { namespace arg = thrust::placeholders; // prerequisite @@ -78,102 +75,123 @@ namespace libcloudphxx ); } + 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 - // TODO: these two function calls only differ in the way pressure is calculated, - // find a way to reuse it (e.g. std::bind?) + // TODO: these function calls only slightly differ... 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(), + if(turb_cond) + 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 pressure_iter, // particle-specific RH, resolved + SGS - // TODO: if turb_cond_switch=0, ssp has size=0 - it shouldnt be a problem, but maybe? thrust::make_transform_iterator( - thrust::make_zip_iterator( - thrust::make_tuple( - thrust::make_zip_iterator( - thrust::make_tuple( - pressure_iter, - sstp_tmp_rv.begin(), - Tp.begin() - )), - ssp.begin() + thrust::make_zip_iterator(thrust::make_tuple( + pressure_iter, + sstp_tmp_rv.begin(), + Tp.begin(), + ssp.begin() )), - detail::RH_hlpr(opts_init.RH_formula, turb_cond) - ), - // particle-specific eta + detail::RH_sgs(opts_init.RH_formula) + ) + )), + 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(thrust::make_tuple( // input - 2nd arg + hlpr_zip_iter, + // particle-specific p + pressure_iter, + // particle-specific RH, resolved + SGS + // TODO: if turb_cond_switch=0, ssp has size=0 - it shouldnt be a problem, but maybe? thrust::make_transform_iterator( - Tp.begin(), - detail::common__vterm__visc() - ), - rd3.begin(), - kpa.begin(), - vt.begin() - ) - ), - rw2.begin(), // output - detail::advance_rw2(dt, RH_max) - ); + thrust::make_zip_iterator(thrust::make_tuple( + pressure_iter, + sstp_tmp_rv.begin(), + Tp.begin() + )), + detail::RH(opts_init.RH_formula) + ) + )), + 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( - sstp_tmp_rh.begin(), - sstp_tmp_rv.begin(), - Tp.begin(), + if(turb_cond) + 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 sstp_tmp_p.begin(), // particle-specific RH, resolved + SGS // TODO: if turb_cond_switch=0, ssp has size=0 - it shouldnt be a problem, but maybe? thrust::make_transform_iterator( - thrust::make_zip_iterator( - thrust::make_tuple( - thrust::make_zip_iterator( - thrust::make_tuple( - sstp_tmp_p.begin(), - sstp_tmp_rv.begin(), - Tp.begin() - )), - ssp.begin() + thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_p.begin(), + sstp_tmp_rv.begin(), + Tp.begin(), + ssp.begin() )), - detail::RH_hlpr(opts_init.RH_formula, turb_cond) - ), - // particle-specific eta + detail::RH_sgs(opts_init.RH_formula) + ) + )), + 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(thrust::make_tuple( // input - 2nd arg + hlpr_zip_iter, + // particle-specific p + sstp_tmp_p.begin(), + // particle-specific RH, resolved + SGS + // TODO: if turb_cond_switch=0, ssp has size=0 - it shouldnt be a problem, but maybe? thrust::make_transform_iterator( - Tp.begin(), - detail::common__vterm__visc() - ), - rd3.begin(), - kpa.begin(), - vt.begin() - ) - ), - rw2.begin(), // output - detail::advance_rw2(dt, RH_max) - ); + thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_p.begin(), + sstp_tmp_rv.begin(), + Tp.begin() + )), + detail::RH(opts_init.RH_formula) + ) + )), + rw2.begin(), // output + detail::advance_rw2(dt, RH_max) + ); } // calc rw3_new - rw3_old diff --git a/src/impl/particles_impl_hskpng_turb_ss.ipp b/src/impl/particles_impl_hskpng_turb_ss.ipp index 3202b5f4c..67429b774 100644 --- a/src/impl/particles_impl_hskpng_turb_ss.ipp +++ b/src/impl/particles_impl_hskpng_turb_ss.ipp @@ -44,6 +44,7 @@ namespace libcloudphxx void particles_t::impl::hskpng_turb_dot_ss() { // calc tau_relax + moms_all(); moms_calc(rw2.begin(), real_t(1./2), false); thrust_device::vector &tau_rlx(count_mom); thrust::transform( diff --git a/src/impl/particles_impl_hskpng_turb_vel.ipp b/src/impl/particles_impl_hskpng_turb_vel.ipp index 84ad04483..8ec33ec47 100644 --- a/src/impl/particles_impl_hskpng_turb_vel.ipp +++ b/src/impl/particles_impl_hskpng_turb_vel.ipp @@ -60,7 +60,7 @@ namespace libcloudphxx 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) + 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( diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 63009f228..41c8a515d 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -71,17 +71,11 @@ 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 && diss_rate.is_null()) - throw std::runtime_error("turbulent advection was not switched off and diss_rate is empty"); + if ( (pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_switch) && diss_rate.is_null()) + throw std::runtime_error("turbulent advection and condesation are not switched off and diss_rate is empty"); - if (!pimpl->opts_init.turb_adve_switch && !diss_rate.is_null()) - throw std::runtime_error("turbulent advection was switched off and diss_rate is not empty"); - - if (pimpl->opts_init.turb_cond_switch && diss_rate.is_null()) - throw std::runtime_error("turbulent condensation was not switched off and diss_rate is empty"); - - if (!pimpl->opts_init.turb_cond_switch && !diss_rate.is_null()) - throw std::runtime_error("turbulent condensation was switched off and diss_rate is not empty"); + if ( !(pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_switch) && !diss_rate.is_null()) + throw std::runtime_error("turbulent advection and condesation are switched off and diss_rate is not empty"); // if (pimpl->l2e[&pimpl->courant_x].size() == 0) // TODO: y, z,... @@ -118,12 +112,14 @@ namespace libcloudphxx 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"); - nancheck(pimpl->diss_rate, " diss_rate after sync-in"); + if(pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_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); - assert(*thrust::min_element(pimpl->diss_rate.begin(), pimpl->diss_rate.end()) >= 0); + if(pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_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) assert(pimpl->opts_init.adve_scheme != as_t::pred_corr || (courant_x.is_null() || ((*(thrust::min_element(pimpl->courant_x.begin(), pimpl->courant_x.end()))) >= real_t(-2.) )) ); @@ -153,6 +149,12 @@ 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"); + + if(opts.turb_cond && pimpl->n_dims<2) + throw std::runtime_error("turbulent condensation works only in 2D and 3D"); + pimpl->should_now_run_cond = false; // condensation/evaporation @@ -305,15 +307,9 @@ namespace libcloudphxx 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_cond && !pimpl->opts_init.turb_cond_switch) - throw std::runtime_error("turb_cond_swtich=False, but turb_cond==True"); - if(opts.turb_adve && pimpl->n_dims==0) throw std::runtime_error("turbulent advection does not work in 0D"); - if(opts.turb_cond && pimpl->n_dims<2) - throw std::runtime_error("turbulent condensation works only in 2D and 3D"); - if (opts.chem_dsl) { // saving rv to be used as rv_old diff --git a/tests/python/unit/sstp_cond.py b/tests/python/unit/sstp_cond.py index 5f10f4f51..72ecfeadc 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 @@ -25,6 +25,7 @@ def lognormal(lnr): 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 = True spinup = 20 backend = lgrngn.backend_t.serial @@ -33,40 +34,47 @@ def lognormal(lnr): 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.turb_cond=1 opts_init.nx = 2 opts_init.dx = 1 opts_init.x1 = opts_init.nx * opts_init.dx +opts_init.nz = 1 +opts_init.dz = 1 +opts_init.z1 = opts_init.nz * opts_init.dz + +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 = arr_t([300., 300.]) - rv = arr_t([ .0025, .0095]) # first cell subsaturated, second cell supersaturated - prtcls.init(th, rv, rhod, Cx=C) + 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 = frombuffer(prtcls.outbuf()).copy() + 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 - prtcls.step_sync(opts, th, rv) + prtcls.step_sync(opts, th, rv)#, diss_rate=eps) for it in range(spinup): prtcls.step_async(opts) - prtcls.step_sync(opts, th, rv) + prtcls.step_sync(opts, th, rv)#, diss_rate=eps) #equilibrium wet moment post spinup prtcls.diag_all() prtcls.diag_wet_mom(3); - wet_post_spin = frombuffer(prtcls.outbuf()).copy() + 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 @@ -77,27 +85,25 @@ def lognormal(lnr): #equilibrium wet moment post adve prtcls.diag_all() prtcls.diag_wet_mom(3); - wet_post_adve = frombuffer(prtcls.outbuf()).copy() + 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] = tmp[1] - rv[1] = tmp[0] + rv[0,0] = tmp[1,0] + rv[1,0] = tmp[0,0] #advect th tmp = th.copy() - th[0] = tmp[1] - th[1] = tmp[0] + th[0,0] = tmp[1,0] + th[1,0] = tmp[0,0] #condensation with advected SDs and rv - prtcls.step_sync(opts, th, rv) + prtcls.step_sync(opts, th, rv)#, diss_rate=eps) #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 + 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) From 9af29da19e272c2fa731e5783cbb893c0a0b6799 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 9 Oct 2018 16:52:06 +0200 Subject: [PATCH 23/44] run sstp_cond test also with turb_cond --- tests/python/unit/sstp_cond.py | 183 ++++++++++++++++++--------------- 1 file changed, 98 insertions(+), 85 deletions(-) diff --git a/tests/python/unit/sstp_cond.py b/tests/python/unit/sstp_cond.py index 72ecfeadc..777f060c6 100644 --- a/tests/python/unit/sstp_cond.py +++ b/tests/python/unit/sstp_cond.py @@ -15,95 +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 -#opts_init.turb_cond_switch = True -spinup = 20 - -backend = lgrngn.backend_t.serial - -opts = lgrngn.opts_t() -opts.sedi=0 -opts.coal=0 -opts.cond=1 -#opts.turb_cond=1 - -opts_init.nx = 2 -opts_init.dx = 1 -opts_init.x1 = opts_init.nx * opts_init.dx - -opts_init.nz = 1 -opts_init.dz = 1 -opts_init.z1 = opts_init.nz * opts_init.dz - -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 +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)#, diss_rate=eps) - for it in range(spinup): - prtcls.step_async(opts) - prtcls.step_sync(opts, th, rv)#, diss_rate=eps) + backend = lgrngn.backend_t.serial - #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 + 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 = 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) + opts_init.nz = 1 + opts_init.dz = 1 + opts_init.z1 = opts_init.nz * opts_init.dz - #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] - - #condensation with advected SDs and rv - prtcls.step_sync(opts, th, rv)#, diss_rate=eps) + 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 = 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) + #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) From 0a944b91ac885304ff81f8e49b9b0dd92f2f4889 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 10 Oct 2018 12:59:10 +0200 Subject: [PATCH 24/44] lgrngn_subsidence test: increase tolerance a bit --- tests/python/unit/lgrngn_subsidence.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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) From ae2a10e8ac2649eef0ded12684a9d67998774f4c Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 10 Oct 2018 13:00:18 +0200 Subject: [PATCH 25/44] multi_gpu copy: add an ssp alias --- .../particles_multi_gpu_impl_step_async_and_copy.ipp | 1 + 1 file changed, 1 insertion(+) 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 1332db391..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 @@ -69,6 +69,7 @@ namespace libcloudphxx 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); From b46ea0237e6da129e1292e1374d9d5d06c984b7d Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 10 Oct 2018 15:09:28 +0200 Subject: [PATCH 26/44] enable turb cond in 0D and 1D --- src/impl/particles_impl_hskpng_remove.ipp | 42 +++++++++---------- src/impl/particles_impl_hskpng_resize.ipp | 10 ++--- src/impl/particles_impl_init_hskpng_npart.ipp | 10 ++--- src/impl/particles_impl_rcyc.ipp | 10 ++--- src/particles_step.ipp | 3 -- tests/python/unit/api_lgrngn.py | 34 +++++++++++++++ 6 files changed, 62 insertions(+), 47 deletions(-) diff --git a/src/impl/particles_impl_hskpng_remove.ipp b/src/impl/particles_impl_hskpng_remove.ipp index 08e20c25e..bcc9b00d7 100644 --- a/src/impl/particles_impl_hskpng_remove.ipp +++ b/src/impl/particles_impl_hskpng_remove.ipp @@ -20,43 +20,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.push_back(&up); - if (opts_init.ny != 0) real_t_vctrs.push_back(&vp); - if (opts_init.nz != 0) real_t_vctrs.push_back(&wp); + 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); } - else if(opts_init.turb_cond_switch) - if (opts_init.nz != 0) real_t_vctrs.push_back(&wp); if(opts_init.turb_cond_switch) { - if (opts_init.nz != 0) - { - real_t_vctrs.push_back(&ssp); - real_t_vctrs.push_back(&dot_ssp); - } + 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 37e71d1a6..0a74ac83d 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -39,16 +39,12 @@ namespace libcloudphxx if (opts_init.ny != 0) vp.resize(n_part); if (opts_init.nz != 0) wp.resize(n_part); } - else if(opts_init.turb_cond_switch) - if (opts_init.nz != 0) wp.resize(n_part); if(opts_init.turb_cond_switch) { - if (opts_init.nz != 0) - { - ssp.resize(n_part); - dot_ssp.resize(n_part); - } + 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) diff --git a/src/impl/particles_impl_init_hskpng_npart.ipp b/src/impl/particles_impl_init_hskpng_npart.ipp index 5c25afab4..fdb1eb9a8 100644 --- a/src/impl/particles_impl_init_hskpng_npart.ipp +++ b/src/impl/particles_impl_init_hskpng_npart.ipp @@ -30,16 +30,12 @@ namespace libcloudphxx 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.); } - else if(opts_init.turb_cond_switch) - if (opts_init.nz != 0) wp.resize(opts_init.n_sd_max, 0.); if(opts_init.turb_cond_switch) { - if (opts_init.nz != 0) - { - ssp.resize(opts_init.n_sd_max, 0.); - dot_ssp.resize(opts_init.n_sd_max, 0.); - } + 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 diff --git a/src/impl/particles_impl_rcyc.ipp b/src/impl/particles_impl_rcyc.ipp index f340d72ae..fe5d5e3c8 100644 --- a/src/impl/particles_impl_rcyc.ipp +++ b/src/impl/particles_impl_rcyc.ipp @@ -119,16 +119,12 @@ namespace libcloudphxx 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); } - else if(opts_init.turb_cond_switch) - if (opts_init.nz > 0) detail::copy_prop(wp.begin(), sorted_id, n_flagged); if(opts_init.turb_cond_switch) { - if (opts_init.nz > 0) - { - detail::copy_prop(ssp.begin(), sorted_id, n_flagged); - detail::copy_prop(dot_ssp.begin(), sorted_id, n_flagged); - } + 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); } { diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 41c8a515d..ce3b92c17 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -152,9 +152,6 @@ namespace libcloudphxx if(opts.turb_cond && !pimpl->opts_init.turb_cond_switch) throw std::runtime_error("turb_cond_swtich=False, but turb_cond==True"); - if(opts.turb_cond && pimpl->n_dims<2) - throw std::runtime_error("turbulent condensation works only in 2D and 3D"); - pimpl->should_now_run_cond = false; // condensation/evaporation diff --git a/tests/python/unit/api_lgrngn.py b/tests/python/unit/api_lgrngn.py index cdccdd6f7..ef366f5c6 100644 --- a/tests/python/unit/api_lgrngn.py +++ b/tests/python/unit/api_lgrngn.py @@ -327,6 +327,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) @@ -370,6 +392,8 @@ def lognormal(lnr): 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) @@ -384,6 +408,8 @@ def lognormal(lnr): opts_init.turb_adve_switch=False opts.turb_adve=False +opts_init.turb_cond_switch=False +opts.turb_cond=False # ---------- # 2D (periodic horizontal domain) @@ -437,6 +463,8 @@ def lognormal(lnr): 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) @@ -452,6 +480,8 @@ def lognormal(lnr): opts_init.turb_adve_switch=False opts.turb_adve=False +opts_init.turb_cond_switch=False +opts.turb_cond=False # ---------- @@ -493,6 +523,8 @@ def lognormal(lnr): 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) @@ -507,6 +539,8 @@ def lognormal(lnr): 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" From a81c95c37b2940f60b3d2f235b43b84b875be62b Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 11 Oct 2018 13:02:47 +0200 Subject: [PATCH 27/44] hskpng_remove: include --- src/impl/particles_impl_hskpng_remove.ipp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/impl/particles_impl_hskpng_remove.ipp b/src/impl/particles_impl_hskpng_remove.ipp index bcc9b00d7..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 { From ae6339b3baab717e307e99318b5029fa571fb7bc Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 11 Oct 2018 15:12:35 +0200 Subject: [PATCH 28/44] hlpr function for cond_sstp to shorten the if/else a bit --- src/impl/particles_impl_cond_sstp.ipp | 194 +++++++++++++------------- 1 file changed, 96 insertions(+), 98 deletions(-) diff --git a/src/impl/particles_impl_cond_sstp.ipp b/src/impl/particles_impl_cond_sstp.ipp index 64dc562a5..4d9a3a53f 100644 --- a/src/impl/particles_impl_cond_sstp.ipp +++ b/src/impl/particles_impl_cond_sstp.ipp @@ -28,6 +28,45 @@ namespace libcloudphxx }; }; + 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, @@ -75,22 +114,8 @@ namespace libcloudphxx ); } - 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 - // TODO: these function calls only slightly differ... if(!const_p) { // particle-specific pressure iterator, used twice @@ -105,93 +130,66 @@ namespace libcloudphxx ); if(turb_cond) - 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 - 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(), - ssp.begin() - )), - detail::RH_sgs(opts_init.RH_formula) - ) - )), - 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(thrust::make_tuple( // input - 2nd arg - hlpr_zip_iter, - // particle-specific p - pressure_iter, - // particle-specific RH, resolved + SGS - // TODO: if turb_cond_switch=0, ssp has size=0 - it shouldnt be a problem, but maybe? - 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) - ) - )), - rw2.begin(), // output - detail::advance_rw2(dt, RH_max) - ); + 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(), + 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 { if(turb_cond) - 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 - sstp_tmp_p.begin(), - // particle-specific RH, resolved + SGS - // TODO: if turb_cond_switch=0, ssp has size=0 - it shouldnt be a problem, but maybe? - thrust::make_transform_iterator( - thrust::make_zip_iterator(thrust::make_tuple( - sstp_tmp_p.begin(), - sstp_tmp_rv.begin(), - Tp.begin(), - ssp.begin() - )), - detail::RH_sgs(opts_init.RH_formula) - ) - )), - 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(thrust::make_tuple( // input - 2nd arg - hlpr_zip_iter, - // particle-specific p - sstp_tmp_p.begin(), - // particle-specific RH, resolved + SGS - // TODO: if turb_cond_switch=0, ssp has size=0 - it shouldnt be a problem, but maybe? - 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) - ) - )), - rw2.begin(), // output - detail::advance_rw2(dt, RH_max) - ); + 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(), + 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 From 1007a91179a8a36fa115b37c054ddb15d0c98704 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 23 Oct 2018 15:25:47 +0200 Subject: [PATCH 29/44] cond_sstp_hlpr declaration in impl --- src/impl/particles_impl.ipp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 7f9f86e6d..4bc0254ed 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -456,6 +456,8 @@ namespace libcloudphxx void cond_dm3_helper(); 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 &); From 1d4734fd1d41600354bcae68619214dddc9c4a12 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 8 Apr 2019 12:26:35 +0200 Subject: [PATCH 30/44] icicle_chem: empty arrinfo for diss_rate --- models/kinematic_2D/src/kin_cloud_2d_lgrngn_chem.hpp | 1 + 1 file changed, 1 insertion(+) 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 ); From cde9d09ab0a0ba401ad0f9a5add943d482c29a35 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 9 Apr 2019 13:33:02 +0200 Subject: [PATCH 31/44] add turb_coal_switch and its sanity checks --- include/libcloudph++/lgrngn/opts_init.hpp | 3 ++- src/impl/particles_impl_init_kernel.ipp | 4 ++++ src/impl/particles_impl_init_sync.ipp | 2 +- src/particles_step.ipp | 14 +++++++------- 4 files changed, 14 insertions(+), 9 deletions(-) diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index 0ae38b549..1c7cedb17 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -104,7 +104,8 @@ namespace libcloudphxx 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 turb_adve_switch, // if true, turbulent motion of SDs is modeled - turb_cond_switch; // if true, turbulent condensation 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; diff --git a/src/impl/particles_impl_init_kernel.ipp b/src/impl/particles_impl_init_kernel.ipp index 7e7efdfb2..59884629e 100644 --- a/src/impl/particles_impl_init_kernel.ipp +++ b/src/impl/particles_impl_init_kernel.ipp @@ -186,6 +186,8 @@ namespace libcloudphxx { throw std::runtime_error("Please supply two kernel parameters: rate of dissipation epsilon [m^2/s^3] and 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); @@ -209,6 +211,8 @@ namespace libcloudphxx { throw std::runtime_error("Please supply two kernel parameters: rate of dissipation epsilon [m^2/s^3] and 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 2c4743b51..616792eab 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/src/impl/particles_impl_init_sync.ipp @@ -20,7 +20,7 @@ namespace libcloudphxx 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) + 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) diff --git a/src/particles_step.ipp b/src/particles_step.ipp index d9356c0ca..e1b166818 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -71,11 +71,11 @@ 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) && diss_rate.is_null()) - throw std::runtime_error("turbulent advection 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 not switched off and diss_rate is empty"); - if ( !(pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_switch) && !diss_rate.is_null()) - throw std::runtime_error("turbulent advection and condesation are switched off and diss_rate 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 switched off and diss_rate is not empty"); // if (pimpl->l2e[&pimpl->courant_x].size() == 0) // TODO: y, z,... @@ -112,13 +112,13 @@ namespace libcloudphxx 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) + 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) + 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) @@ -357,7 +357,7 @@ namespace libcloudphxx if (opts.turb_adve || opts.turb_cond) { - // calc tke (diss_rate now holds TKE, not dissipation rate!) + // 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) From 953579e95ce2dc1c766a543afb730ed9d029b430 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 9 Apr 2019 14:13:13 +0200 Subject: [PATCH 32/44] pass tke diss rate to onishis kernel --- src/detail/kernels.hpp | 7 ++++--- src/detail/tpl_calc_wrapper.hpp | 3 ++- src/impl/particles_impl_coal.ipp | 16 +++++++++++----- src/impl/particles_impl_init_kernel.ipp | 8 ++++---- 4 files changed, 21 insertions(+), 13 deletions(-) diff --git a/src/detail/kernels.hpp b/src/detail/kernels.hpp index 553694dea..3af796375 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,11 +218,12 @@ 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 ); + printf("onishi calc diss rate: %g\n", thrust::get(tpl_wrap.get_ro_calc())); real_t res = kernel_geometric::interpolated_efficiency(rwa, rwb) * // stagnant air collision efficiency wang_collision_enhancement(rwa, rwb, kernel_base::k_params[0]) * // Wang turbulent collision efficiency enhancement, k_params[0] - epsilon diff --git a/src/detail/tpl_calc_wrapper.hpp b/src/detail/tpl_calc_wrapper.hpp index da8d6a305..839c564e2 100644 --- a/src/detail/tpl_calc_wrapper.hpp +++ b/src/detail/tpl_calc_wrapper.hpp @@ -21,7 +21,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/impl/particles_impl_coal.ipp b/src/impl/particles_impl_coal.ipp index 84a800f5a..5f5ab0055 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; @@ -334,9 +335,10 @@ namespace libcloudphxx typedef thrust::zip_iterator< thrust::tuple< pi_real_t, // rhod - pi_real_t // eta + pi_real_t, // eta + pi_real_t // tke dissipation rate > - > zip_ro_calc_t; //read-only parameters passed to the calc() function, later also epsilon and Re_lambda + > zip_ro_calc_t; //read-only parameters passed to the calc() function, TODO: add Re_lambda for Onishi kernel zip_ro_t zip_ro_it( thrust::make_tuple( @@ -360,7 +362,11 @@ namespace libcloudphxx // rhod thrust::make_permutation_iterator(rhod.begin(), sorted_ijk.begin()), // eta - thrust::make_permutation_iterator(eta.begin(), sorted_ijk.begin()) + thrust::make_permutation_iterator(eta.begin(), sorted_ijk.begin()), + // tke dissipation rate + opts_init.turb_coal_switch ? + thrust::make_permutation_iterator(diss_rate.begin(), sorted_ijk.begin()) : + thrust::make_permutation_iterator(thrust::make_constant_iterator(0), sorted_ijk.begin()) ) ); diff --git a/src/impl/particles_impl_init_kernel.ipp b/src/impl/particles_impl_init_kernel.ipp index 59884629e..9b4ab8eee 100644 --- a/src/impl/particles_impl_init_kernel.ipp +++ b/src/impl/particles_impl_init_kernel.ipp @@ -182,9 +182,9 @@ 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"); @@ -207,9 +207,9 @@ 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"); From cab2c42ea0dcc9155aed10c79cf2d1034b136005 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 10 Apr 2019 17:29:25 +0200 Subject: [PATCH 33/44] actually include libcloudphxx_INCLUDE_DIRS when building libcloud... --- CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) 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) From 80dd14647349ce2a224c33ed1c4096077454f91f Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 11 Apr 2019 14:04:34 +0200 Subject: [PATCH 34/44] WIP --- src/detail/kernels.hpp | 2 +- src/detail/tpl_calc_wrapper.hpp | 2 ++ src/impl/.particles_impl_coal.ipp.swp | Bin 0 -> 36864 bytes src/impl/particles_impl_coal.ipp | 38 +++++++++++++++++++------- 4 files changed, 31 insertions(+), 11 deletions(-) create mode 100644 src/impl/.particles_impl_coal.ipp.swp diff --git a/src/detail/kernels.hpp b/src/detail/kernels.hpp index 3af796375..fa8b8bdfb 100644 --- a/src/detail/kernels.hpp +++ b/src/detail/kernels.hpp @@ -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[0], thrust::get(tpl_wrap.get_ro_calc()) // k_params[0] - 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 839c564e2..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 diff --git a/src/impl/.particles_impl_coal.ipp.swp b/src/impl/.particles_impl_coal.ipp.swp new file mode 100644 index 0000000000000000000000000000000000000000..bbd432fc2d96a357b82ae9fbbb357295cb16083b GIT binary patch literal 36864 zcmeI536LCDdB?}*G++q_bDQExfS8q5yIKhnSV@cV5yuXe93K#4@X$NcyW42CdpzB< zDHBykRMlwBA|1sot{Ly83BP=UlT5aNJ|1E#P!li-5n_r0V0bk5y@(=&>2p5x zsM!wFhWu5pAE`Wlp#2@#?}yp~TsaT2_fNLpQRVqV?frA@_u8`C)R2S_{<|}wWot8fs?_b zzz-i@tKAI#4txN-17zT(;GZ8>t9=vbBL$~x^yZpCS)Vt#JpX(LW1>@yqw- zEx0_tQmPwfy>6`Hxf%5wRW}9J!^zI2Yt?IFS&#GN}mG>g$r(p=1S(L=; znNu@Udmxv(U@2M-RG2BkQ`GKVyPloD<2q6bm3FD56Ly0v>DDLeUMI$`|CAcmgD$5u z>-L(NN}EADY%By#^0K9^lbCZOsN+*gsq!zHLDwhkmfEpjWlP;&n$6BGhS|P7b9ESV zVy}`eYPEE-`J4($e|XgEw4}FPCc#2R-)a6%A!sGNd8_YleL%Tths%>H_#evD%xV+N zH#*&LB}#f}%}+J-+AnICX`JoTiS4P+gjE9$~ty3_ri7>$`zTy+)3P07o4Gz1-dv?ognHm)>5?y3rl}_W9>-I zyl?)llUH^R;FX}1k*>%w=M`J-C|(JU0k0IDXIWijivQ{Fz4h>l%eKlZ1#9}eQnWWO zE$7Uk%Ph5VPV&oJKvc?v)Qt8c zi<{9tzlG0moU~&{GlIm^s+}Yqm8|$uB~|lbwin=H{xz@_JU&4S6WLo4y zqgn}5tq6f3UC9JINGd$1Chwb%o%agz)Z7zYf|JK;C203TkvXmI5!DS^Q7?58D2YS0 z8f8l=T@KpqiUeZ4lBre16OSmXb;M3GEt#0@(#h}e3i~?x?%S_+P1u!zp2b>;ij!<9 ziWgP1z-2K_QLnBgh}5jZk|b=IDhRXdng+4z#o-D}7x$LusWO!jIyS?0J10QF=X#yz zEe)7U%~Eglx`58Ki>Q6gHdT=YVk2|wJy(QXHLQj#m$IASoi&-NMdRx+oJO6;3~R z-oeOP~n@YNi#K$+oI1STjkqFz>sZIxZ-r1+0dEkOG z&CO=gIntF8t}corCc!kyj))L_LC{Tu)!OuSUX(Iwi1w9ec~-7P@!H=cre5^__n?n| z2Z;WUimpE&K$jOi|0M8N-25Z(CXj+N!7Mld90%^^-S68x?vi}QpD_i-6c|%rOo1^4 z#uOM+U`&BA1;!K@Q(#Pi|9uK*_4ac7B`A-Dqo*+9$3fg&3Y&-OhdTJ~w(=L-WGBI- zEBmm+6os@IZk^bxUzaSrOJ%j>+1Wmm(>`lmsu#8`UY)?|9yesaH;T9O_erKo3@=mj z;bIimCw9!#^ey|xds2SInN*g8Ltz7t+j1`p#AVlraA0Wn8}>+->7+`Nt{BClmb*nk z$2h}{W(T!hWsXh#Mr1gt_5UZql^??&E&6|8IrM|*{eK4D2`&bwfZNdZuLjPS|01q; zfoX6PdjBiIGk^jQ1Yg0Y|A(Lh_J9eXz;WPT@#%jSd=k7Kybe4OJOCUI?!?D`19%lU z4NQY6a6I@le*f2l9>~CUFadsy-~X?{YrzWWfnNd726qw%@CV?9U^BRmIDjr#0?!8* zff?|4a0f8~*MT?M_<%{iw}MB4&EO8=1HJ(M9J~u00s(Nb0-qx$;6val@J4VHOoCIu zBfyV{2lxhfH#i5(f^S3r+rXE=7lHK0he74zj{|od^e~UUEb276Ndy0I6L!atC59*Y zc*D~ldZWimYo(tRUwX*h^Yapy!Aeb+rvHr2J6Tw885Ml)WnLbLKYf0aM%tMwQr%|M z^6LAtu$z=ilso>|<(25L_r~ZNrOjM?4QyexMx^)h@9mn&>Eqp|8a&Bwlxjg>M?0+B zXO#}-jmclMCjv8dJ`QVqm-O`swa|-o9K@u$ z@ZcE-XSLt|{IJmumgif6!f~e##!mNhcK^3X&QC3 zxW*5Zg|;E3)RQ_jcG%$$yGu#YvyMKRH2dsJV)<-^;^a1Q3_>b5KkTazX|<6R^PJxR zSAX14MN6b9;x1s?l&^W8=kgX$7LA{j7fkSteED9&W0kyX;&bGylD)*&SmL@%YNO+I zBy6Wh?dW8P=-RXqDkj6In72|YvXKupEeRYlN~4#B(j@P49XXN_qKMFqZTh=tb>Ts3 zCcQYrA6PJnp@=PUPE*k@iaevk9sz&qtb$sReNM8zVv95buc%d7#5l&?Xhy=cgM`pM z;t!ew_Z4F-j|a0fZV|TA3}HVVK|eslOtGqmD<(Nwb zm$-;j8(dYp+$~DWCB@!C$mG3wpHZN%TK9+ubHy;b&uU&Cbp>5F)&Mg$U|}RQ;gyVQ z*Bjj+ZZ%-+y3t4SWfN($gt>``0P|#eHxpFcN|vpYKh=|7Gg)TllvFwrpwms}gZZd! zU2(cJ33f{*Bp@)+RXXY+Vb;zFU8ETO$<;6966zO!QNLPUos4%|;bC2X&PBos3Qa7p zu@YJ4mJj*R@a|3O1Yz8&Pvq3Ivtf9++3uxy5gUXW#&sv|y}bQ}pO%e{s93m}HkAI- zFmFKR)hShSs0`v|LI;EVEyDz}v%@t!DVw*FUe@Vl`Fp0KczzvO&^zrtZz0x~Du3q6 z6<4cUW8qb|FKpQ#H>59!fN}}(v9{|&NN$5d>$|2Va=KDwdBrf#F1SR3Pb9Lmt%FBW z9&kAy`e80QLI*ZZR)y`PYPP*(*6b0l$7+gD9US2Fw`(pxcFw7K1Xq!UXLY)-9!;Q% zFOoW~@S5f`c|cT+1hLo>axP^~flK*IemXc+nc&pa=S$_-2yGyDj#zO#=;vtNls^23byGW&`D z{|)HlKLDctyZ!gKq3{1axDs^1Q^E1z9q9TJ2k;v3YH%^w4|aerpGaIkxEx#p9s+)B z^Zp0E7r=?&>*)8_fJ?!t;JfJbSAjQ!2G{}~3pRlh!4J^w-wCF`7trND1KtY~a4L8d zxDB2DYv6kDd*Iz511|*!z{%hSbo&kvUH{o&0^Er{|6cIRU>7(U{3p8nE#OAb0Sn+# za1J;g{0JTXR`7msD!3E<{q^7%z>Vna*Mi>%ZQy?#rxapB%Q0UA$%8EQlafv+Z%^a` z+vLT(;zq$5MzW&Zjb4c)k%>~KF1Q=*j`*=^%(i3F;X?U7H@itYsk^a4igdJ6spb=4 zbiGGd&&oD~O0Rm8UamTbP1W^RVHGIsk6JWPZygm$OkW}rpD47j5GR#V*3CPvv2n7=xC|0QbJJ`eEKbG71!#OqM1Ut2tusI!2 z?85B}<8pB&F;vFT++a1e7}ha5FE{2F<4d9!qdylvJPWa!VXl!!h}g9((j}vp^tg%L znw@pAU?n58j4H2?iV1ivLv-nxTS++%a|9ngD!fJxVumVrReF9Pe~rkpT*p!k*Ud=M zzN#215uzf9O3%s_I0h}c-xA3Yw>sWT{GeeHp6Mue9poPE)1AkPZ?D?1BXZ+u;P})j zQxuKnl1pdD-Lixw>kw*56JfKrB2%^yOB}Q^S+futHYUW0J?i6ANz^8mSxhdJ1cWUF(A)-V-}Oy(QM7hOts|)$4|w zG|UFOhTiFdWiM``EvINgY1Yw)5Ib5|_;%FLk>{A^Bg)g4tHot*urG-q)0qnClphkylq3|wQ6FQ7;eIl#*?C=#8rqf zPylXdrQrxLTttt^Zn98pbcMZ549(vB@L79c)mJi)J_CJ|$}eto<%rFSzy;`c19!!B~cqUt=Mps~i8C$n*?pxk9ArtA^#Yr94YG zU4@m9Pr4$9%4Q3AG^7s0wW|+Z}?Pso=ME^(a*ZQhlxbyw~7G3`>;LYF?uooN$ zK98P%6ZkmzDEJ5vTfiY8Hh}HmG2jk#|8Ihife(WZf!BZpEP}H@9ee|Oz%}5NU=ExF zer)XnSMz-ZcnR1FzJe{_x4^~VZfpU!gTr7I>;{;(Fh^SHDwV6QH|0$4;Jxdxgo`h^@Y_ucrG9p|hf8bKpE%ob2@6Q$YVB{lD^vXWvm zrrvfh7S?iRN4=^jF|yX1`04_+Wj`SSIPF+j=9;W2T$Ul?7PlsNF~rLDKuwBN4pi{X zqUIGXf~+Oj`l>jvL1$myJ}#<|&lFg)o6A4qNp`1DIkvA^V^OZqmbau$gmWM-SHeKqfZVTV|-#LQ3pu)9@@jPm`aTEUPbPrCMu3^Rz; zHc`}!Bx%f}k2wjZwK%oCP~i-{<`Wd#er-_>X%>2J8H$uCU)}hefs$!D2$8&9|qEhMVG1FK(~oo>Z*F{2pce?CL0{nW+FQv z{lPja2Wcc(ARXn@ixyihOBpzyzmpjz)S>1YLklQiJ>gS%@Khg}3B1wNCGtjp>%g7- z%G{GsA$?1Nh0HTKoF~6B_uLF=gOGSSpqEO4O|FA-3R`yY+BZm7S9PPGbM}g!2u+Lz zNA*Z}$d;TsLvqgx6Y60u($u#H0BGHLjsB~2CploA;v5wsexo>OtCg^wa3s)LEYP2Y z7@_K`oaSRIkzUpdx?#(&25Y9g+7PGmlWN4=t@{V16mPGShgQ49wD>W2vX;^37D)zu zYy+ir_v78u-*&b`B$=g!_Vim8am$HGVy6akMJZ6r^&3qZ^r5w-5G*q`@-VE*Pf}An zZf8G$K~_!$ft(#;jVwwbn3W!^)R6+k9bLG7YDI#PZ27AV z=Q23EqN~Cix;Wo@QRSz$lQZsv)Ui**QlIE~Xqe`YX>$|6>_e{v|N%_-B%IhQ@0BLGC z>E75g#bICi71ySLui1uGUc+HjX^(_h<^+}Ub*nsq(p+Dwd?GO`*erJzlrIsO`NE_y zedfy^7N#B=>M}QoWCJ@C;IU=KJ6Nc{gz;J3kd(d+*W z+zdVg-VL4!P6ppaxBnct8oUWy0zQCVFR}iczyrbE==0wP-vieJISX(FG{9rQ=h5qL z1aATvm;kq<&tC_w1!sd_0(I~J@NsndrvcIDBd{540&>2e(DAK8%X|DWf{KxSOap@2 z!4Vn9OxdM49h6)2bV7hQ(TA9LD~XiahPbY0%bkXaz%)Wx3v7LQ3r|y~N*q`UVWK63 z5$#y>-G6U%nzB_bO|uMl+gyjlb*j7;=AFTRAvLyleasiI|%wAUPgH@MiIoD+l$(+)Y2&6DcD!3)ov#YXaUbIg)n zKAClswqIjBk;AO?T0l-Xq8K^gNb)I7O@+H;8w-HGis|msLFRsBmhLoTHfrnmL`&_# zJ!7Qnoa+*Z-VGPRuIvHR#vWM2MCjREoFmXV^>v6>Dfp=`9tu%vXEGLZy z$q-9L=v6tGvsy{f09g*K!!qhlJQ`X}YbG38&kaaTsrMSSGKwmP*%np_d4o@`=gO7n zxix(||D#uzTJO;-P`Rq6kX4yUAv<9Pq|lpPbPBoZPeBigO)BR}QTi%{m=0%cjP_Z# z3uGN*Q0fVz+3UrcSgAP~((flzbpEsm3^t9eoZwYSt>(n^>($EX8S=6m^65P)Y*6yh zuYxpV0}SL&ebl000H5Da7|57M94mQ*kml!XqpvihR8u&%QCfu<1^96LH!Nn@5_W-+LV<6KB|!sZ>(wgBOb}>q^dL zB%cM31Hm?Ln3sS&QujQ9RTI;FVV%OM{$OKHWApzWQ8O)vPQttV7mUvhsz4 z4kFhIJuZ>545@azF*RBCOD5DlC5uR0R?;q{H$iLW%4a-!q}UiE9n$OQ^;}E8>^8{g z9Wz0`kT6GQii6RRRJw)Q^kUO-LpxxAbzdA+-CpQdsnYzaIbUttW?Ljois8pA4%zRX zaqMj)s|yA63_57m?mK{XuVQ~tCa0u&dCIYV&fTzcPTaKG87uDA$`CG}rD-gFd3Sg@ znkBpGo&ESWDuO~A|BClS|GyDkRs7JR|DS91*8fJ=|2p^^a1C&00e%NP|0-}k*aN)2HJoq%Y8XN>)L$7}wcsV!< z4uB_vyU^`#0B-}o3A*5Xa1MAXcmntd`u-hy z-R{Mdc4G5v}eg%E@QHgwgJImZtR?nsV zLF4@jNB}qV>nMT8Zofj!$(GVW@7~+AP^kOsn#OJD?mgSKaIchA7+!{PYsuL0f0eAb z`_ay`p<6(0), sorted_ijk.begin()) + thrust::make_permutation_iterator(thrust::make_constant_iterator(0), sorted_ijk.begin()) ) - ); + ; + + zip_ro_calc_t zip_ro_calc_turb_it = + 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( @@ -390,11 +399,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(opts.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"); From a41f5287236a5618392175e9149eedd4b5dda2af Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 11 Apr 2019 14:14:50 +0200 Subject: [PATCH 35/44] opts_init.rd_min to set minimal dry radius --- include/libcloudph++/lgrngn/opts_init.hpp | 5 ++++- src/impl/particles_impl_dist_analysis.ipp | 8 ++++++++ 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index 0ae38b549..b560dd9af 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -124,6 +124,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), @@ -156,7 +158,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/src/impl/particles_impl_dist_analysis.ipp b/src/impl/particles_impl_dist_analysis.ipp index 23797f6a8..a9cdda184 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, 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, log(opts_init.rd_min)); // user-defined lower limit for rd } }; }; From 259009f4c3eee31d37a444f0784833ab023172ee Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 11 Apr 2019 14:36:06 +0200 Subject: [PATCH 36/44] workin turb coal diss_Rate --- include/libcloudph++/lgrngn/opts.hpp | 4 +-- src/impl/.particles_impl_coal.ipp.swp | Bin 36864 -> 0 bytes src/impl/particles_impl.ipp | 2 +- src/impl/particles_impl_coal.ipp | 46 ++++++++++++-------------- src/particles_step.ipp | 2 +- 5 files changed, 25 insertions(+), 29 deletions(-) delete mode 100644 src/impl/.particles_impl_coal.ipp.swp diff --git a/include/libcloudph++/lgrngn/opts.hpp b/include/libcloudph++/lgrngn/opts.hpp index b0ede176e..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, turb_adve, turb_cond; + bool adve, sedi, cond, coal, src, rcyc, turb_adve, turb_cond, turb_coal; // RH limit for drop growth real_t RH_max; @@ -32,7 +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_adve(false), turb_cond(false), turb_coal(false), RH_max(44) // :) (anything greater than 1.1 would be enough { } diff --git a/src/impl/.particles_impl_coal.ipp.swp b/src/impl/.particles_impl_coal.ipp.swp deleted file mode 100644 index bbd432fc2d96a357b82ae9fbbb357295cb16083b..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 36864 zcmeI536LCDdB?}*G++q_bDQExfS8q5yIKhnSV@cV5yuXe93K#4@X$NcyW42CdpzB< zDHBykRMlwBA|1sot{Ly83BP=UlT5aNJ|1E#P!li-5n_r0V0bk5y@(=&>2p5x zsM!wFhWu5pAE`Wlp#2@#?}yp~TsaT2_fNLpQRVqV?frA@_u8`C)R2S_{<|}wWot8fs?_b zzz-i@tKAI#4txN-17zT(;GZ8>t9=vbBL$~x^yZpCS)Vt#JpX(LW1>@yqw- zEx0_tQmPwfy>6`Hxf%5wRW}9J!^zI2Yt?IFS&#GN}mG>g$r(p=1S(L=; znNu@Udmxv(U@2M-RG2BkQ`GKVyPloD<2q6bm3FD56Ly0v>DDLeUMI$`|CAcmgD$5u z>-L(NN}EADY%By#^0K9^lbCZOsN+*gsq!zHLDwhkmfEpjWlP;&n$6BGhS|P7b9ESV zVy}`eYPEE-`J4($e|XgEw4}FPCc#2R-)a6%A!sGNd8_YleL%Tths%>H_#evD%xV+N zH#*&LB}#f}%}+J-+AnICX`JoTiS4P+gjE9$~ty3_ri7>$`zTy+)3P07o4Gz1-dv?ognHm)>5?y3rl}_W9>-I zyl?)llUH^R;FX}1k*>%w=M`J-C|(JU0k0IDXIWijivQ{Fz4h>l%eKlZ1#9}eQnWWO zE$7Uk%Ph5VPV&oJKvc?v)Qt8c zi<{9tzlG0moU~&{GlIm^s+}Yqm8|$uB~|lbwin=H{xz@_JU&4S6WLo4y zqgn}5tq6f3UC9JINGd$1Chwb%o%agz)Z7zYf|JK;C203TkvXmI5!DS^Q7?58D2YS0 z8f8l=T@KpqiUeZ4lBre16OSmXb;M3GEt#0@(#h}e3i~?x?%S_+P1u!zp2b>;ij!<9 ziWgP1z-2K_QLnBgh}5jZk|b=IDhRXdng+4z#o-D}7x$LusWO!jIyS?0J10QF=X#yz zEe)7U%~Eglx`58Ki>Q6gHdT=YVk2|wJy(QXHLQj#m$IASoi&-NMdRx+oJO6;3~R z-oeOP~n@YNi#K$+oI1STjkqFz>sZIxZ-r1+0dEkOG z&CO=gIntF8t}corCc!kyj))L_LC{Tu)!OuSUX(Iwi1w9ec~-7P@!H=cre5^__n?n| z2Z;WUimpE&K$jOi|0M8N-25Z(CXj+N!7Mld90%^^-S68x?vi}QpD_i-6c|%rOo1^4 z#uOM+U`&BA1;!K@Q(#Pi|9uK*_4ac7B`A-Dqo*+9$3fg&3Y&-OhdTJ~w(=L-WGBI- zEBmm+6os@IZk^bxUzaSrOJ%j>+1Wmm(>`lmsu#8`UY)?|9yesaH;T9O_erKo3@=mj z;bIimCw9!#^ey|xds2SInN*g8Ltz7t+j1`p#AVlraA0Wn8}>+->7+`Nt{BClmb*nk z$2h}{W(T!hWsXh#Mr1gt_5UZql^??&E&6|8IrM|*{eK4D2`&bwfZNdZuLjPS|01q; zfoX6PdjBiIGk^jQ1Yg0Y|A(Lh_J9eXz;WPT@#%jSd=k7Kybe4OJOCUI?!?D`19%lU z4NQY6a6I@le*f2l9>~CUFadsy-~X?{YrzWWfnNd726qw%@CV?9U^BRmIDjr#0?!8* zff?|4a0f8~*MT?M_<%{iw}MB4&EO8=1HJ(M9J~u00s(Nb0-qx$;6val@J4VHOoCIu zBfyV{2lxhfH#i5(f^S3r+rXE=7lHK0he74zj{|od^e~UUEb276Ndy0I6L!atC59*Y zc*D~ldZWimYo(tRUwX*h^Yapy!Aeb+rvHr2J6Tw885Ml)WnLbLKYf0aM%tMwQr%|M z^6LAtu$z=ilso>|<(25L_r~ZNrOjM?4QyexMx^)h@9mn&>Eqp|8a&Bwlxjg>M?0+B zXO#}-jmclMCjv8dJ`QVqm-O`swa|-o9K@u$ z@ZcE-XSLt|{IJmumgif6!f~e##!mNhcK^3X&QC3 zxW*5Zg|;E3)RQ_jcG%$$yGu#YvyMKRH2dsJV)<-^;^a1Q3_>b5KkTazX|<6R^PJxR zSAX14MN6b9;x1s?l&^W8=kgX$7LA{j7fkSteED9&W0kyX;&bGylD)*&SmL@%YNO+I zBy6Wh?dW8P=-RXqDkj6In72|YvXKupEeRYlN~4#B(j@P49XXN_qKMFqZTh=tb>Ts3 zCcQYrA6PJnp@=PUPE*k@iaevk9sz&qtb$sReNM8zVv95buc%d7#5l&?Xhy=cgM`pM z;t!ew_Z4F-j|a0fZV|TA3}HVVK|eslOtGqmD<(Nwb zm$-;j8(dYp+$~DWCB@!C$mG3wpHZN%TK9+ubHy;b&uU&Cbp>5F)&Mg$U|}RQ;gyVQ z*Bjj+ZZ%-+y3t4SWfN($gt>``0P|#eHxpFcN|vpYKh=|7Gg)TllvFwrpwms}gZZd! zU2(cJ33f{*Bp@)+RXXY+Vb;zFU8ETO$<;6966zO!QNLPUos4%|;bC2X&PBos3Qa7p zu@YJ4mJj*R@a|3O1Yz8&Pvq3Ivtf9++3uxy5gUXW#&sv|y}bQ}pO%e{s93m}HkAI- zFmFKR)hShSs0`v|LI;EVEyDz}v%@t!DVw*FUe@Vl`Fp0KczzvO&^zrtZz0x~Du3q6 z6<4cUW8qb|FKpQ#H>59!fN}}(v9{|&NN$5d>$|2Va=KDwdBrf#F1SR3Pb9Lmt%FBW z9&kAy`e80QLI*ZZR)y`PYPP*(*6b0l$7+gD9US2Fw`(pxcFw7K1Xq!UXLY)-9!;Q% zFOoW~@S5f`c|cT+1hLo>axP^~flK*IemXc+nc&pa=S$_-2yGyDj#zO#=;vtNls^23byGW&`D z{|)HlKLDctyZ!gKq3{1axDs^1Q^E1z9q9TJ2k;v3YH%^w4|aerpGaIkxEx#p9s+)B z^Zp0E7r=?&>*)8_fJ?!t;JfJbSAjQ!2G{}~3pRlh!4J^w-wCF`7trND1KtY~a4L8d zxDB2DYv6kDd*Iz511|*!z{%hSbo&kvUH{o&0^Er{|6cIRU>7(U{3p8nE#OAb0Sn+# za1J;g{0JTXR`7msD!3E<{q^7%z>Vna*Mi>%ZQy?#rxapB%Q0UA$%8EQlafv+Z%^a` z+vLT(;zq$5MzW&Zjb4c)k%>~KF1Q=*j`*=^%(i3F;X?U7H@itYsk^a4igdJ6spb=4 zbiGGd&&oD~O0Rm8UamTbP1W^RVHGIsk6JWPZygm$OkW}rpD47j5GR#V*3CPvv2n7=xC|0QbJJ`eEKbG71!#OqM1Ut2tusI!2 z?85B}<8pB&F;vFT++a1e7}ha5FE{2F<4d9!qdylvJPWa!VXl!!h}g9((j}vp^tg%L znw@pAU?n58j4H2?iV1ivLv-nxTS++%a|9ngD!fJxVumVrReF9Pe~rkpT*p!k*Ud=M zzN#215uzf9O3%s_I0h}c-xA3Yw>sWT{GeeHp6Mue9poPE)1AkPZ?D?1BXZ+u;P})j zQxuKnl1pdD-Lixw>kw*56JfKrB2%^yOB}Q^S+futHYUW0J?i6ANz^8mSxhdJ1cWUF(A)-V-}Oy(QM7hOts|)$4|w zG|UFOhTiFdWiM``EvINgY1Yw)5Ib5|_;%FLk>{A^Bg)g4tHot*urG-q)0qnClphkylq3|wQ6FQ7;eIl#*?C=#8rqf zPylXdrQrxLTttt^Zn98pbcMZ549(vB@L79c)mJi)J_CJ|$}eto<%rFSzy;`c19!!B~cqUt=Mps~i8C$n*?pxk9ArtA^#Yr94YG zU4@m9Pr4$9%4Q3AG^7s0wW|+Z}?Pso=ME^(a*ZQhlxbyw~7G3`>;LYF?uooN$ zK98P%6ZkmzDEJ5vTfiY8Hh}HmG2jk#|8Ihife(WZf!BZpEP}H@9ee|Oz%}5NU=ExF zer)XnSMz-ZcnR1FzJe{_x4^~VZfpU!gTr7I>;{;(Fh^SHDwV6QH|0$4;Jxdxgo`h^@Y_ucrG9p|hf8bKpE%ob2@6Q$YVB{lD^vXWvm zrrvfh7S?iRN4=^jF|yX1`04_+Wj`SSIPF+j=9;W2T$Ul?7PlsNF~rLDKuwBN4pi{X zqUIGXf~+Oj`l>jvL1$myJ}#<|&lFg)o6A4qNp`1DIkvA^V^OZqmbau$gmWM-SHeKqfZVTV|-#LQ3pu)9@@jPm`aTEUPbPrCMu3^Rz; zHc`}!Bx%f}k2wjZwK%oCP~i-{<`Wd#er-_>X%>2J8H$uCU)}hefs$!D2$8&9|qEhMVG1FK(~oo>Z*F{2pce?CL0{nW+FQv z{lPja2Wcc(ARXn@ixyihOBpzyzmpjz)S>1YLklQiJ>gS%@Khg}3B1wNCGtjp>%g7- z%G{GsA$?1Nh0HTKoF~6B_uLF=gOGSSpqEO4O|FA-3R`yY+BZm7S9PPGbM}g!2u+Lz zNA*Z}$d;TsLvqgx6Y60u($u#H0BGHLjsB~2CploA;v5wsexo>OtCg^wa3s)LEYP2Y z7@_K`oaSRIkzUpdx?#(&25Y9g+7PGmlWN4=t@{V16mPGShgQ49wD>W2vX;^37D)zu zYy+ir_v78u-*&b`B$=g!_Vim8am$HGVy6akMJZ6r^&3qZ^r5w-5G*q`@-VE*Pf}An zZf8G$K~_!$ft(#;jVwwbn3W!^)R6+k9bLG7YDI#PZ27AV z=Q23EqN~Cix;Wo@QRSz$lQZsv)Ui**QlIE~Xqe`YX>$|6>_e{v|N%_-B%IhQ@0BLGC z>E75g#bICi71ySLui1uGUc+HjX^(_h<^+}Ub*nsq(p+Dwd?GO`*erJzlrIsO`NE_y zedfy^7N#B=>M}QoWCJ@C;IU=KJ6Nc{gz;J3kd(d+*W z+zdVg-VL4!P6ppaxBnct8oUWy0zQCVFR}iczyrbE==0wP-vieJISX(FG{9rQ=h5qL z1aATvm;kq<&tC_w1!sd_0(I~J@NsndrvcIDBd{540&>2e(DAK8%X|DWf{KxSOap@2 z!4Vn9OxdM49h6)2bV7hQ(TA9LD~XiahPbY0%bkXaz%)Wx3v7LQ3r|y~N*q`UVWK63 z5$#y>-G6U%nzB_bO|uMl+gyjlb*j7;=AFTRAvLyleasiI|%wAUPgH@MiIoD+l$(+)Y2&6DcD!3)ov#YXaUbIg)n zKAClswqIjBk;AO?T0l-Xq8K^gNb)I7O@+H;8w-HGis|msLFRsBmhLoTHfrnmL`&_# zJ!7Qnoa+*Z-VGPRuIvHR#vWM2MCjREoFmXV^>v6>Dfp=`9tu%vXEGLZy z$q-9L=v6tGvsy{f09g*K!!qhlJQ`X}YbG38&kaaTsrMSSGKwmP*%np_d4o@`=gO7n zxix(||D#uzTJO;-P`Rq6kX4yUAv<9Pq|lpPbPBoZPeBigO)BR}QTi%{m=0%cjP_Z# z3uGN*Q0fVz+3UrcSgAP~((flzbpEsm3^t9eoZwYSt>(n^>($EX8S=6m^65P)Y*6yh zuYxpV0}SL&ebl000H5Da7|57M94mQ*kml!XqpvihR8u&%QCfu<1^96LH!Nn@5_W-+LV<6KB|!sZ>(wgBOb}>q^dL zB%cM31Hm?Ln3sS&QujQ9RTI;FVV%OM{$OKHWApzWQ8O)vPQttV7mUvhsz4 z4kFhIJuZ>545@azF*RBCOD5DlC5uR0R?;q{H$iLW%4a-!q}UiE9n$OQ^;}E8>^8{g z9Wz0`kT6GQii6RRRJw)Q^kUO-LpxxAbzdA+-CpQdsnYzaIbUttW?Ljois8pA4%zRX zaqMj)s|yA63_57m?mK{XuVQ~tCa0u&dCIYV&fTzcPTaKG87uDA$`CG}rD-gFd3Sg@ znkBpGo&ESWDuO~A|BClS|GyDkRs7JR|DS91*8fJ=|2p^^a1C&00e%NP|0-}k*aN)2HJoq%Y8XN>)L$7}wcsV!< z4uB_vyU^`#0B-}o3A*5Xa1MAXcmntd`u-hy z-R{Mdc4G5v}eg%E@QHgwgJImZtR?nsV zLF4@jNB}qV>nMT8Zofj!$(GVW@7~+AP^kOsn#OJD?mgSKaIchA7+!{PYsuL0f0eAb z`_ay`p<6 &, 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(); diff --git a/src/impl/particles_impl_coal.ipp b/src/impl/particles_impl_coal.ipp index 8708901b1..6e39ddddf 100644 --- a/src/impl/particles_impl_coal.ipp +++ b/src/impl/particles_impl_coal.ipp @@ -237,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 @@ -332,14 +332,6 @@ namespace libcloudphxx > > zip_rw_t; - typedef thrust::zip_iterator< - thrust::tuple< - pi_real_t, // rhod - pi_real_t, // eta - pi_real_t // tke dissipation rate - > - > zip_ro_calc_t; //read-only parameters passed to the calc() function, TODO: add Re_lambda for Onishi kernel - zip_ro_t zip_ro_it( thrust::make_tuple( // u01 @@ -358,24 +350,28 @@ namespace libcloudphxx ); auto 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()), - // tke dissipation rate - thrust::make_permutation_iterator(thrust::make_constant_iterator(0), sorted_ijk.begin()) + 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()) + ) ) ; - zip_ro_calc_t zip_ro_calc_turb_it = - 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()) + 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()) + ) ) ; @@ -400,7 +396,7 @@ namespace libcloudphxx ); - if(opts.turb_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, diff --git a/src/particles_step.ipp b/src/particles_step.ipp index e1b166818..4f5b4451f 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -339,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) From f8555ee8116b84d10f9920beee3e32ff562ea1ba Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 11 Apr 2019 14:52:11 +0200 Subject: [PATCH 37/44] opts_init switches and bindigs for turb coal --- bindings/python/lib.cpp | 2 ++ include/libcloudph++/lgrngn/opts_init.hpp | 1 + tests/python/physics/coalescence_onishi_hall.py | 8 ++++++-- tests/python/unit/col_kernels.py | 13 +++++++++++-- 4 files changed, 20 insertions(+), 4 deletions(-) diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index d12c85d4b..bbb3b7510 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -246,6 +246,7 @@ BOOST_PYTHON_MODULE(libcloudphxx) .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) @@ -273,6 +274,7 @@ BOOST_PYTHON_MODULE(libcloudphxx) .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) diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index 1c7cedb17..ed13755f9 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -145,6 +145,7 @@ namespace libcloudphxx 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), 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/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) From d68b3ec2276f7fc0b2fa06f6ba90416363946578 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 11 Apr 2019 16:06:55 +0200 Subject: [PATCH 38/44] dist_analysis - cast output to real_t, why does log(float) produce double? --- src/impl/particles_impl_dist_analysis.ipp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/impl/particles_impl_dist_analysis.ipp b/src/impl/particles_impl_dist_analysis.ipp index a9cdda184..74cafaafd 100644 --- a/src/impl/particles_impl_dist_analysis.ipp +++ b/src/impl/particles_impl_dist_analysis.ipp @@ -58,7 +58,7 @@ namespace libcloudphxx #if !defined(__NVCC__) using std::max; #endif - log_rd_min = max(log_rd_min, log(opts_init.rd_min)); // user-defined lower limit for rd + log_rd_min = max(log_rd_min, real_t(log(opts_init.rd_min))); // user-defined lower limit for rd }; template @@ -96,7 +96,7 @@ namespace libcloudphxx #if !defined(__NVCC__) using std::max; #endif - log_rd_min = max(log_rd_min, log(opts_init.rd_min)); // user-defined lower limit for rd + log_rd_min = max(log_rd_min, real_t(log(opts_init.rd_min))); // user-defined lower limit for rd } }; }; From 0f82da23be5db641b874d6afb15289e136060054 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 11 Apr 2019 16:06:55 +0200 Subject: [PATCH 39/44] dist_analysis - cast output to real_t, why does log(float) produce double? --- src/impl/particles_impl_dist_analysis.ipp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/impl/particles_impl_dist_analysis.ipp b/src/impl/particles_impl_dist_analysis.ipp index a9cdda184..74cafaafd 100644 --- a/src/impl/particles_impl_dist_analysis.ipp +++ b/src/impl/particles_impl_dist_analysis.ipp @@ -58,7 +58,7 @@ namespace libcloudphxx #if !defined(__NVCC__) using std::max; #endif - log_rd_min = max(log_rd_min, log(opts_init.rd_min)); // user-defined lower limit for rd + log_rd_min = max(log_rd_min, real_t(log(opts_init.rd_min))); // user-defined lower limit for rd }; template @@ -96,7 +96,7 @@ namespace libcloudphxx #if !defined(__NVCC__) using std::max; #endif - log_rd_min = max(log_rd_min, log(opts_init.rd_min)); // user-defined lower limit for rd + log_rd_min = max(log_rd_min, real_t(log(opts_init.rd_min))); // user-defined lower limit for rd } }; }; From e988dfa8edb5b4bb2ec728451db772e68d65c7bc Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 12 Apr 2019 14:39:17 +0200 Subject: [PATCH 40/44] fix tau_relax calculation in cells that dont contain any SDs --- src/impl/particles_impl_hskpng_turb_ss.ipp | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/src/impl/particles_impl_hskpng_turb_ss.ipp b/src/impl/particles_impl_hskpng_turb_ss.ipp index 67429b774..28b257433 100644 --- a/src/impl/particles_impl_hskpng_turb_ss.ipp +++ b/src/impl/particles_impl_hskpng_turb_ss.ipp @@ -43,15 +43,27 @@ namespace libcloudphxx template void particles_t::impl::hskpng_turb_dot_ss() { - // calc tau_relax + // calc tau_relax (only in cells that contain any SDs) moms_all(); moms_calc(rw2.begin(), real_t(1./2), false); thrust_device::vector &tau_rlx(count_mom); thrust::transform( - tau_rlx.begin(), - tau_rlx.end(), - dv.begin(), - tau_rlx.begin(), + thrust::make_permutation_iterator( + tau_rlx.begin(), + count_ijk.begin() + ), + thrust::make_permutation_iterator( + tau_rlx.begin(), + count_ijk.begin() + ) + count_n, + thrust::make_permutation_iterator( + dv.begin(), + count_ijk.begin() + ), + thrust::make_permutation_iterator( + tau_rlx.begin(), + count_ijk.begin() + ), detail::common__turbulence__tau_relax() ); From 102e1a4c92378baac90ad145e019df075c56c3a3 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 12 Apr 2019 15:48:55 +0200 Subject: [PATCH 41/44] overflow of buffers: resize instead of thrhowing errors --- src/impl/particles_impl.ipp | 3 +++ src/impl/particles_impl_bcnd.ipp | 16 ++++++++++++++-- src/impl/particles_impl_init_hskpng_npart.ipp | 8 ++++---- 3 files changed, 21 insertions(+), 6 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 91745a7fd..1c9ba77ab 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -86,6 +86,9 @@ namespace libcloudphxx 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 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_init_hskpng_npart.ipp b/src/impl/particles_impl_init_hskpng_npart.ipp index f76b37937..c6f0e811b 100644 --- a/src/impl/particles_impl_init_hskpng_npart.ipp +++ b/src/impl/particles_impl_init_hskpng_npart.ipp @@ -68,11 +68,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); } } }; From 34416a5655b0e1eaec824d015af358c0a55a14c7 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Sat, 13 Apr 2019 14:23:09 +0200 Subject: [PATCH 42/44] onishi kernel: return 0 if eps == 0 --- src/detail/kernel_onishi_nograv.hpp | 2 ++ src/detail/kernels.hpp | 1 - 2 files changed, 2 insertions(+), 1 deletion(-) 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 fa8b8bdfb..b231b2f9a 100644 --- a/src/detail/kernels.hpp +++ b/src/detail/kernels.hpp @@ -223,7 +223,6 @@ namespace libcloudphxx common::moist_air::rho_w() / si::kilograms * si::cubic_metres / thrust::get(tpl_wrap.get_ro_calc()) // ratio of water to air density ); - printf("onishi calc diss rate: %g\n", thrust::get(tpl_wrap.get_ro_calc())); real_t res = kernel_geometric::interpolated_efficiency(rwa, rwb) * // stagnant air collision efficiency wang_collision_enhancement(rwa, rwb, kernel_base::k_params[0]) * // Wang turbulent collision efficiency enhancement, k_params[0] - epsilon From 26605c7f81326c8f7299fa0884f770ece7f7b111 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 24 Apr 2019 12:57:09 +0200 Subject: [PATCH 43/44] fix tau_relax calculation in cells that dont contain any SDs - cd --- src/impl/particles_impl_hskpng_turb_ss.ipp | 33 +++++++++++++--------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/src/impl/particles_impl_hskpng_turb_ss.ipp b/src/impl/particles_impl_hskpng_turb_ss.ipp index 28b257433..429402088 100644 --- a/src/impl/particles_impl_hskpng_turb_ss.ipp +++ b/src/impl/particles_impl_hskpng_turb_ss.ipp @@ -30,6 +30,7 @@ namespace libcloudphxx 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, @@ -43,28 +44,34 @@ namespace libcloudphxx template void particles_t::impl::hskpng_turb_dot_ss() { - // calc tau_relax (only in cells that contain any SDs) + 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_device::vector &tau_rlx(count_mom); thrust::transform( - thrust::make_permutation_iterator( - tau_rlx.begin(), - count_ijk.begin() - ), - thrust::make_permutation_iterator( - tau_rlx.begin(), - count_ijk.begin() - ) + count_n, + 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(), + tau_rlx.begin(), count_ijk.begin() - ), - detail::common__turbulence__tau_relax() + ) ); // calc ss perturb for each SD From 4ed8660c189f47cb77b2f8d12b30975d473a9f70 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 13 May 2019 15:08:12 +0200 Subject: [PATCH 44/44] SGS length scale = dz --- include/libcloudph++/common/turbulence.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/include/libcloudph++/common/turbulence.hpp b/include/libcloudph++/common/turbulence.hpp index f7bd26118..3ac1f7748 100644 --- a/include/libcloudph++/common/turbulence.hpp +++ b/include/libcloudph++/common/turbulence.hpp @@ -66,7 +66,8 @@ namespace libcloudphxx quantity dx, quantity dz ) - { return quantity(sqrt(dx*dz));} + //{ return quantity(sqrt(dx*dz));} + { return quantity(dz);} template quantity length_scale( @@ -74,7 +75,8 @@ namespace libcloudphxx quantity dy, quantity dz ) - { return quantity(pow((dx/si::metres)*(dy/si::metres)*(dz/si::metres), real_t(1./3.)) * si::metres);} + //{ 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