diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index c012d8823..08f120d88 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -104,6 +104,9 @@ namespace libcloudphxx // do we want to track the time SDs spend inside clouds bool diag_incloud_time; + // do we want to calculate mass flux caused by mass teleportation in the coalescence algorithm + bool diag_coal_tele_mass_flux; + // RH threshold for calculating equilibrium condition at t=0 real_t RH_max; @@ -232,6 +235,7 @@ namespace libcloudphxx supstp_rlx(1), rd_min(0.), diag_incloud_time(false), + diag_coal_tele_mass_flux(true), no_ccn_at_init(false), open_side_walls(false), periodic_topbot_walls(false), diff --git a/include/libcloudph++/lgrngn/particles.hpp b/include/libcloudph++/lgrngn/particles.hpp index f494b8d9a..51774b5ce 100644 --- a/include/libcloudph++/lgrngn/particles.hpp +++ b/include/libcloudph++/lgrngn/particles.hpp @@ -108,6 +108,7 @@ namespace libcloudphxx virtual void diag_incloud_time_mom(const int&) { assert(false); } // requires opts_init.diag_incloud_time==true virtual void diag_max_rw() { assert(false); } virtual void diag_vel_div() { assert(false); } + virtual void diag_coal_tele_mass_flux() { assert(false); } virtual std::map diag_puddle() { assert(false); return std::map(); } virtual real_t *outbuf() { assert(false); return NULL; } @@ -197,6 +198,7 @@ namespace libcloudphxx void diag_precip_rate(); void diag_max_rw(); void diag_vel_div(); + void diag_coal_tele_mass_flux(); std::map diag_puddle(); real_t *outbuf(); diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index b07d3332c..d764f3d1e 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -147,7 +147,8 @@ namespace libcloudphxx p, // pressure [Pa] RH, // relative humisity eta,// dynamic viscosity - diss_rate; // turbulent kinetic energy dissipation rate + diss_rate, // turbulent kinetic energy dissipation rate + coal_tele_mass_flux; // vertical mass flux due to mass "teleportation" at collision (diagnostic, accumulated over time) thrust_device::vector w_LS; // large-scale subsidence velocity profile thrust_device::vector SGS_mix_len; // SGS mixing length profile diff --git a/src/impl/particles_impl_coal.ipp b/src/impl/particles_impl_coal.ipp index 205ccacb2..ccd1d96ec 100644 --- a/src/impl/particles_impl_coal.ipp +++ b/src/impl/particles_impl_coal.ipp @@ -102,15 +102,22 @@ namespace libcloudphxx typename tup_t > BOOST_GPU_ENABLED - void collide(tup_t tpl, const n_t &col_no) + void collide(tup_t tpl, const n_t &col_no, const real_t &z_a, const real_t &z_b, real_t &coal_tele_mass_flux_pp) { #if !defined(__NVCC__) using std::cbrt; using std::sqrt; #endif + // multiplicity change (eq. 12 in Shima et al. 2009) thrust::get(tpl) -= col_no * thrust::get(tpl); + // calculate vertical mass flux from _a to _b + // teleported mass = col_no * n_b * mass_a + // distance = z_b - z_a, >0 mean upward flux + coal_tele_mass_flux_pp = (z_b - z_a) * col_no * thrust::get(tpl) * thrust::get(tpl) * sqrt(thrust::get(tpl)); + // done later: * 4/3 PI rho_w / dt + // wet radius change (eq. 13 in Shima et al. 2009) const real_t rw_b = cbrt( col_no * thrust::get(tpl) * sqrt(thrust::get(tpl)) + @@ -132,15 +139,17 @@ namespace libcloudphxx template struct collider { - // read-only parameters + // read-only parameters (not really, coal_tele_mass_flux_pp is overwritten) typedef thrust::tuple< real_t, // random number (u01) real_t, // scaling factor thrust_size_t, thrust_size_t, // ix thrust_size_t, thrust_size_t, // off (index within cell) - real_t // dv + real_t, // dv + real_t, real_t, // z + real_t // coal_tele_mass_flux_pp > tpl_ro_t; - enum { u01_ix, scl_ix, ix_a_ix, ix_b_ix, off_a_ix, off_b_ix, dv_ix }; + enum { u01_ix, scl_ix, ix_a_ix, ix_b_ix, off_a_ix, off_b_ix, dv_ix, z_a_ix, z_b_ix, coal_tele_mass_flux_pp_ix }; // read-write parameters = return type typedef thrust::tuple< @@ -236,7 +245,7 @@ namespace libcloudphxx rw2_a_ix, rw2_b_ix, rd3_a_ix, rd3_b_ix, vt_a_ix, vt_b_ix - >(thrust::get<1>(tpl_ro_rw), col_no); + >(thrust::get<1>(tpl_ro_rw), col_no, thrust::get(thrust::get<0>(tpl_ro_rw)), thrust::get(thrust::get<0>(tpl_ro_rw)), thrust::get(thrust::get<0>(tpl_ro_rw))); thrust::get(thrust::get<1>(tpl_ro_rw)) = real_t(na_ge_nb); // col vector for the second in a pair stores info on which one has greater multiplicity } else @@ -248,7 +257,7 @@ namespace libcloudphxx rw2_b_ix, rw2_a_ix, rd3_b_ix, rd3_a_ix, vt_b_ix, vt_a_ix - >(thrust::get<1>(tpl_ro_rw), col_no); + >(thrust::get<1>(tpl_ro_rw), col_no, thrust::get(thrust::get<0>(tpl_ro_rw)), thrust::get(thrust::get<0>(tpl_ro_rw)), thrust::get(thrust::get<0>(tpl_ro_rw))); thrust::get(thrust::get<1>(tpl_ro_rw)) = real_t(nb_gt_na); // col vector for the second in a pair stores info on which one has greater multiplicity } thrust::get(thrust::get<1>(tpl_ro_rw)) = real_t(col_no); // col vector for the first in a pair stores info on number of collisions @@ -274,14 +283,19 @@ namespace libcloudphxx // references to tmp data thrust_device::vector &scl(tmp_device_real_cell), // scale factor for probablility - &col(tmp_device_real_part); // number of collisions, used in chemistry, NOTE: it's the same as u01, so it overwrites already used random numbers + &col(tmp_device_real_part), // number of collisions, used in chemistry, NOTE: it's the same as u01, so it overwrites already used random numbers // 1st one of a pair stores number of collisions, 2nd one stores info on which one has greater multiplicity + &coal_tele_mass_flux_pp(tmp_device_real_part2); // mass flux caused by teleportation of droplets at coalescence + // 1st one of a pair - data, 2nd - empty thrust_device::vector &off(tmp_device_size_cell); // offset for getting index of particle within a cell // laying out scale factor onto ijk grid // fill with 0s if not all cells will be updated in the following copy if(count_n!=n_cell) thrust::fill(scl.begin(), scl.end(), real_t(0.)); + + if(opts_init.diag_coal_tele_mass_flux) + thrust::fill(coal_tele_mass_flux_pp.begin(), coal_tele_mass_flux_pp.end(), real_t(0)); thrust::copy( count_mom.begin(), // input - begin @@ -338,7 +352,9 @@ namespace libcloudphxx thrust::counting_iterator, // ix_a thrust::counting_iterator, // ix_b pi_size_t, pi_size_t, // off_a & off_b - pi_real_t // dv + pi_real_t, // dv + pi_real_t, pi_real_t, // z_a & z_b + i_real_t // col_tele_mass_flux > > zip_ro_t; @@ -365,7 +381,13 @@ namespace libcloudphxx thrust::make_permutation_iterator(off.begin(), sorted_ijk.begin()), thrust::make_permutation_iterator(off.begin(), sorted_ijk.begin())+1, // dv - thrust::make_permutation_iterator(dv.begin(), sorted_ijk.begin()) + thrust::make_permutation_iterator(dv.begin(), sorted_ijk.begin()), + // vertical position + thrust::make_permutation_iterator(z.begin(), sorted_id.begin()), + thrust::make_permutation_iterator(z.begin(), sorted_id.begin())+1, + // mass flux caused by teleportation in the coalescence algorithm - its overwritten, not read-only! + // TODO: pass it only if opts_init.diag_coal_tele_mass_flux==true + coal_tele_mass_flux_pp.begin() ) ); @@ -503,6 +525,37 @@ namespace libcloudphxx ); nancheck(incloud_time, "incloud_time - post coalescence"); } + + // per-cell sum of coalescence teleportation mass flux + if(opts_init.diag_coal_tele_mass_flux) + { + thrust::fill(tmp_device_real_cell.begin(), tmp_device_real_cell.end(), 0); // TODO: not needed? + + // store sum from this step in tmp_device_real_cell + auto it_pair = thrust::reduce_by_key( + // input - keys + sorted_ijk.begin(), sorted_ijk.begin()+n_part, + // input - values + coal_tele_mass_flux_pp.begin(), + // output - keys + count_ijk.begin(), + // output - values + tmp_device_real_cell.begin() + ); + + count_n = it_pair.first - count_ijk.begin(); + assert(count_n <= n_cell); + + // add to accumulated values from previous steps +// thrust::transform(tmp_device_real_cell.begin(), tmp_device_real_cell.end(), coal_tele_mass_flux.begin(), coal_tele_mass_flux.begin(), thrust::plus()); + thrust::transform( + tmp_device_real_cell.begin(), + tmp_device_real_cell.begin() + count_n, + thrust::make_permutation_iterator(coal_tele_mass_flux.begin(), count_ijk.begin()), + thrust::make_permutation_iterator(coal_tele_mass_flux.begin(), count_ijk.begin()), + thrust::plus() + ); + } } }; }; diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index d9f494bb7..f3a561c98 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -53,7 +53,7 @@ namespace libcloudphxx { tmp_device_real_part1.resize(n_part); } - if((allow_sstp_cond && opts_init.exact_sstp_cond) || n_dims==3 || opts_init.turb_cond_switch) + if((allow_sstp_cond && opts_init.exact_sstp_cond) || n_dims==3 || opts_init.turb_cond_switch || opts_init.diag_coal_tele_mass_flux) { tmp_device_real_part2.resize(n_part); } diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index 2b92d6601..8be126c84 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -35,6 +35,8 @@ namespace libcloudphxx sstp_tmp_th.resize(n_cell); sstp_tmp_rh.resize(n_cell); } + if(opts_init.diag_coal_tele_mass_flux) + coal_tele_mass_flux.resize(n_cell, 0); } }; }; diff --git a/src/impl/particles_impl_reserve_hskpng_npart.ipp b/src/impl/particles_impl_reserve_hskpng_npart.ipp index 063969e7f..5ab28e6b1 100644 --- a/src/impl/particles_impl_reserve_hskpng_npart.ipp +++ b/src/impl/particles_impl_reserve_hskpng_npart.ipp @@ -56,7 +56,7 @@ namespace libcloudphxx { tmp_device_real_part1.reserve(opts_init.n_sd_max); } - if((allow_sstp_cond && opts_init.exact_sstp_cond) || n_dims==3 || opts_init.turb_cond_switch) + if((allow_sstp_cond && opts_init.exact_sstp_cond) || n_dims==3 || opts_init.turb_cond_switch || opts_init.diag_coal_tele_mass_flux) { tmp_device_real_part2.reserve(opts_init.n_sd_max); } diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index a948143a3..790e84869 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -349,6 +349,18 @@ namespace libcloudphxx assert(0 && "diag_incloud_time_mom called, but opts_init.diag_incloud_time==false"); } + // output accumulated mass flux due to coalescence teleportation + template + void particles_t::diag_coal_tele_mass_flux() + { + assert(pimpl->opts_init.diag_coal_tele_mass_flux && "diag_coal_tele_mass_flux called, but opts_init.diag_coal_tele_mass_flux==false"); + thrust::copy(pimpl->coal_tele_mass_flux.begin(), pimpl->coal_tele_mass_flux.end(), pimpl->count_mom.begin()); + + // mass flux defined in all cells + pimpl->count_n = pimpl->n_cell; + thrust::sequence(pimpl->count_ijk.begin(), pimpl->count_ijk.end()); + } + // computes mass density function for wet radii using estimator from Shima et al. (2009) template void particles_t::diag_wet_mass_dens(const real_t &rad, const real_t &sig0)