From 814c7f8de66608479bdf87b5bc673120909546d9 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 5 May 2022 14:57:56 +0200 Subject: [PATCH 01/38] grid refinement: ijk_ref and its values + syncing in refined arrays --- include/libcloudph++/lgrngn/opts_init.hpp | 4 +- src/impl/particles_impl.ipp | 9 ++- src/impl/particles_impl_hskpng_ijk.ipp | 63 ++++++++++++++----- src/impl/particles_impl_hskpng_remove.ipp | 1 + src/impl/particles_impl_hskpng_resize.ipp | 1 + .../particles_impl_init_SD_with_distros.ipp | 17 ++--- .../particles_impl_init_SD_with_sizes.ipp | 15 +++-- src/impl/particles_impl_init_e2l.ipp | 16 +++-- src/impl/particles_impl_init_sanity_check.ipp | 11 ++++ .../particles_impl_reserve_hskpng_npart.ipp | 1 + src/impl/particles_impl_rlx_dry_distros.ipp | 1 + ...articles_impl_src_dry_distros_matching.ipp | 8 ++- .../particles_impl_src_dry_distros_simple.ipp | 10 +-- src/particles_init.ipp | 15 ++--- src/particles_step.ipp | 11 +--- 15 files changed, 122 insertions(+), 61 deletions(-) diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index c012d8823..cfdab7a7c 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -46,6 +46,7 @@ namespace libcloudphxx // Eulerian component parameters int nx, ny, nz; real_t dx, dy, dz, dt; + unsigned int n_ref; // number of grid refinements per 'normal' cell; refined cell size is (nx, ny, nz) / n_ref; tht and rv are on refined cells, rest on normal cells; condensation is resolved on refined cells // no. of substeps int sstp_cond, sstp_coal; @@ -236,7 +237,8 @@ namespace libcloudphxx open_side_walls(false), periodic_topbot_walls(false), variable_dt_switch(false), - rng_seed_init_switch(false) + rng_seed_init_switch(false), + n_ref(1) {} // dtor (just to silence -Winline warnings) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index b07d3332c..e809b9ef9 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -46,6 +46,7 @@ namespace libcloudphxx opts_init_t opts_init; // a copy const int n_dims; const thrust_size_t n_cell; + const int ny_ref, nz_ref; // number of refined cells thrust_size_t n_part, // total number of SDs n_part_old, // total number of SDs before source n_part_to_init; // number of SDs to be initialized by source @@ -107,6 +108,7 @@ namespace libcloudphxx // housekeeping data (per particle) thrust_device::vector i, j, k, ijk, // Eulerian grid cell indices (always zero for 0D) + ijk_ref, // ijk in refined cells sorted_id, sorted_ijk; // Arakawa-C grid helper vars @@ -333,6 +335,8 @@ namespace libcloudphxx n_dims == 2 ? halo_size * (_opts_init.nz + 1): // 2D halo_size * (_opts_init.nz + 1) * _opts_init.ny // 3D ), + ny_ref(n_dims == 3 ? (opts_init.ny - 1) * opts_init.n_ref + 1 : 0), + nz_ref(n_dims >= 2 ? (opts_init.nz - 1) * opts_init.n_ref + 1 : 0), w_LS(_opts_init.w_LS), SGS_mix_len(_opts_init.SGS_mix_len), adve_scheme(_opts_init.adve_scheme), @@ -467,7 +471,7 @@ namespace libcloudphxx void init_count_num_src(const thrust_size_t &); template void conc_to_number(arr_t &arr); - void init_e2l(const arrinfo_t &, thrust_device::vector*, const int = 0, const int = 0, const int = 0, const long int = 0); + void init_e2l(const arrinfo_t &, thrust_device::vector*, const bool = false, const int = 0, const int = 0, const int = 0, const long int = 0); void init_wet(); void init_sync(); void init_grid(); @@ -487,9 +491,10 @@ namespace libcloudphxx void hskpng_sort(); void hskpng_shuffle_and_sort(); void hskpng_count(); - void ravel_ijk(const thrust_size_t begin_shift = 0); + void ravel_ijk(const thrust_size_t begin_shift = 0, const bool refined = false); void unravel_ijk(const thrust_size_t begin_shift = 0); void hskpng_ijk(); + void hskpng_ijk_ref(const thrust_size_t begin_shift = 0, const bool unravel_ijk = true); void hskpng_Tpr(); void hskpng_mfp(); diff --git a/src/impl/particles_impl_hskpng_ijk.ipp b/src/impl/particles_impl_hskpng_ijk.ipp index 2320462f2..16a99da78 100644 --- a/src/impl/particles_impl_hskpng_ijk.ipp +++ b/src/impl/particles_impl_hskpng_ijk.ipp @@ -30,21 +30,22 @@ namespace libcloudphxx // calc ijk from i, j and k template - void particles_t::impl::ravel_ijk(const thrust_size_t begin_shift) // default = 0 + void particles_t::impl::ravel_ijk(const thrust_size_t begin_shift, const bool refined) { + auto & _ijk = refined ? ijk_ref : ijk; switch (n_dims) { case 0: break; case 1: - thrust::copy(i.begin()+begin_shift, i.end(), ijk.begin()+begin_shift); + thrust::copy(i.begin()+begin_shift, i.end(), _ijk.begin()+begin_shift); break; case 2: namespace arg = thrust::placeholders; thrust::transform( i.begin()+begin_shift, i.end(), // input - first arg k.begin()+begin_shift, // input - second arg - ijk.begin()+begin_shift, // output + _ijk.begin()+begin_shift, // output arg::_1 * opts_init.nz + arg::_2 // assuming z varies first ); break; @@ -53,14 +54,14 @@ namespace libcloudphxx thrust::transform( i.begin()+begin_shift, i.end(), // input - first arg j.begin()+begin_shift, // input - second arg - ijk.begin()+begin_shift, // output + _ijk.begin()+begin_shift, // output arg::_1 * (opts_init.nz * opts_init.ny) + arg::_2 * opts_init.nz ); thrust::transform( - ijk.begin()+begin_shift, ijk.end(), + _ijk.begin()+begin_shift, _ijk.end(), k.begin()+begin_shift, - ijk.begin()+begin_shift, // in-place! + _ijk.begin()+begin_shift, // in-place! arg::_1 + arg::_2 ); // TODO: replace these two transforms with single one @@ -122,23 +123,57 @@ namespace libcloudphxx } } - template - void particles_t::impl::hskpng_ijk() - { - // helper functor - struct { + namespace { + // helper functor to calculate cell index + template + struct cell_idx_hlpr_t{ + + const thrust_size_t begin_shift; + void operator()( thrust_device::vector &vx, thrust_device::vector &vi, const real_t &vd ) { thrust::transform( - vx.begin(), vx.end(), // input - vi.begin(), // output + vx.begin() + begin_shift, vx.end(), // input + vi.begin() + begin_shift, // output detail::divide_by_constant_and_cast(vd) // has to be done on doubles to avoid i==nx due to low precision of nvcc math; TODO: now that rand uniform has range [0,1), float might be good here? ); } - } helper; + + cell_idx_hlpr_t(const thrust_size_t begin_shift = 0): + begin_shift(begin_shift) + {} + }; + }; + + // ijk_ref - cell number in the refined grid + template + void particles_t::impl::hskpng_ijk_ref(const thrust_size_t begin_shift, const bool _unravel_ijk) + { + if(opts_init.n_ref <= 1) return; + + cell_idx_hlpr_t helper(begin_shift); + + if (opts_init.nx != 0) helper(x, i, opts_init.dx / opts_init.n_ref); + if (opts_init.ny != 0) helper(y, j, opts_init.dy / opts_init.n_ref); + if (opts_init.nz != 0) helper(z, k, opts_init.dz / opts_init.n_ref); + + // raveling i, j & k into ijk_ref + ravel_ijk(begin_shift, true); + + // go back to i, j, k in the normal cell range + if(_unravel_ijk) + unravel_ijk(begin_shift); + } + + template + void particles_t::impl::hskpng_ijk() + { + hskpng_ijk_ref(0, false); + + cell_idx_hlpr_t helper; if (opts_init.nx != 0) helper(x, i, opts_init.dx); if (opts_init.ny != 0) helper(y, j, opts_init.dy); diff --git a/src/impl/particles_impl_hskpng_remove.ipp b/src/impl/particles_impl_hskpng_remove.ipp index ae31eee6b..1b783bf71 100644 --- a/src/impl/particles_impl_hskpng_remove.ipp +++ b/src/impl/particles_impl_hskpng_remove.ipp @@ -26,6 +26,7 @@ namespace libcloudphxx 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.n_ref > 1) n_t_vctrs.insert(&ijk_ref); namespace arg = thrust::placeholders; diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index d9f494bb7..f8dc28a56 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -24,6 +24,7 @@ namespace libcloudphxx n.resize(n_part); tmp_device_n_part.resize(n_part); tmp_device_size_part.resize(n_part); + if (opts_init.n_ref > 1) ijk_ref.resize(n_part); vt.resize(n_part, detail::invalid); diff --git a/src/impl/particles_impl_init_SD_with_distros.ipp b/src/impl/particles_impl_init_SD_with_distros.ipp index 0a20ba6c1..d6ad1275f 100644 --- a/src/impl/particles_impl_init_SD_with_distros.ipp +++ b/src/impl/particles_impl_init_SD_with_distros.ipp @@ -52,6 +52,16 @@ namespace libcloudphxx template void particles_t::impl::init_SD_with_distros_finalize(const real_t &kappa, const bool unravel_ijk_switch) { + // ijk -> i, j, k + if(unravel_ijk_switch) + unravel_ijk(n_part_old); + + // initialising particle positions + init_xyz(); + + // ijk of the refined cells + hskpng_ijk_ref(n_part_old); + // init kappa init_kappa(kappa); @@ -77,13 +87,6 @@ namespace libcloudphxx if (opts_init.chem_switch){ chem_vol_ante(); } - - // ijk -> i, j, k - if(unravel_ijk_switch) - unravel_ijk(n_part_old); - - // initialising particle positions - init_xyz(); if(opts_init.diag_incloud_time) init_incloud_time(); diff --git a/src/impl/particles_impl_init_SD_with_sizes.ipp b/src/impl/particles_impl_init_SD_with_sizes.ipp index b7e89cccb..f0e3563a0 100644 --- a/src/impl/particles_impl_init_SD_with_sizes.ipp +++ b/src/impl/particles_impl_init_SD_with_sizes.ipp @@ -41,6 +41,15 @@ namespace libcloudphxx // init ijk vector using count_num, also n_part and resize n_part vectors init_ijk(); + + // ijk -> i, j, k + unravel_ijk(n_part_old); + + // initialising particle positions + init_xyz(); + + // ijk of the refined cells + hskpng_ijk_ref(n_part_old); // initialising dry radii (needs ijk) init_dry_dry_sizes(sni->first); @@ -73,12 +82,6 @@ namespace libcloudphxx if (opts_init.chem_switch){ chem_vol_ante(); } - - // ijk -> i, j, k - unravel_ijk(n_part_old); - - // initialising particle positions - init_xyz(); } } } diff --git a/src/impl/particles_impl_init_e2l.ipp b/src/impl/particles_impl_init_e2l.ipp index a2fd3d13e..1d658a7e9 100644 --- a/src/impl/particles_impl_init_e2l.ipp +++ b/src/impl/particles_impl_init_e2l.ipp @@ -34,6 +34,7 @@ namespace libcloudphxx void particles_t::impl::init_e2l( const arrinfo_t &arr, thrust_device::vector * key, + const bool refined, const int ext_x, const int ext_y, const int ext_z, const long int offset ) @@ -49,6 +50,9 @@ namespace libcloudphxx long int max_stride = 0; // stride of the dimension with highest stride long int max_stride_n_cell = 0; // number of cells in that direction in this process (with halo) + const auto ny = refined ? ny_ref : opts_init.ny, + nz = refined ? nz_ref : opts_init.nz; + switch (n_dims) { namespace arg = thrust::placeholders; @@ -79,8 +83,8 @@ namespace libcloudphxx // output l2e[key].begin(), // op - arr.strides[0] * /* i = */ (arg::_1 / (opts_init.nz + ext_z)) + - arr.strides[1] * /* j = */ (arg::_1 % (opts_init.nz + ext_z)) // module of negative value might not work in 2003 standard? + arr.strides[0] * /* i = */ (arg::_1 / (nz + ext_z)) + + arr.strides[1] * /* j = */ (arg::_1 % (nz + ext_z)) // module of negative value might not work in 2003 standard? ); max_stride = arr.strides[0]; max_stride_n_cell = n_x_tot + ext_x; @@ -94,12 +98,12 @@ namespace libcloudphxx // output l2e[key].begin(), // op - arr.strides[0] * /* i = */ (arg::_1 / ((opts_init.nz + ext_z) * (opts_init.ny + ext_y))) + - arr.strides[1] * /* j = */ ((arg::_1 / (opts_init.nz + ext_z)) % (opts_init.ny + ext_y)) + - arr.strides[2] * /* k = */ (arg::_1 % ((opts_init.nz + ext_z))) + arr.strides[0] * /* i = */ (arg::_1 / ((nz + ext_z) * (ny + ext_y))) + + arr.strides[1] * /* j = */ ((arg::_1 / (nz + ext_z)) % (ny + ext_y)) + + arr.strides[2] * /* k = */ (arg::_1 % ((nz + ext_z))) ); max_stride = std::max(arr.strides[0], arr.strides[1]); // handles kji (C-style) and kij storage orders - max_stride_n_cell = max_stride == arr.strides[0] ? n_x_tot + ext_x : opts_init.ny + ext_y; + max_stride_n_cell = max_stride == arr.strides[0] ? n_x_tot + ext_x : ny + ext_y; break; default: assert(false); } diff --git a/src/impl/particles_impl_init_sanity_check.ipp b/src/impl/particles_impl_init_sanity_check.ipp index 5d309dafa..c4b763f28 100644 --- a/src/impl/particles_impl_init_sanity_check.ipp +++ b/src/impl/particles_impl_init_sanity_check.ipp @@ -141,6 +141,17 @@ namespace libcloudphxx throw std::runtime_error("rlx_timescale <= 0"); if(opts_init.rlx_switch && opts_init.chem_switch) throw std::runtime_error("CCN relaxation does not work with chemistry"); + +/* + if(n_dims == 3 && opts_init.n_ref > 1 && (opts_init.ny - 1) % opts_init.n_ref != 0) + throw std::runtime_error("opts_init.ny is not correctly divisible by n_ref ((opts_init.ny - 1) % opts_init.n_ref != 0)"); + if(n_dims >= 2 && opts_init.n_ref > 1 && (opts_init.nz - 1) % opts_init.n_ref != 0) + throw std::runtime_error("opts_init.nz is not correctly divisible by n_ref ((opts_init.nz - 1) % opts_init.n_ref != 0)"); + */ + if(n_dims <= 1 && opts_init.n_ref > 1) + throw std::runtime_error("opts_init.n_ref > 1 works only in 2D and 3D"); + if(opts_init.n_ref == 0) + throw std::runtime_error("opts_init.n_ref==0; set opts_init.n_ref=1 for no refinement"); } }; }; diff --git a/src/impl/particles_impl_reserve_hskpng_npart.ipp b/src/impl/particles_impl_reserve_hskpng_npart.ipp index 063969e7f..bbe272b65 100644 --- a/src/impl/particles_impl_reserve_hskpng_npart.ipp +++ b/src/impl/particles_impl_reserve_hskpng_npart.ipp @@ -17,6 +17,7 @@ namespace libcloudphxx if (opts_init.ny != 0) j.reserve(opts_init.n_sd_max); // > TODO: are they needed at all? if (opts_init.nz != 0) k.reserve(opts_init.n_sd_max); // ijk.reserve(opts_init.n_sd_max); + if(opts_init.n_ref > 1) ijk_ref.reserve(n_part); if (opts_init.nx != 0) x.reserve(opts_init.n_sd_max); if (opts_init.ny != 0) y.reserve(opts_init.n_sd_max); diff --git a/src/impl/particles_impl_rlx_dry_distros.ipp b/src/impl/particles_impl_rlx_dry_distros.ipp index 528b9eca7..e5a0de7db 100644 --- a/src/impl/particles_impl_rlx_dry_distros.ipp +++ b/src/impl/particles_impl_rlx_dry_distros.ipp @@ -222,6 +222,7 @@ namespace libcloudphxx if(n_dims==3) j.resize(n_part); // we dont check in i and k because relax works only in 2D and 3D rd3.resize(n_part); n.resize(n_part); + if(opts_init.n_ref > 1) ijk_ref.resize(n_part); // --- init k --- thrust_device::vector &ptr(tmp_device_size_cell); diff --git a/src/impl/particles_impl_src_dry_distros_matching.ipp b/src/impl/particles_impl_src_dry_distros_matching.ipp index 61faa20bb..f0c201fd9 100644 --- a/src/impl/particles_impl_src_dry_distros_matching.ipp +++ b/src/impl/particles_impl_src_dry_distros_matching.ipp @@ -218,9 +218,6 @@ namespace libcloudphxx init_n_sd_conc( *(opts_init.src_dry_distros.begin()->second) ); // TODO: document that n_of_lnrd_stp is expected! - - // init rw - init_wet(); // ijk -> i, j, k unravel_ijk(n_part_old); @@ -228,6 +225,11 @@ namespace libcloudphxx // init x, y, z, i, j, k init_xyz(); + hskpng_ijk_ref(n_part_old); + + // init rw + init_wet(); + // TODO: init chem { diff --git a/src/impl/particles_impl_src_dry_distros_simple.ipp b/src/impl/particles_impl_src_dry_distros_simple.ipp index 8b7d733b8..df90acfd6 100644 --- a/src/impl/particles_impl_src_dry_distros_simple.ipp +++ b/src/impl/particles_impl_src_dry_distros_simple.ipp @@ -53,16 +53,18 @@ namespace libcloudphxx init_n_sd_conc( *(opts_init.src_dry_distros.begin()->second) ); // TODO: document that n_of_lnrd_stp is expected! - - // init rw - init_wet(); // ijk -> i, j, k unravel_ijk(n_part_old); - // init x, y, z, i, j, k + // init x, y, z init_xyz(); + hskpng_ijk_ref(n_part_old); + + // init rw + init_wet(); + // TODO: init chem } }; diff --git a/src/particles_init.ipp b/src/particles_init.ipp index 5fcf9554c..9f95e0c55 100644 --- a/src/particles_init.ipp +++ b/src/particles_init.ipp @@ -39,8 +39,8 @@ namespace libcloudphxx // initialising Eulerian-Lagrangian coupling pimpl->init_sync(); // also, init of ambient_chem vectors - pimpl->init_e2l(th, &pimpl->th); - pimpl->init_e2l(rv, &pimpl->rv); + pimpl->init_e2l(th, &pimpl->th, true); + pimpl->init_e2l(rv, &pimpl->rv, true); pimpl->init_e2l(rhod, &pimpl->rhod); if(pimpl->const_p) pimpl->init_e2l(p, &pimpl->p); @@ -48,14 +48,9 @@ namespace libcloudphxx #if !defined(__NVCC__) using std::max; #endif - if (!courant_x.is_null()) pimpl->init_e2l(courant_x, &pimpl->courant_x, 1, 0, 0, - pimpl->halo_x); - if (!courant_y.is_null()) pimpl->init_e2l(courant_y, &pimpl->courant_y, 0, 1, 0, pimpl->n_x_bfr * pimpl->opts_init.nz - pimpl->halo_y); - 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 (!courant_x.is_null()) pimpl->init_e2l(courant_x, &pimpl->courant_x, 1, 0, 0); - if (!courant_y.is_null()) pimpl->init_e2l(courant_y, &pimpl->courant_y, 0, 1, 0, pimpl->n_x_bfr * pimpl->opts_init.nz ); - 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) ); - */ + if (!courant_x.is_null()) pimpl->init_e2l(courant_x, &pimpl->courant_x, false, 1, 0, 0, - pimpl->halo_x); + if (!courant_y.is_null()) pimpl->init_e2l(courant_y, &pimpl->courant_y, false, 0, 1, 0, pimpl->n_x_bfr * pimpl->opts_init.nz - pimpl->halo_y); + if (!courant_z.is_null()) pimpl->init_e2l(courant_z, &pimpl->courant_z, false, 0, 0, 1, pimpl->n_x_bfr * max(1, pimpl->opts_init.ny) - pimpl->halo_z); if (pimpl->opts_init.chem_switch) for (int i = 0; i < chem_gas_n; ++i) diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 0ffee2b0e..0bf108f47 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -81,14 +81,9 @@ namespace libcloudphxx if (pimpl->l2e[&pimpl->courant_x].size() == 0) // TODO: y, z,... { // TODO: copy-pasted from init - if (!courant_x.is_null()) pimpl->init_e2l(courant_x, &pimpl->courant_x, 1, 0, 0, - pimpl->halo_x); - if (!courant_y.is_null()) pimpl->init_e2l(courant_y, &pimpl->courant_y, 0, 1, 0, pimpl->n_x_bfr * pimpl->opts_init.nz - pimpl->halo_y); - if (!courant_z.is_null()) pimpl->init_e2l(courant_z, &pimpl->courant_z, 0, 0, 1, pimpl->n_x_bfr * std::max(1, pimpl->opts_init.ny) - pimpl->halo_z); - /* - if (!courant_x.is_null()) pimpl->init_e2l(courant_x, &pimpl->courant_x, 1, 0, 0); - if (!courant_y.is_null()) pimpl->init_e2l(courant_y, &pimpl->courant_y, 0, 1, 0, pimpl->n_x_bfr * pimpl->opts_init.nz ); - if (!courant_z.is_null()) pimpl->init_e2l(courant_z, &pimpl->courant_z, 0, 0, 1, pimpl->n_x_bfr * std::max(1, pimpl->opts_init.ny) ); - */ + if (!courant_x.is_null()) pimpl->init_e2l(courant_x, &pimpl->courant_x, false, 1, 0, 0, - pimpl->halo_x); + if (!courant_y.is_null()) pimpl->init_e2l(courant_y, &pimpl->courant_y, false, 0, 1, 0, pimpl->n_x_bfr * pimpl->opts_init.nz - pimpl->halo_y); + if (!courant_z.is_null()) pimpl->init_e2l(courant_z, &pimpl->courant_z, false, 0, 0, 1, pimpl->n_x_bfr * std::max(1, pimpl->opts_init.ny) - pimpl->halo_z); } if (pimpl->l2e[&pimpl->diss_rate].size() == 0) From daf9c919227a1ae60ce105ac9e47ec183e6f6cac Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 20 Jun 2022 15:18:33 +0200 Subject: [PATCH 02/38] hskpng ijk_ref: assume MPI boundary cuts a refined cell in half --- src/impl/particles_impl_hskpng_ijk.ipp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/impl/particles_impl_hskpng_ijk.ipp b/src/impl/particles_impl_hskpng_ijk.ipp index 16a99da78..d9184239b 100644 --- a/src/impl/particles_impl_hskpng_ijk.ipp +++ b/src/impl/particles_impl_hskpng_ijk.ipp @@ -15,15 +15,15 @@ namespace libcloudphxx { /// @brief returns ret_t(x/c) template - struct divide_by_constant_and_cast + struct add_and_divide_by_constant_and_cast { - arg_t c; - divide_by_constant_and_cast(arg_t c) : c(c) {} + arg_t c, s; + add_and_divide_by_constant_and_cast(arg_t c, arg_t s) : c(c), s(s) {} BOOST_GPU_ENABLED ret_t operator()(arg_t x) { - return ret_t(x/c); + return ret_t((x+s)/c); } }; }; @@ -133,12 +133,13 @@ namespace libcloudphxx void operator()( thrust_device::vector &vx, thrust_device::vector &vi, - const real_t &vd + const real_t &vd, + const real_t pos_shift = 0 ) { thrust::transform( vx.begin() + begin_shift, vx.end(), // input vi.begin() + begin_shift, // output - detail::divide_by_constant_and_cast(vd) // has to be done on doubles to avoid i==nx due to low precision of nvcc math; TODO: now that rand uniform has range [0,1), float might be good here? + detail::add_and_divide_by_constant_and_cast(vd, pos_shift) // has to be done on doubles to avoid i==nx due to low precision of nvcc math; TODO: now that rand uniform has range [0,1), float might be good here? ); } @@ -149,6 +150,7 @@ namespace libcloudphxx }; // ijk_ref - cell number in the refined grid + // begin_shift > 0 - used when we want to calculate ijk_ref for a subset of all particles (e.g. when initializing added particles) template void particles_t::impl::hskpng_ijk_ref(const thrust_size_t begin_shift, const bool _unravel_ijk) { @@ -156,7 +158,9 @@ namespace libcloudphxx cell_idx_hlpr_t helper(begin_shift); - if (opts_init.nx != 0) helper(x, i, opts_init.dx / opts_init.n_ref); + const real_t x_shift = mpi_rank == 0 ? 0 : opts_init.dx / opts_init.n_ref / 2.; // we shift position by half of the refined cell, because the MPI boundary cuts the refined cell in half + + if (opts_init.nx != 0) helper(x, i, opts_init.dx / opts_init.n_ref, x_shift); if (opts_init.ny != 0) helper(y, j, opts_init.dy / opts_init.n_ref); if (opts_init.nz != 0) helper(z, k, opts_init.dz / opts_init.n_ref); From 6e1fed51f07406dc30e2d677d202d997af6005fb Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 7 Jul 2022 12:38:06 +0200 Subject: [PATCH 03/38] sanity checks: dont allow untested combinations of grid refinement with non-default options --- include/libcloudph++/lgrngn/opts_init.hpp | 4 ++-- src/impl/particles_impl_init_sanity_check.ipp | 20 +++++++++++++++++++ src/particles_init.ipp | 6 +++--- src/particles_step.ipp | 2 +- 4 files changed, 26 insertions(+), 6 deletions(-) diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index cfdab7a7c..c293c86e2 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -46,7 +46,7 @@ namespace libcloudphxx // Eulerian component parameters int nx, ny, nz; real_t dx, dy, dz, dt; - unsigned int n_ref; // number of grid refinements per 'normal' cell; refined cell size is (nx, ny, nz) / n_ref; tht and rv are on refined cells, rest on normal cells; condensation is resolved on refined cells + unsigned int n_ref; // number of grid refinements per 'normal' cell; refined cell size is (dx, dy, dz) / n_ref; tht and rv are on refined cells, courants on normal cells; condensation is resolved on refined cells, rest (coalescence, advection, etc.) on normal cells // no. of substeps int sstp_cond, sstp_coal; @@ -66,7 +66,7 @@ namespace libcloudphxx // is it allowed to change dt during simulation through opts.dt bool variable_dt_switch; - // or, alternatively to sd_conc_mean, multiplicity of all SDs = const + // alternatively to sd_conc, multiplicity of all SDs = const unsigned long long sd_const_multi; // max no. of super-droplets in the system diff --git a/src/impl/particles_impl_init_sanity_check.ipp b/src/impl/particles_impl_init_sanity_check.ipp index c4b763f28..eba510ff2 100644 --- a/src/impl/particles_impl_init_sanity_check.ipp +++ b/src/impl/particles_impl_init_sanity_check.ipp @@ -152,6 +152,26 @@ namespace libcloudphxx throw std::runtime_error("opts_init.n_ref > 1 works only in 2D and 3D"); if(opts_init.n_ref == 0) throw std::runtime_error("opts_init.n_ref==0; set opts_init.n_ref=1 for no refinement"); + // guard from using grid refinement with non-default options + // this is done because these combinations were not tested, some of them might actually work + if(opts_init.n_ref > 1 && opts_init.chem_switch) + throw std::runtime_error("grid refinement is not compatible with chemistry (opts_init.n_ref > 1 && chem_switch == 1)"); + if(opts_init.n_ref > 1 && opts_init.rlx_switch) + throw std::runtime_error("grid refinement is not compatible with aerosol relaxation (opts_init.n_ref > 1 && rlx_switch == 1)"); + if(opts_init.n_ref > 1 && opts_init.src_type!=src_t::off) + throw std::runtime_error("grid refinement is not compatible with aerosol sources (opts_init.n_ref > 1 && opts_init.src_type != src_t::off)"); + if(opts_init.n_ref > 1 && opts_init.turb_adve_switch) + throw std::runtime_error("grid refinement is not compatible with SGS advection (opts_init.n_ref > 1 && opts_init.turb_adve_switch == 1)"); + if(opts_init.n_ref > 1 && opts_init.turb_cond_switch) + throw std::runtime_error("grid refinement is not compatible with SGS condensation (opts_init.n_ref > 1 && opts_init.turb_cond_switch == 1)"); + if(opts_init.n_ref > 1 && opts_init.turb_coal_switch) + throw std::runtime_error("grid refinement is not compatible with SGS coalescence kernels (opts_init.n_ref > 1 && opts_init.turb_coal_switch == 1)"); + if(opts_init.n_ref > 1 && opts_init.exact_sstp_cond) + throw std::runtime_error("grid refinement is not compatible with per-particle condensation substepping (opts_init.n_ref > 1 && opts_init.exact_sstp_cond == 1)"); + if(opts_init.n_ref > 1 && opts_init.diag_incloud_time) + throw std::runtime_error("grid refinement is not compatible with diagnosing incloud time (opts_init.n_ref > 1 && opts_init.diag_incloud_time == 1)"); + if(opts_init.n_ref > 1 && const_p==0) + throw std::runtime_error("grid refinement is not compatible with variable pressure profile (opts_init.n_ref > 1, but pressure profile was not passed to particles.init())"); } }; }; diff --git a/src/particles_init.ipp b/src/particles_init.ipp index 9f95e0c55..18ab2e830 100644 --- a/src/particles_init.ipp +++ b/src/particles_init.ipp @@ -25,14 +25,14 @@ namespace libcloudphxx const std::map > ambient_chem ) { + // is a constant pressure profile used? + pimpl->const_p = !p.is_null(); + pimpl->init_sanity_check(th, rv, rhod, p, courant_x, courant_y, courant_z, ambient_chem); // set rng seed to be used during init if(pimpl->opts_init.rng_seed_init_switch) pimpl->rng.reseed(pimpl->opts_init.rng_seed_init); - - // is a constant pressure profile used? - pimpl->const_p = !p.is_null(); // if pressure comes from a profile, sstp_tmp_p also needs to be copied between distributed memories if(pimpl->const_p && pimpl->allow_sstp_cond && pimpl->opts_init.exact_sstp_cond) pimpl->distmem_real_vctrs.insert(&pimpl->sstp_tmp_p); diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 0bf108f47..608b85748 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -176,7 +176,7 @@ namespace libcloudphxx { // calculate mean free path needed for molecular correction // NOTE: this should be don per each substep, but right now there is no logic - // that would make it easy to do in exact (per-cell) substepping + // that would make it easy to do in exact (per-particle) substepping pimpl->hskpng_mfp(); if(pimpl->opts_init.exact_sstp_cond && pimpl->sstp_cond > 1) From 3e781c733d04a9e027e232db8e441479408432c3 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 7 Jul 2022 12:46:54 +0200 Subject: [PATCH 04/38] dont allow grid refinement in multi_CUDA --- src/impl_multi_gpu/particles_multi_gpu_impl.ipp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/impl_multi_gpu/particles_multi_gpu_impl.ipp b/src/impl_multi_gpu/particles_multi_gpu_impl.ipp index 3dd334527..3df228627 100644 --- a/src/impl_multi_gpu/particles_multi_gpu_impl.ipp +++ b/src/impl_multi_gpu/particles_multi_gpu_impl.ipp @@ -44,9 +44,8 @@ namespace libcloudphxx // TODO: move these sanity checks to sanity_checks? if(glob_opts_init.chem_switch) throw std::runtime_error("multi_CUDA is not yet compatible with chemistry. Use other backend or turn off opts_init.chem_switch."); - - if(glob_opts_init.nx == 0) - throw std::runtime_error("multi_CUDA doesn't work for 0D setup."); + if(glob_opts_init.nx == 0) throw std::runtime_error("multi_CUDA doesn't work for 0D setup."); + if(glob_opts_init.n_ref > 1) throw std::runtime_error("multi_CUDA doesn't work with grid refinement."); if (!(glob_opts_init.x1 > glob_opts_init.x0 && glob_opts_init.x1 <= glob_opts_init.nx * glob_opts_init.dx)) throw std::runtime_error("!(x1 > x0 & x1 <= min(1,nx)*dx)"); From 84c3b66c540546d3d457a2e7ea8b166b3f83c860 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 7 Jul 2022 13:10:59 +0200 Subject: [PATCH 05/38] opts_init.hpp: comments about required refined arrays --- include/libcloudph++/lgrngn/opts_init.hpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index c293c86e2..121c3648b 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -46,7 +46,12 @@ namespace libcloudphxx // Eulerian component parameters int nx, ny, nz; real_t dx, dy, dz, dt; - unsigned int n_ref; // number of grid refinements per 'normal' cell; refined cell size is (dx, dy, dz) / n_ref; tht and rv are on refined cells, courants on normal cells; condensation is resolved on refined cells, rest (coalescence, advection, etc.) on normal cells + unsigned int n_ref; // number of grid refinements per 'normal' cell; refined cell size is (dx, dy, dz) / n_ref; + // condensation (and sedimentation?) is resolved on refined cells, rest (coalescence, advection, etc.) on normal cells + // input arrays on the refined grid: pressure p (if supplied), rv, tht + // input arrays that require two versions, on the normal and on the refined grid: density (rhod, rhod_ref) + // internal arrays on the refined grid: RH, T + // internal arrays that require two versions, on the normal and on the refined grid: dynamic viscosity eta (eta, eta_ref), cell volume dv (not sure...) // no. of substeps int sstp_cond, sstp_coal; From f107434c10eeed01fdf603a86385d8ba6e5fe9b9 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 7 Jul 2022 13:24:28 +0200 Subject: [PATCH 06/38] hskpng_mfp on refined grid --- src/impl/particles_impl.ipp | 8 ++++++-- src/impl/particles_impl_cond.ipp | 5 +++-- src/impl/particles_impl_hskpng_mfp.ipp | 5 +++-- src/impl/particles_impl_init_hskpng_ncell.ipp | 5 +++++ 4 files changed, 17 insertions(+), 6 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index e809b9ef9..bc1c59ff0 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -45,8 +45,8 @@ namespace libcloudphxx // member fields opts_init_t opts_init; // a copy const int n_dims; - const thrust_size_t n_cell; - const int ny_ref, nz_ref; // number of refined cells + const int nx_ref, ny_ref, nz_ref; // number of refined cells in each direction + const thrust_size_t n_cell, n_cell_ref; // total number of normal and refined cells thrust_size_t n_part, // total number of SDs n_part_old, // total number of SDs before source n_part_to_init; // number of SDs to be initialized by source @@ -221,6 +221,8 @@ namespace libcloudphxx tmp_device_real_cell, tmp_device_real_cell1, tmp_device_real_cell2, + tmp_device_real_cell_ref, + tmp_device_real_cell_ref1, &u01; // uniform random numbers between 0 and 1 // TODO: use the tmp array as rand argument? thrust_device::vector tmp_device_n_part, @@ -335,8 +337,10 @@ namespace libcloudphxx n_dims == 2 ? halo_size * (_opts_init.nz + 1): // 2D halo_size * (_opts_init.nz + 1) * _opts_init.ny // 3D ), + nx_ref(n_dims >= 1 ? (opts_init.nx - 1) * opts_init.n_ref + 1 : 0), // TODO: with MPI this is probably not true ny_ref(n_dims == 3 ? (opts_init.ny - 1) * opts_init.n_ref + 1 : 0), nz_ref(n_dims >= 2 ? (opts_init.nz - 1) * opts_init.n_ref + 1 : 0), + n_cell_ref(m1(nx_ref) * m1(ny_ref) * m1(nz_ref)), w_LS(_opts_init.w_LS), SGS_mix_len(_opts_init.SGS_mix_len), adve_scheme(_opts_init.adve_scheme), diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index c3fdfa618..b9e622c39 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -18,8 +18,9 @@ namespace libcloudphxx namespace arg = thrust::placeholders; - thrust_device::vector &lambda_D(tmp_device_real_cell1); // real_cell used in cond.ipp - thrust_device::vector &lambda_K(tmp_device_real_cell2); // real_cell used in cond.ipp + // needs to be the same as in hskpng_mfp + thrust_device::vector &lambda_D(opts_init.n_ref > 1 ? tmp_device_real_cell_ref : tmp_device_real_cell1); + thrust_device::vector &lambda_K(opts_init.n_ref > 1 ? tmp_device_real_cell_ref1 : tmp_device_real_cell2); // --- calc liquid water content before cond --- hskpng_sort(); diff --git a/src/impl/particles_impl_hskpng_mfp.ipp b/src/impl/particles_impl_hskpng_mfp.ipp index 1cf02535b..9ac12abaa 100644 --- a/src/impl/particles_impl_hskpng_mfp.ipp +++ b/src/impl/particles_impl_hskpng_mfp.ipp @@ -41,8 +41,9 @@ namespace libcloudphxx template void particles_t::impl::hskpng_mfp() { - thrust_device::vector &lambda_D(tmp_device_real_cell1); // real_cell used in cond.ipp - thrust_device::vector &lambda_K(tmp_device_real_cell2); // real_cell used in cond.ipp + // tmp_device_real_cell cant be used here, it is used in cond.ipp + thrust_device::vector &lambda_D(opts_init.n_ref > 1 ? tmp_device_real_cell_ref : tmp_device_real_cell1); + thrust_device::vector &lambda_K(opts_init.n_ref > 1 ? tmp_device_real_cell_ref1 : tmp_device_real_cell2); thrust::transform(T.begin(), T.end(), lambda_D.begin(), detail::common__mean_free_path__lambda_D()); thrust::transform(T.begin(), T.end(), p.begin(), lambda_K.begin(), detail::common__mean_free_path__lambda_K()); diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index 2b92d6601..7191b269a 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -35,6 +35,11 @@ namespace libcloudphxx sstp_tmp_th.resize(n_cell); sstp_tmp_rh.resize(n_cell); } + if(opts_init.n_ref > 1) + { + tmp_device_real_cell_ref.resize(n_cell_ref); + tmp_device_real_cell_ref1.resize(n_cell_ref); + } } }; }; From 881adf5386f129d49cd9a9828a58ce53a4810ecc Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 7 Jul 2022 13:47:33 +0200 Subject: [PATCH 07/38] ijk_ref_hlpr --- src/impl/particles_impl.ipp | 3 +++ src/impl/particles_impl_hskpng_ijk.ipp | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index bc1c59ff0..20de99d58 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -110,6 +110,8 @@ namespace libcloudphxx i, j, k, ijk, // Eulerian grid cell indices (always zero for 0D) ijk_ref, // ijk in refined cells sorted_id, sorted_ijk; + + thrust_device::vector &ijk_ref_hlpr; // ijk_ref if refinement is done, ijk otherwise // Arakawa-C grid helper vars thrust_device::vector @@ -341,6 +343,7 @@ namespace libcloudphxx ny_ref(n_dims == 3 ? (opts_init.ny - 1) * opts_init.n_ref + 1 : 0), nz_ref(n_dims >= 2 ? (opts_init.nz - 1) * opts_init.n_ref + 1 : 0), n_cell_ref(m1(nx_ref) * m1(ny_ref) * m1(nz_ref)), + ijk_ref_hlpr(opts_init.n_ref > 1 ? ijk_ref : ijk), w_LS(_opts_init.w_LS), SGS_mix_len(_opts_init.SGS_mix_len), adve_scheme(_opts_init.adve_scheme), diff --git a/src/impl/particles_impl_hskpng_ijk.ipp b/src/impl/particles_impl_hskpng_ijk.ipp index d9184239b..3b18e9cc1 100644 --- a/src/impl/particles_impl_hskpng_ijk.ipp +++ b/src/impl/particles_impl_hskpng_ijk.ipp @@ -154,7 +154,7 @@ namespace libcloudphxx template void particles_t::impl::hskpng_ijk_ref(const thrust_size_t begin_shift, const bool _unravel_ijk) { - if(opts_init.n_ref <= 1) return; + if(opts_init.n_ref == 1) return; cell_idx_hlpr_t helper(begin_shift); From 6dacb40e00ecf8cec0ca3f3f447840525e789814 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 7 Jul 2022 13:58:09 +0200 Subject: [PATCH 08/38] init size of refined internal arrays --- src/impl/particles_impl.ipp | 12 +++++++++++- src/impl/particles_impl_init_hskpng_ncell.ipp | 13 +++++++++++-- 2 files changed, 22 insertions(+), 3 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 20de99d58..3a3ae07e9 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -126,6 +126,15 @@ namespace libcloudphxx count_mom; // statistical moment // TODO (perhaps tmp_device_real_cell could be referenced?) thrust_size_t count_n; + // moment-counting stuff on the refined grid + thrust_device::vector + count_ijk_ref; + thrust_device::vector + count_num_ref; + thrust_device::vector + count_mom_ref; + thrust_size_t count_n_ref; + // Eulerian-Lagrangian interface vars thrust_device::vector rhod, // dry air density @@ -151,6 +160,7 @@ namespace libcloudphxx p, // pressure [Pa] RH, // relative humisity eta,// dynamic viscosity + eta_ref,// dynamic viscosity on the refined grid diss_rate; // turbulent kinetic energy dissipation rate thrust_device::vector w_LS; // large-scale subsidence velocity profile @@ -342,7 +352,7 @@ namespace libcloudphxx nx_ref(n_dims >= 1 ? (opts_init.nx - 1) * opts_init.n_ref + 1 : 0), // TODO: with MPI this is probably not true ny_ref(n_dims == 3 ? (opts_init.ny - 1) * opts_init.n_ref + 1 : 0), nz_ref(n_dims >= 2 ? (opts_init.nz - 1) * opts_init.n_ref + 1 : 0), - n_cell_ref(m1(nx_ref) * m1(ny_ref) * m1(nz_ref)), + n_cell_ref(m1(nx_ref) * m1(ny_ref) * m1(nz_ref)), // NOTE: needs to be equal to n_cell for n_ref == 1 ijk_ref_hlpr(opts_init.n_ref > 1 ? ijk_ref : ijk), w_LS(_opts_init.w_LS), SGS_mix_len(_opts_init.SGS_mix_len), diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index 7191b269a..b3f94264c 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -13,8 +13,6 @@ namespace libcloudphxx void particles_t::impl::init_hskpng_ncell() { // memory allocation, p already resized in init_sync() - T.resize(n_cell); - RH.resize(n_cell); eta.resize(n_cell); count_ijk.resize(n_cell); @@ -35,10 +33,21 @@ namespace libcloudphxx sstp_tmp_th.resize(n_cell); sstp_tmp_rh.resize(n_cell); } + + // arrays on the refined grid if(opts_init.n_ref > 1) { + T.resize(n_cell_ref); + RH.resize(n_cell_ref); + eta_ref.resize(n_cell); + tmp_device_real_cell_ref.resize(n_cell_ref); tmp_device_real_cell_ref1.resize(n_cell_ref); + + count_ijk_ref.resize(n_cell_ref); + count_num_ref.resize(n_cell_ref); + count_mom_ref.resize(n_cell_ref); + count_n_ref = 0; } } }; From b5db9df72151d8b7c2ab1e76fe27f29a44a847a7 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 7 Jul 2022 14:09:14 +0200 Subject: [PATCH 09/38] refined grid: diag_RH --- src/impl/particles_impl.ipp | 34 ++++++++++++++++------------------ src/particles_diag.ipp | 10 +++++----- 2 files changed, 21 insertions(+), 23 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 3a3ae07e9..4f5eb0a88 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -110,30 +110,19 @@ namespace libcloudphxx i, j, k, ijk, // Eulerian grid cell indices (always zero for 0D) ijk_ref, // ijk in refined cells sorted_id, sorted_ijk; - - thrust_device::vector &ijk_ref_hlpr; // ijk_ref if refinement is done, ijk otherwise // Arakawa-C grid helper vars thrust_device::vector lft, rgt, abv, blw, fre, hnd; // TODO: could be reused after advection! - // moment-counting stuff + // moment-counting stuff on normal and refined grid thrust_device::vector - count_ijk; // key-value pair for sorting particles by cell index + count_ijk, count_ijk_ref; // key-value pair for sorting particles by cell index thrust_device::vector - count_num; // number of particles in a given grid cell + count_num, count_num_ref; // number of particles in a given grid cell thrust_device::vector - count_mom; // statistical moment // TODO (perhaps tmp_device_real_cell could be referenced?) - thrust_size_t count_n; - - // moment-counting stuff on the refined grid - thrust_device::vector - count_ijk_ref; - thrust_device::vector - count_num_ref; - thrust_device::vector - count_mom_ref; - thrust_size_t count_n_ref; + count_mom, count_mom_ref; // statistical moment // TODO (perhaps tmp_device_real_cell could be referenced?) + thrust_size_t count_n, count_n_ref; // Eulerian-Lagrangian interface vars thrust_device::vector @@ -195,6 +184,11 @@ namespace libcloudphxx thrust::host_vector > l2e; + // helpers that reference refined arrays if refinement is done and non-refined otherwise + thrust_device::vector &ijk_ref_hlpr, &count_ijk_ref_hlpr; + thrust_device::vector &count_num_ref_hlpr; + thrust_device::vector &count_mom_ref_hlpr, &eta_ref_hlpr; + // chem stuff // TODO: consider changing the unit to AMU or alike (very small numbers!) std::vector::iterator > @@ -353,13 +347,17 @@ namespace libcloudphxx ny_ref(n_dims == 3 ? (opts_init.ny - 1) * opts_init.n_ref + 1 : 0), nz_ref(n_dims >= 2 ? (opts_init.nz - 1) * opts_init.n_ref + 1 : 0), n_cell_ref(m1(nx_ref) * m1(ny_ref) * m1(nz_ref)), // NOTE: needs to be equal to n_cell for n_ref == 1 - ijk_ref_hlpr(opts_init.n_ref > 1 ? ijk_ref : ijk), w_LS(_opts_init.w_LS), SGS_mix_len(_opts_init.SGS_mix_len), adve_scheme(_opts_init.adve_scheme), allow_sstp_cond(_opts_init.sstp_cond > 1 || _opts_init.variable_dt_switch), allow_sstp_chem(_opts_init.sstp_chem > 1 || _opts_init.variable_dt_switch), - pure_const_multi (((_opts_init.sd_conc) == 0) && (_opts_init.sd_const_multi > 0 || _opts_init.dry_sizes.size() > 0)) // coal prob can be greater than one only in sd_conc simulations + pure_const_multi (((_opts_init.sd_conc) == 0) && (_opts_init.sd_const_multi > 0 || _opts_init.dry_sizes.size() > 0)), // coal prob can be greater than one only in sd_conc simulations + ijk_ref_hlpr (opts_init.n_ref > 1 ? ijk_ref : ijk), + count_ijk_ref_hlpr(opts_init.n_ref > 1 ? count_ijk_ref : count_ijk), + count_num_ref_hlpr(opts_init.n_ref > 1 ? count_num_ref : count_num), + count_mom_ref_hlpr(opts_init.n_ref > 1 ? count_mom_ref : count_mom), + eta_ref_hlpr (opts_init.n_ref > 1 ? eta_ref : eta) { // set 0 dev_count to mark that its not a multi_CUDA spawn diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index a948143a3..bf6588f60 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -147,14 +147,14 @@ namespace libcloudphxx pimpl->hskpng_Tpr(); thrust::copy( - pimpl->RH.begin(), - pimpl->RH.end(), - pimpl->count_mom.begin() + pimpl->RH_ref_hlpr.begin(), + pimpl->RH_ref_hlpr.end(), + pimpl->count_mom_ref_hlpr.begin() ); // RH defined in all cells - pimpl->count_n = pimpl->n_cell; - thrust::sequence(pimpl->count_ijk.begin(), pimpl->count_ijk.end()); + pimpl->count_n_ref = pimpl->n_cell_ref; + thrust::sequence(pimpl->count_ijk_ref_hlpr.begin(), pimpl->count_ijk_ref_hlpr.end()); } // records super-droplet concentration per grid cell From e1eeedf1ee807c2a1c86eb2d81df46cb49fbf8fe Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 7 Jul 2022 14:26:22 +0200 Subject: [PATCH 10/38] refinement: T and RH on refined grid --- src/impl/particles_impl.ipp | 11 +++++--- src/impl/particles_impl_chem_dissoc.ipp | 4 +-- src/impl/particles_impl_chem_henry.ipp | 2 +- src/impl/particles_impl_chem_react.ipp | 2 +- src/impl/particles_impl_chem_strength.ipp | 4 +-- src/impl/particles_impl_cond.ipp | 6 ++--- src/impl/particles_impl_hskpng_Tpr.ipp | 16 ++++++------ src/impl/particles_impl_hskpng_mfp.ipp | 4 +-- src/impl/particles_impl_hskpng_vterm.ipp | 4 +-- src/impl/particles_impl_init_hskpng_ncell.ipp | 4 +-- src/impl/particles_impl_init_wet.ipp | 4 +-- .../particles_impl_update_incloud_time.ipp | 4 +-- src/impl/particles_impl_update_th_rv.ipp | 2 +- src/particles_diag.ipp | 26 +++++++++---------- 14 files changed, 48 insertions(+), 45 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 4f5eb0a88..acfd5a5ec 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -144,12 +144,13 @@ namespace libcloudphxx // map of the accumulated volume/volume/mass of water/dry/chem that fell out of the domain std::map output_puddle; + // _ref means that the array is on the refined grid (which is the same as the normal grid if n_ref==1), some only need to be on one of the grids, some need to be on both thrust_device::vector - T, // temperature [K] - p, // pressure [Pa] - RH, // relative humisity + T_ref, // temperature [K] + p_ref, // pressure [Pa] + RH_ref, // relative humisity eta,// dynamic viscosity - eta_ref,// dynamic viscosity on the refined grid + eta_ref,// dynamic viscosity diss_rate; // turbulent kinetic energy dissipation rate thrust_device::vector w_LS; // large-scale subsidence velocity profile @@ -185,6 +186,8 @@ namespace libcloudphxx > l2e; // helpers that reference refined arrays if refinement is done and non-refined otherwise + // used to conserve memory if n_ref==1 (no refinement) + // needed only for arrays that have both a refined and a normal version, not needed for arrays only on the refined grid thrust_device::vector &ijk_ref_hlpr, &count_ijk_ref_hlpr; thrust_device::vector &count_num_ref_hlpr; thrust_device::vector &count_mom_ref_hlpr, &eta_ref_hlpr; diff --git a/src/impl/particles_impl_chem_dissoc.ipp b/src/impl/particles_impl_chem_dissoc.ipp index cfa0a691b..1da66ec9b 100644 --- a/src/impl/particles_impl_chem_dissoc.ipp +++ b/src/impl/particles_impl_chem_dissoc.ipp @@ -172,11 +172,11 @@ namespace libcloudphxx zip_it_t(thrust::make_tuple( chem_bgn[SO2], chem_bgn[CO2], chem_bgn[HNO3], chem_bgn[NH3], chem_bgn[S_VI], V.begin(), - thrust::make_permutation_iterator(T.begin(), ijk.begin()))), // input - begin + thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()))), // input - begin zip_it_t(thrust::make_tuple( chem_end[SO2], chem_end[CO2], chem_end[HNO3], chem_end[NH3], chem_end[S_VI], V.end(), - thrust::make_permutation_iterator(T.end(), ijk.end()))), // input - end + thrust::make_permutation_iterator(T_ref.end(), ijk_ref_hlpr.end()))), // input - end chem_flag.begin(), // stencil chem_bgn[H], // output detail::chem_electroneutral(), // op diff --git a/src/impl/particles_impl_chem_henry.ipp b/src/impl/particles_impl_chem_henry.ipp index fabaa8f44..a496d1b72 100644 --- a/src/impl/particles_impl_chem_henry.ipp +++ b/src/impl/particles_impl_chem_henry.ipp @@ -335,7 +335,7 @@ namespace libcloudphxx V.begin(), V.end(), // input - 1st arg thrust::make_zip_iterator(thrust::make_tuple( // input - 2nd arg thrust::make_permutation_iterator(p.begin(), ijk.begin()), - thrust::make_permutation_iterator(T.begin(), ijk.begin()), + thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()), thrust::make_permutation_iterator(ambient_chem[(chem_species_t)i].begin(), ijk.begin()), chem_bgn[i], rw2.begin(), diff --git a/src/impl/particles_impl_chem_react.ipp b/src/impl/particles_impl_chem_react.ipp index b95e750fe..c48341d88 100644 --- a/src/impl/particles_impl_chem_react.ipp +++ b/src/impl/particles_impl_chem_react.ipp @@ -278,7 +278,7 @@ namespace libcloudphxx detail::chem_rhs( dt, V, - thrust::make_permutation_iterator(T.begin(), ijk.begin()), + thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()), chem_bgn[H], chem_flag ), // TODO: make it an impl member field diff --git a/src/impl/particles_impl_chem_strength.ipp b/src/impl/particles_impl_chem_strength.ipp index 43338ea9f..3d8972611 100644 --- a/src/impl/particles_impl_chem_strength.ipp +++ b/src/impl/particles_impl_chem_strength.ipp @@ -97,10 +97,10 @@ namespace libcloudphxx thrust::transform( zip_it_t(thrust::make_tuple( chem_bgn[SO2], chem_bgn[CO2], chem_bgn[HNO3], chem_bgn[NH3], chem_bgn[S_VI], chem_bgn[H], - V.begin(), thrust::make_permutation_iterator(T.begin(), ijk.begin()), rw2.begin())), // input - begin + V.begin(), thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()), rw2.begin())), // input - begin zip_it_t(thrust::make_tuple( chem_end[SO2], chem_end[CO2], chem_end[HNO3], chem_end[NH3], chem_end[S_VI], chem_end[H], - V.end(), thrust::make_permutation_iterator(T.end(), ijk.end()), rw2.end())), // input - end + V.end(), thrust::make_permutation_iterator(T_ref.end(), ijk_ref_hlpr.end()), rw2.end())), // input - end chem_flag.begin(), // output detail::set_chem_flag() // op ); diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index b9e622c39..495c6926e 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -43,7 +43,7 @@ namespace libcloudphxx 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(T_ref.begin(), ijk_ref_hlpr.begin()), thrust::make_permutation_iterator(eta.begin(), ijk.begin()), rd3.begin(), kpa.begin(), @@ -59,7 +59,7 @@ namespace libcloudphxx thrust_device::vector &RH_plus_ssp(tmp_device_real_part2); thrust::transform( ssp.begin(), ssp.end(), - thrust::make_permutation_iterator(RH.begin(), ijk.begin()), + thrust::make_permutation_iterator(RH_ref.begin(), ijk_ref_hlpr.begin()), RH_plus_ssp.begin(), arg::_1 + arg::_2 ); @@ -84,7 +84,7 @@ namespace libcloudphxx thrust::make_tuple( hlpr_zip_iter, thrust::make_permutation_iterator(p.begin(), ijk.begin()), - thrust::make_permutation_iterator(RH.begin(), ijk.begin()) + thrust::make_permutation_iterator(RH_ref.begin(), ijk_ref_hlpr.begin()) ) ), rw2.begin(), // output diff --git a/src/impl/particles_impl_hskpng_Tpr.ipp b/src/impl/particles_impl_hskpng_Tpr.ipp index 38433bcad..8957ff86b 100644 --- a/src/impl/particles_impl_hskpng_Tpr.ipp +++ b/src/impl/particles_impl_hskpng_Tpr.ipp @@ -179,7 +179,7 @@ namespace libcloudphxx thrust::transform( th.begin(), th.end(), // input - first arg rhod.begin(), // input - second arg - T.begin(), // output + T_ref.begin(), // output detail::common__theta_dry__T_rhod() ); } @@ -189,7 +189,7 @@ namespace libcloudphxx thrust::transform( th.begin(), th.end(), // input - first arg thrust::make_zip_iterator(thrust::make_tuple(rv.begin(), p.begin())), // input - second and third args - T.begin(), // output + T_ref.begin(), // output detail::common__theta_dry__T_p() ); } @@ -207,17 +207,17 @@ namespace libcloudphxx { // p = common::theta_dry::p(rhod, r, T); thrust::transform( - zip_it_t(thrust::make_tuple(rhod.begin(), rv.begin(), T.begin())), // input - begin - zip_it_t(thrust::make_tuple(rhod.end(), rv.end(), T.end() )), // input - end + zip_it_t(thrust::make_tuple(rhod.begin(), rv.begin(), T_ref.begin())), // input - begin + zip_it_t(thrust::make_tuple(rhod.end(), rv.end(), T_ref.end() )), // input - end p.begin(), // output detail::common__theta_dry__p() ); } thrust::transform( - zip_it_t(thrust::make_tuple(p.begin(), rv.begin(), T.begin())), // input - begin - zip_it_t(thrust::make_tuple(p.end(), rv.end(), T.end() )), // input - end - RH.begin(), // output + zip_it_t(thrust::make_tuple(p.begin(), rv.begin(), T_ref.begin())), // input - begin + zip_it_t(thrust::make_tuple(p.end(), rv.end(), T_ref.end() )), // input - end + RH_ref.begin(), // output detail::RH(opts_init.RH_formula) ); } @@ -225,7 +225,7 @@ namespace libcloudphxx // dynamic viscosity { thrust::transform( - T.begin(), T.end(), // 1st arg + T_ref.begin(), T_ref.end(), // 1st arg eta.begin(), // output detail::common__vterm__visc() ); diff --git a/src/impl/particles_impl_hskpng_mfp.ipp b/src/impl/particles_impl_hskpng_mfp.ipp index 9ac12abaa..b3586682b 100644 --- a/src/impl/particles_impl_hskpng_mfp.ipp +++ b/src/impl/particles_impl_hskpng_mfp.ipp @@ -45,8 +45,8 @@ namespace libcloudphxx thrust_device::vector &lambda_D(opts_init.n_ref > 1 ? tmp_device_real_cell_ref : tmp_device_real_cell1); thrust_device::vector &lambda_K(opts_init.n_ref > 1 ? tmp_device_real_cell_ref1 : tmp_device_real_cell2); - thrust::transform(T.begin(), T.end(), lambda_D.begin(), detail::common__mean_free_path__lambda_D()); - thrust::transform(T.begin(), T.end(), p.begin(), lambda_K.begin(), detail::common__mean_free_path__lambda_K()); + thrust::transform(T_ref.begin(), T_ref.end(), lambda_D.begin(), detail::common__mean_free_path__lambda_D()); + thrust::transform(T_ref.begin(), T_ref.end(), p.begin(), lambda_K.begin(), detail::common__mean_free_path__lambda_K()); } }; }; diff --git a/src/impl/particles_impl_hskpng_vterm.ipp b/src/impl/particles_impl_hskpng_vterm.ipp index 765d7768c..aaae7e3e2 100644 --- a/src/impl/particles_impl_hskpng_vterm.ipp +++ b/src/impl/particles_impl_hskpng_vterm.ipp @@ -172,7 +172,7 @@ namespace libcloudphxx thrust::transform_if( rw2.begin(), rw2.end(), // input - 1st arg zip_it_t(thrust::make_tuple( - thrust::make_permutation_iterator(T.begin(), ijk.begin()), + thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()), thrust::make_permutation_iterator(p.begin(), ijk.begin()), thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin(), ijk.begin()) @@ -219,7 +219,7 @@ namespace libcloudphxx thrust::transform( rw2.begin(), rw2.end(), // input - 1st arg zip_it_t(thrust::make_tuple( - thrust::make_permutation_iterator(T.begin(), ijk.begin()), + thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()), thrust::make_permutation_iterator(p.begin(), ijk.begin()), thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin(), ijk.begin()) diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index b3f94264c..4a40cb775 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -37,8 +37,8 @@ namespace libcloudphxx // arrays on the refined grid if(opts_init.n_ref > 1) { - T.resize(n_cell_ref); - RH.resize(n_cell_ref); + T_ref.resize(n_cell_ref); + RH_ref.resize(n_cell_ref); eta_ref.resize(n_cell); tmp_device_real_cell_ref.resize(n_cell_ref); diff --git a/src/impl/particles_impl_init_wet.ipp b/src/impl/particles_impl_init_wet.ipp index 77646cd47..96be622c2 100644 --- a/src/impl/particles_impl_init_wet.ipp +++ b/src/impl/particles_impl_init_wet.ipp @@ -63,8 +63,8 @@ namespace libcloudphxx zip_it_t zip_it(thrust::make_tuple( rd3.begin() + n_part_old, kpa.begin() + n_part_old, - pi_t(RH.begin(), ijk.begin() + n_part_old), - pi_t(T.begin(), ijk.begin() + n_part_old) + pi_t(RH_ref.begin(), ijk_ref_hlpr.begin() + n_part_old), + pi_t(T_ref.begin(), ijk_ref_hlpr.begin() + n_part_old) )); thrust::transform( diff --git a/src/impl/particles_impl_update_incloud_time.ipp b/src/impl/particles_impl_update_incloud_time.ipp index 75dbec831..369f56c39 100644 --- a/src/impl/particles_impl_update_incloud_time.ipp +++ b/src/impl/particles_impl_update_incloud_time.ipp @@ -44,8 +44,8 @@ namespace libcloudphxx thrust::make_zip_iterator(make_tuple( kpa.begin(), thrust::make_permutation_iterator( - T.begin(), - ijk.begin() + T_ref.begin(), + ijk_ref_hlpr.begin() ) )), // input - 2nd arg rc2.begin(), // output diff --git a/src/impl/particles_impl_update_th_rv.ipp b/src/impl/particles_impl_update_th_rv.ipp index f8b70bd94..af575e038 100644 --- a/src/impl/particles_impl_update_th_rv.ipp +++ b/src/impl/particles_impl_update_th_rv.ipp @@ -77,7 +77,7 @@ namespace libcloudphxx thrust::make_transform_iterator( zip_it_t(thrust::make_tuple( drv.begin(), // - T.begin(), // dth = drv * d_th_d_rv(T, th) + T_ref.begin(), // dth = drv * d_th_d_rv(T, th) th.begin() // )), detail::dth() diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index bf6588f60..b3f1aff61 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -130,14 +130,14 @@ namespace libcloudphxx pimpl->hskpng_Tpr(); thrust::copy( - pimpl->T.begin(), - pimpl->T.end(), - pimpl->count_mom.begin() + pimpl->T_ref.begin(), + pimpl->T_ref.end(), + pimpl->count_mom_ref_hlpr.begin() ); // T defined in all cells - pimpl->count_n = pimpl->n_cell; - thrust::sequence(pimpl->count_ijk.begin(), pimpl->count_ijk.end()); + pimpl->count_n_ref = pimpl->n_cell_ref; + thrust::sequence(pimpl->count_ijk_ref_hlpr.begin(), pimpl->count_ijk_ref_hlpr.end()); } // records relative humidity @@ -147,8 +147,8 @@ namespace libcloudphxx pimpl->hskpng_Tpr(); thrust::copy( - pimpl->RH_ref_hlpr.begin(), - pimpl->RH_ref_hlpr.end(), + pimpl->RH_ref.begin(), + pimpl->RH_ref.end(), pimpl->count_mom_ref_hlpr.begin() ); @@ -256,12 +256,12 @@ namespace libcloudphxx thrust::make_zip_iterator(make_tuple( pimpl->kpa.begin(), thrust::make_permutation_iterator( - pimpl->T.begin(), - pimpl->ijk.begin() + pimpl->T_ref.begin(), + pimpl->ijk_ref_hlpr.begin() ), thrust::make_permutation_iterator( - pimpl->RH.begin(), - pimpl->ijk.begin() + pimpl->RH_ref.begin(), + pimpl->ijk_ref_hlpr.begin() ) )), // input - 2nd arg RH_minus_Sc.begin(), // output @@ -285,8 +285,8 @@ namespace libcloudphxx thrust::make_zip_iterator(make_tuple( pimpl->kpa.begin(), thrust::make_permutation_iterator( - pimpl->T.begin(), - pimpl->ijk.begin() + pimpl->T_ref.begin(), + pimpl->ijk_ref_hlpr.begin() ) )), // input - 2nd arg rc2.begin(), // output From 784318231e87054c3b5b9f67d373397de8593155 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 7 Jul 2022 14:44:49 +0200 Subject: [PATCH 11/38] refinement: eta p rhod on refined grids --- src/impl/particles_impl.ipp | 9 ++++---- src/impl/particles_impl_hskpng_Tpr.ipp | 23 +++++++++++-------- src/impl/particles_impl_init_hskpng_ncell.ipp | 2 +- src/impl/particles_impl_init_sync.ipp | 7 +++--- 4 files changed, 23 insertions(+), 18 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index acfd5a5ec..5b88cf99f 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -126,9 +126,10 @@ namespace libcloudphxx // Eulerian-Lagrangian interface vars thrust_device::vector - rhod, // dry air density - th, // potential temperature (dry) - rv, // water vapour mixing ratio + rhod, // dry air density + rhod_ref, // same on the refined grid + th_ref, // potential temperature (dry) + rv_ref, // water vapour mixing ratio sstp_tmp_chem_0, // ditto for trace gases sstp_tmp_chem_1, // ditto for trace gases sstp_tmp_chem_2, // ditto for trace gases @@ -149,7 +150,7 @@ namespace libcloudphxx T_ref, // temperature [K] p_ref, // pressure [Pa] RH_ref, // relative humisity - eta,// dynamic viscosity + //eta,// dynamic viscosity, commented to show that we dont know how to calculate it! require input of non-refined tht? average from refined cells? eta_ref,// dynamic viscosity diss_rate; // turbulent kinetic energy dissipation rate diff --git a/src/impl/particles_impl_hskpng_Tpr.ipp b/src/impl/particles_impl_hskpng_Tpr.ipp index 8957ff86b..6b8a49af8 100644 --- a/src/impl/particles_impl_hskpng_Tpr.ipp +++ b/src/impl/particles_impl_hskpng_Tpr.ipp @@ -177,8 +177,8 @@ namespace libcloudphxx { // T = common::theta_dry::T(th, rhod); thrust::transform( - th.begin(), th.end(), // input - first arg - rhod.begin(), // input - second arg + th_ref.begin(), th_ref.end(), // input - first arg + rhod_ref.begin(), // input - second arg T_ref.begin(), // output detail::common__theta_dry__T_rhod() ); @@ -187,8 +187,8 @@ namespace libcloudphxx { // T = th * exner(p_tot) thrust::transform( - th.begin(), th.end(), // input - first arg - thrust::make_zip_iterator(thrust::make_tuple(rv.begin(), p.begin())), // input - second and third args + th_ref.begin(), th_ref.end(), // input - first arg + thrust::make_zip_iterator(thrust::make_tuple(rv_ref.begin(), p_ref.begin())), // input - second and third args T_ref.begin(), // output detail::common__theta_dry__T_p() ); @@ -207,16 +207,16 @@ namespace libcloudphxx { // p = common::theta_dry::p(rhod, r, T); thrust::transform( - zip_it_t(thrust::make_tuple(rhod.begin(), rv.begin(), T_ref.begin())), // input - begin - zip_it_t(thrust::make_tuple(rhod.end(), rv.end(), T_ref.end() )), // input - end - p.begin(), // output + zip_it_t(thrust::make_tuple(rhod_ref.begin(), rv_ref.begin(), T_ref.begin())), // input - begin + zip_it_t(thrust::make_tuple(rhod_ref.end(), rv_ref.end(), T_ref.end() )), // input - end + p_ref.begin(), // output detail::common__theta_dry__p() ); } thrust::transform( - zip_it_t(thrust::make_tuple(p.begin(), rv.begin(), T_ref.begin())), // input - begin - zip_it_t(thrust::make_tuple(p.end(), rv.end(), T_ref.end() )), // input - end + zip_it_t(thrust::make_tuple(p_ref.begin(), rv_ref.begin(), T_ref.begin())), // input - begin + zip_it_t(thrust::make_tuple(p_ref.end(), rv_ref.end(), T_ref.end() )), // input - end RH_ref.begin(), // output detail::RH(opts_init.RH_formula) ); @@ -224,11 +224,14 @@ namespace libcloudphxx // dynamic viscosity { + // on refined grid thrust::transform( T_ref.begin(), T_ref.end(), // 1st arg - eta.begin(), // output + eta_ref.begin(), // output detail::common__vterm__visc() ); + // on normal grid: how if we dont know T, because we dont know th? + // TODO: use averaging from refined to normal? } // adjusting dv if using a parcel set-up (1kg of dry air) diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index 4a40cb775..25e84c732 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -39,7 +39,7 @@ namespace libcloudphxx { T_ref.resize(n_cell_ref); RH_ref.resize(n_cell_ref); - eta_ref.resize(n_cell); + eta_ref.resize(n_cell_ref); tmp_device_real_cell_ref.resize(n_cell_ref); tmp_device_real_cell_ref1.resize(n_cell_ref); diff --git a/src/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index 616792eab..3c7d4c405 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/src/impl/particles_impl_init_sync.ipp @@ -14,9 +14,10 @@ namespace libcloudphxx { // memory allocation for scalar fields rhod.resize(n_cell); - p.resize(n_cell); - th.resize(n_cell); - rv.resize(n_cell); + rhod_ref.resize(n_cell_ref); + p_ref.resize(n_cell_ref); + th_ref.resize(n_cell_ref); + rv_ref.resize(n_cell_ref); if(opts_init.chem_switch) for (int i = 0; i < chem_gas_n; ++i) ambient_chem[(chem_species_t)i].resize(n_cell); From 219ea049d13371886fcb191103559e0ecc3cf34d Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 7 Jul 2022 16:53:46 +0200 Subject: [PATCH 12/38] refinement wip --- src/impl/particles_impl.ipp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 5b88cf99f..adca4095f 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -147,11 +147,11 @@ namespace libcloudphxx // _ref means that the array is on the refined grid (which is the same as the normal grid if n_ref==1), some only need to be on one of the grids, some need to be on both thrust_device::vector - T_ref, // temperature [K] - p_ref, // pressure [Pa] - RH_ref, // relative humisity + T_ref, // temperature [K] + p_ref, // pressure [Pa] + RH_ref, // relative humisity //eta,// dynamic viscosity, commented to show that we dont know how to calculate it! require input of non-refined tht? average from refined cells? - eta_ref,// dynamic viscosity + eta_ref, // dynamic viscosity diss_rate; // turbulent kinetic energy dissipation rate thrust_device::vector w_LS; // large-scale subsidence velocity profile @@ -191,7 +191,7 @@ namespace libcloudphxx // needed only for arrays that have both a refined and a normal version, not needed for arrays only on the refined grid thrust_device::vector &ijk_ref_hlpr, &count_ijk_ref_hlpr; thrust_device::vector &count_num_ref_hlpr; - thrust_device::vector &count_mom_ref_hlpr, &eta_ref_hlpr; + thrust_device::vector &count_mom_ref_hlpr, &eta_ref_hlpr, &rhod_ref_hlpr; // chem stuff // TODO: consider changing the unit to AMU or alike (very small numbers!) @@ -361,7 +361,8 @@ namespace libcloudphxx count_ijk_ref_hlpr(opts_init.n_ref > 1 ? count_ijk_ref : count_ijk), count_num_ref_hlpr(opts_init.n_ref > 1 ? count_num_ref : count_num), count_mom_ref_hlpr(opts_init.n_ref > 1 ? count_mom_ref : count_mom), - eta_ref_hlpr (opts_init.n_ref > 1 ? eta_ref : eta) + eta_ref_hlpr (opts_init.n_ref > 1 ? eta_ref : eta), + rhod_ref_hlpr (opts_init.n_ref > 1 ? rhod_ref : rhod) { // set 0 dev_count to mark that its not a multi_CUDA spawn From e115c3001f1e0c560449d8ef7bc1d481c08da9d2 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 8 Jul 2022 15:15:44 +0200 Subject: [PATCH 13/38] wip on grid refinement classes --- src/detail/grid_refinement_helpers.hpp | 87 ++++++++++++++++++++++++++ src/impl/particles_impl.ipp | 30 +++++---- src/particles.tpp | 1 + 3 files changed, 102 insertions(+), 16 deletions(-) create mode 100644 src/detail/grid_refinement_helpers.hpp diff --git a/src/detail/grid_refinement_helpers.hpp b/src/detail/grid_refinement_helpers.hpp new file mode 100644 index 000000000..68367f15f --- /dev/null +++ b/src/detail/grid_refinement_helpers.hpp @@ -0,0 +1,87 @@ +// helper classes to store arrays associated with grid refinement + +#pragma once + +namespace libcloudphxx +{ + namespace lgrngn + { + template + class ref_common + { + private: + const bool ref_flag; // true if refinement is actually done + thrust_device::vector arr, arr_ref; // actual data, arr_ref initialized only if refinement is actually done + + protected: + // ctor + ref_common(const bool ref_flag): + ref_flag(ref_flag) + {} + + auto begin() + { + return arr.begin(); + } + auto begin_ref() + { + return ref_flag ? arr_ref.begin() : arr.begin(); + } + auto end() + { + return arr.end(); + } + auto end_ref() + { + return ref_flag ? arr_ref.end() : arr.end(); + } + void resize(thrust_size_t n, thrust_size_t n_ref) + { + arr.resize(n); + if(ref_flag) arr_ref.resize(n_ref); + } + void reserve(thrust_size_t n, thrust_size_t n_ref) + { + arr.reserve(n); + if(ref_flag) arr_ref.reserve(n_ref); + } + } + + // for arrays of the size of the grid + template + class grid_ref : ref_common + { + using parent_t = ref_common; + + public: + // ctor + grid_ref(const int &n_cell, const int &n_cell_ref): + parent_t(n_cell != n_cell_ref) + { + parent_t::resize(n_cell, n_cell_ref); + } + }; + + // for arrays of the size of the number of particles + template + class part_ref + { + using parent_t = ref_common; + + public: + // ctor + part_ref(const int &n_ref): + parent_t(n_ref > 1) + {} + + void resize(int n) + { + parent_t::resize(n, n); + } + void reserve(int n) + { + parent_t::reserve(n, n); + } + }; + }; +}; diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index adca4095f..b11e619b1 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -26,8 +26,7 @@ namespace libcloudphxx namespace detail { enum { invalid = -1 }; - - }; + }; // pimpl stuff template @@ -107,9 +106,20 @@ namespace libcloudphxx // housekeeping data (per particle) thrust_device::vector - i, j, k, ijk, // Eulerian grid cell indices (always zero for 0D) + i, j, k, // Eulerian grid cell indices (always zero for 0D) ijk_ref, // ijk in refined cells sorted_id, sorted_ijk; + part_ref ijk; // ijk in the normal and in the refined grid + + // helpers that reference refined arrays if refinement is done and non-refined otherwise + // used to conserve memory if n_ref==1 (no refinement) + // needed only for arrays that have both a refined and a normal version, not needed for arrays only on the refined grid +// thrust_device::vector &ijk_ref_hlpr; + /* + , &count_ijk_ref_hlpr; + thrust_device::vector &count_num_ref_hlpr; + thrust_device::vector &count_mom_ref_hlpr, &eta_ref_hlpr, &rhod_ref_hlpr; + */ // Arakawa-C grid helper vars thrust_device::vector @@ -186,13 +196,6 @@ namespace libcloudphxx thrust::host_vector > l2e; - // helpers that reference refined arrays if refinement is done and non-refined otherwise - // used to conserve memory if n_ref==1 (no refinement) - // needed only for arrays that have both a refined and a normal version, not needed for arrays only on the refined grid - thrust_device::vector &ijk_ref_hlpr, &count_ijk_ref_hlpr; - thrust_device::vector &count_num_ref_hlpr; - thrust_device::vector &count_mom_ref_hlpr, &eta_ref_hlpr, &rhod_ref_hlpr; - // chem stuff // TODO: consider changing the unit to AMU or alike (very small numbers!) std::vector::iterator > @@ -357,12 +360,7 @@ namespace libcloudphxx allow_sstp_cond(_opts_init.sstp_cond > 1 || _opts_init.variable_dt_switch), allow_sstp_chem(_opts_init.sstp_chem > 1 || _opts_init.variable_dt_switch), pure_const_multi (((_opts_init.sd_conc) == 0) && (_opts_init.sd_const_multi > 0 || _opts_init.dry_sizes.size() > 0)), // coal prob can be greater than one only in sd_conc simulations - ijk_ref_hlpr (opts_init.n_ref > 1 ? ijk_ref : ijk), - count_ijk_ref_hlpr(opts_init.n_ref > 1 ? count_ijk_ref : count_ijk), - count_num_ref_hlpr(opts_init.n_ref > 1 ? count_num_ref : count_num), - count_mom_ref_hlpr(opts_init.n_ref > 1 ? count_mom_ref : count_mom), - eta_ref_hlpr (opts_init.n_ref > 1 ? eta_ref : eta), - rhod_ref_hlpr (opts_init.n_ref > 1 ? rhod_ref : rhod) + ijk(opts_init.n_ref) { // set 0 dev_count to mark that its not a multi_CUDA spawn diff --git a/src/particles.tpp b/src/particles.tpp index 477f7a295..ad1327f9f 100644 --- a/src/particles.tpp +++ b/src/particles.tpp @@ -28,6 +28,7 @@ #include "detail/kernel_interpolation.hpp" #include "detail/functors_host.hpp" #include "detail/ran_with_mpi.hpp" +#include "detail/grid_refinement_helpers.hpp" //kernel definitions #include "detail/kernel_definitions/hall_efficiencies.hpp" From 21a7819addf4f3cd76c9bd4cc71c428010024ff5 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 8 Jul 2022 15:52:28 +0200 Subject: [PATCH 14/38] ijk as part_ref class --- src/impl/particles_impl.ipp | 1 - src/impl/particles_impl_chem_dissoc.ipp | 4 +-- src/impl/particles_impl_chem_henry.ipp | 2 +- src/impl/particles_impl_chem_react.ipp | 2 +- src/impl/particles_impl_chem_strength.ipp | 4 +-- src/impl/particles_impl_cond.ipp | 6 ++--- src/impl/particles_impl_hskpng_ijk.ipp | 16 ++++++------ src/impl/particles_impl_hskpng_remove.ipp | 25 ++++++++++++++++--- src/impl/particles_impl_hskpng_resize.ipp | 4 +-- src/impl/particles_impl_hskpng_vterm.ipp | 4 +-- src/impl/particles_impl_init_wet.ipp | 4 +-- .../particles_impl_reserve_hskpng_npart.ipp | 1 - src/impl/particles_impl_rlx_dry_distros.ipp | 1 - .../particles_impl_update_incloud_time.ipp | 2 +- src/particles_diag.ipp | 10 ++++---- 15 files changed, 51 insertions(+), 35 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index b11e619b1..47f048582 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -107,7 +107,6 @@ namespace libcloudphxx // housekeeping data (per particle) thrust_device::vector i, j, k, // Eulerian grid cell indices (always zero for 0D) - ijk_ref, // ijk in refined cells sorted_id, sorted_ijk; part_ref ijk; // ijk in the normal and in the refined grid diff --git a/src/impl/particles_impl_chem_dissoc.ipp b/src/impl/particles_impl_chem_dissoc.ipp index 1da66ec9b..74ebbf7d5 100644 --- a/src/impl/particles_impl_chem_dissoc.ipp +++ b/src/impl/particles_impl_chem_dissoc.ipp @@ -172,11 +172,11 @@ namespace libcloudphxx zip_it_t(thrust::make_tuple( chem_bgn[SO2], chem_bgn[CO2], chem_bgn[HNO3], chem_bgn[NH3], chem_bgn[S_VI], V.begin(), - thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()))), // input - begin + thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()))), // input - begin zip_it_t(thrust::make_tuple( chem_end[SO2], chem_end[CO2], chem_end[HNO3], chem_end[NH3], chem_end[S_VI], V.end(), - thrust::make_permutation_iterator(T_ref.end(), ijk_ref_hlpr.end()))), // input - end + thrust::make_permutation_iterator(T_ref.end(), ijk.end_ref()))), // input - end chem_flag.begin(), // stencil chem_bgn[H], // output detail::chem_electroneutral(), // op diff --git a/src/impl/particles_impl_chem_henry.ipp b/src/impl/particles_impl_chem_henry.ipp index a496d1b72..65ee5f4bb 100644 --- a/src/impl/particles_impl_chem_henry.ipp +++ b/src/impl/particles_impl_chem_henry.ipp @@ -335,7 +335,7 @@ namespace libcloudphxx V.begin(), V.end(), // input - 1st arg thrust::make_zip_iterator(thrust::make_tuple( // input - 2nd arg thrust::make_permutation_iterator(p.begin(), ijk.begin()), - thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()), + thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()), thrust::make_permutation_iterator(ambient_chem[(chem_species_t)i].begin(), ijk.begin()), chem_bgn[i], rw2.begin(), diff --git a/src/impl/particles_impl_chem_react.ipp b/src/impl/particles_impl_chem_react.ipp index c48341d88..f79f7fc76 100644 --- a/src/impl/particles_impl_chem_react.ipp +++ b/src/impl/particles_impl_chem_react.ipp @@ -278,7 +278,7 @@ namespace libcloudphxx detail::chem_rhs( dt, V, - thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()), + thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()), chem_bgn[H], chem_flag ), // TODO: make it an impl member field diff --git a/src/impl/particles_impl_chem_strength.ipp b/src/impl/particles_impl_chem_strength.ipp index 3d8972611..e485c90c0 100644 --- a/src/impl/particles_impl_chem_strength.ipp +++ b/src/impl/particles_impl_chem_strength.ipp @@ -97,10 +97,10 @@ namespace libcloudphxx thrust::transform( zip_it_t(thrust::make_tuple( chem_bgn[SO2], chem_bgn[CO2], chem_bgn[HNO3], chem_bgn[NH3], chem_bgn[S_VI], chem_bgn[H], - V.begin(), thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()), rw2.begin())), // input - begin + V.begin(), thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()), rw2.begin())), // input - begin zip_it_t(thrust::make_tuple( chem_end[SO2], chem_end[CO2], chem_end[HNO3], chem_end[NH3], chem_end[S_VI], chem_end[H], - V.end(), thrust::make_permutation_iterator(T_ref.end(), ijk_ref_hlpr.end()), rw2.end())), // input - end + V.end(), thrust::make_permutation_iterator(T_ref.end(), ijk.end_ref()), rw2.end())), // input - end chem_flag.begin(), // output detail::set_chem_flag() // op ); diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index 495c6926e..32a632246 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -43,7 +43,7 @@ namespace libcloudphxx 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_ref.begin(), ijk_ref_hlpr.begin()), + thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()), thrust::make_permutation_iterator(eta.begin(), ijk.begin()), rd3.begin(), kpa.begin(), @@ -59,7 +59,7 @@ namespace libcloudphxx thrust_device::vector &RH_plus_ssp(tmp_device_real_part2); thrust::transform( ssp.begin(), ssp.end(), - thrust::make_permutation_iterator(RH_ref.begin(), ijk_ref_hlpr.begin()), + thrust::make_permutation_iterator(RH_ref.begin(), ijk.begin_ref()), RH_plus_ssp.begin(), arg::_1 + arg::_2 ); @@ -84,7 +84,7 @@ namespace libcloudphxx thrust::make_tuple( hlpr_zip_iter, thrust::make_permutation_iterator(p.begin(), ijk.begin()), - thrust::make_permutation_iterator(RH_ref.begin(), ijk_ref_hlpr.begin()) + thrust::make_permutation_iterator(RH_ref.begin(), ijk.begin_ref()) ) ), rw2.begin(), // output diff --git a/src/impl/particles_impl_hskpng_ijk.ipp b/src/impl/particles_impl_hskpng_ijk.ipp index 3b18e9cc1..081d4ec2d 100644 --- a/src/impl/particles_impl_hskpng_ijk.ipp +++ b/src/impl/particles_impl_hskpng_ijk.ipp @@ -32,20 +32,20 @@ namespace libcloudphxx template void particles_t::impl::ravel_ijk(const thrust_size_t begin_shift, const bool refined) { - auto & _ijk = refined ? ijk_ref : ijk; + auto &ijk_begin = refined ? ijk.begin_ref() : ijk.begin(); switch (n_dims) { case 0: break; case 1: - thrust::copy(i.begin()+begin_shift, i.end(), _ijk.begin()+begin_shift); + thrust::copy(i.begin()+begin_shift, i.end(), ijk_begin+begin_shift); break; case 2: namespace arg = thrust::placeholders; thrust::transform( i.begin()+begin_shift, i.end(), // input - first arg k.begin()+begin_shift, // input - second arg - _ijk.begin()+begin_shift, // output + ijk_begin+begin_shift, // output arg::_1 * opts_init.nz + arg::_2 // assuming z varies first ); break; @@ -54,14 +54,14 @@ namespace libcloudphxx thrust::transform( i.begin()+begin_shift, i.end(), // input - first arg j.begin()+begin_shift, // input - second arg - _ijk.begin()+begin_shift, // output + ijk_begin+begin_shift, // output arg::_1 * (opts_init.nz * opts_init.ny) + arg::_2 * opts_init.nz ); thrust::transform( - _ijk.begin()+begin_shift, _ijk.end(), - k.begin()+begin_shift, - _ijk.begin()+begin_shift, // in-place! + k.begin()+begin_shift, k.end() + ijk_begin+begin_shift, + ijk_begin+begin_shift, // in-place! arg::_1 + arg::_2 ); // TODO: replace these two transforms with single one @@ -149,7 +149,7 @@ namespace libcloudphxx }; }; - // ijk_ref - cell number in the refined grid + // cell number in the refined grid // begin_shift > 0 - used when we want to calculate ijk_ref for a subset of all particles (e.g. when initializing added particles) template void particles_t::impl::hskpng_ijk_ref(const thrust_size_t begin_shift, const bool _unravel_ijk) diff --git a/src/impl/particles_impl_hskpng_remove.ipp b/src/impl/particles_impl_hskpng_remove.ipp index 1b783bf71..6ef943dc5 100644 --- a/src/impl/particles_impl_hskpng_remove.ipp +++ b/src/impl/particles_impl_hskpng_remove.ipp @@ -21,12 +21,12 @@ namespace libcloudphxx { // TODO: init these vctrs once per run std::set*> n_t_vctrs; - n_t_vctrs.insert(&ijk); - 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.n_ref > 1) n_t_vctrs.insert(&ijk_ref); + + std::set*> n_t_vctrs_ref; + n_t_vctrs_ref.insert(&ijk); namespace arg = thrust::placeholders; @@ -76,6 +76,25 @@ namespace libcloudphxx ); } + // remove from n_t efined vectors + for(auto vec: n_t_vctrs_ref) + { + thrust::remove_if( + vec->begin(), + vec->begin() + n_part, + n.begin(), + arg::_1 == 0 + ); + // add conditional remove to the part_ref class? + if(opts_init.n_ref > 1) + thrust::remove_if( + vec->begin_ref(), + vec->begin_ref() + n_part, + n.begin(), + arg::_1 == 0 + ); + } + // remove from n and set new n_part auto new_end = thrust::remove_if( n.begin(), diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index f8dc28a56..757b58791 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -15,7 +15,7 @@ namespace libcloudphxx } } { - thrust_device::vector *vec[] = {&ijk, &sorted_id, &sorted_ijk}; + thrust_device::vector *vec[] = {&sorted_id, &sorted_ijk}; for(int i=0; i<3; ++i) { vec[i]->resize(n_part); @@ -24,7 +24,7 @@ namespace libcloudphxx n.resize(n_part); tmp_device_n_part.resize(n_part); tmp_device_size_part.resize(n_part); - if (opts_init.n_ref > 1) ijk_ref.resize(n_part); + ijk.resize(n_part); vt.resize(n_part, detail::invalid); diff --git a/src/impl/particles_impl_hskpng_vterm.ipp b/src/impl/particles_impl_hskpng_vterm.ipp index aaae7e3e2..cf9c1ff22 100644 --- a/src/impl/particles_impl_hskpng_vterm.ipp +++ b/src/impl/particles_impl_hskpng_vterm.ipp @@ -172,7 +172,7 @@ namespace libcloudphxx thrust::transform_if( rw2.begin(), rw2.end(), // input - 1st arg zip_it_t(thrust::make_tuple( - thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()), + thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()), thrust::make_permutation_iterator(p.begin(), ijk.begin()), thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin(), ijk.begin()) @@ -219,7 +219,7 @@ namespace libcloudphxx thrust::transform( rw2.begin(), rw2.end(), // input - 1st arg zip_it_t(thrust::make_tuple( - thrust::make_permutation_iterator(T_ref.begin(), ijk_ref_hlpr.begin()), + thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()), thrust::make_permutation_iterator(p.begin(), ijk.begin()), thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin(), ijk.begin()) diff --git a/src/impl/particles_impl_init_wet.ipp b/src/impl/particles_impl_init_wet.ipp index 96be622c2..046c606aa 100644 --- a/src/impl/particles_impl_init_wet.ipp +++ b/src/impl/particles_impl_init_wet.ipp @@ -63,8 +63,8 @@ namespace libcloudphxx zip_it_t zip_it(thrust::make_tuple( rd3.begin() + n_part_old, kpa.begin() + n_part_old, - pi_t(RH_ref.begin(), ijk_ref_hlpr.begin() + n_part_old), - pi_t(T_ref.begin(), ijk_ref_hlpr.begin() + n_part_old) + pi_t(RH_ref.begin(), ijk.begin_ref() + n_part_old), + pi_t(T_ref.begin(), ijk.begin_ref() + n_part_old) )); thrust::transform( diff --git a/src/impl/particles_impl_reserve_hskpng_npart.ipp b/src/impl/particles_impl_reserve_hskpng_npart.ipp index bbe272b65..063969e7f 100644 --- a/src/impl/particles_impl_reserve_hskpng_npart.ipp +++ b/src/impl/particles_impl_reserve_hskpng_npart.ipp @@ -17,7 +17,6 @@ namespace libcloudphxx if (opts_init.ny != 0) j.reserve(opts_init.n_sd_max); // > TODO: are they needed at all? if (opts_init.nz != 0) k.reserve(opts_init.n_sd_max); // ijk.reserve(opts_init.n_sd_max); - if(opts_init.n_ref > 1) ijk_ref.reserve(n_part); if (opts_init.nx != 0) x.reserve(opts_init.n_sd_max); if (opts_init.ny != 0) y.reserve(opts_init.n_sd_max); diff --git a/src/impl/particles_impl_rlx_dry_distros.ipp b/src/impl/particles_impl_rlx_dry_distros.ipp index e5a0de7db..528b9eca7 100644 --- a/src/impl/particles_impl_rlx_dry_distros.ipp +++ b/src/impl/particles_impl_rlx_dry_distros.ipp @@ -222,7 +222,6 @@ namespace libcloudphxx if(n_dims==3) j.resize(n_part); // we dont check in i and k because relax works only in 2D and 3D rd3.resize(n_part); n.resize(n_part); - if(opts_init.n_ref > 1) ijk_ref.resize(n_part); // --- init k --- thrust_device::vector &ptr(tmp_device_size_cell); diff --git a/src/impl/particles_impl_update_incloud_time.ipp b/src/impl/particles_impl_update_incloud_time.ipp index 369f56c39..9db16eb35 100644 --- a/src/impl/particles_impl_update_incloud_time.ipp +++ b/src/impl/particles_impl_update_incloud_time.ipp @@ -45,7 +45,7 @@ namespace libcloudphxx kpa.begin(), thrust::make_permutation_iterator( T_ref.begin(), - ijk_ref_hlpr.begin() + ijk.begin_ref() ) )), // input - 2nd arg rc2.begin(), // output diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index b3f1aff61..ceef8aa6b 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -137,7 +137,7 @@ namespace libcloudphxx // T defined in all cells pimpl->count_n_ref = pimpl->n_cell_ref; - thrust::sequence(pimpl->count_ijk_ref_hlpr.begin(), pimpl->count_ijk_ref_hlpr.end()); + thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } // records relative humidity @@ -154,7 +154,7 @@ namespace libcloudphxx // RH defined in all cells pimpl->count_n_ref = pimpl->n_cell_ref; - thrust::sequence(pimpl->count_ijk_ref_hlpr.begin(), pimpl->count_ijk_ref_hlpr.end()); + thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } // records super-droplet concentration per grid cell @@ -257,11 +257,11 @@ namespace libcloudphxx pimpl->kpa.begin(), thrust::make_permutation_iterator( pimpl->T_ref.begin(), - pimpl->ijk_ref_hlpr.begin() + pimpl->ijk.begin_ref() ), thrust::make_permutation_iterator( pimpl->RH_ref.begin(), - pimpl->ijk_ref_hlpr.begin() + pimpl->ijk.begin_ref() ) )), // input - 2nd arg RH_minus_Sc.begin(), // output @@ -286,7 +286,7 @@ namespace libcloudphxx pimpl->kpa.begin(), thrust::make_permutation_iterator( pimpl->T_ref.begin(), - pimpl->ijk_ref_hlpr.begin() + pimpl->ijk.begin_ref() ) )), // input - 2nd arg rc2.begin(), // output From 0d3dce536fffff19fabc137e5ac40d921ba714b9 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 15 Jul 2022 16:47:49 +0200 Subject: [PATCH 15/38] WIP: all ref vars as X_ref class, TODO: ctors --- src/detail/grid_refinement_helpers.hpp | 62 ++++++++++++++++++++------ src/impl/particles_impl.ipp | 58 ++++++++++++------------ src/impl/particles_impl_cond.ipp | 24 +++++----- src/impl/particles_impl_hskpng_mfp.ipp | 8 ++-- src/impl/particles_impl_init_ijk.ipp | 2 +- src/particles_diag.ipp | 20 ++++----- 6 files changed, 103 insertions(+), 71 deletions(-) diff --git a/src/detail/grid_refinement_helpers.hpp b/src/detail/grid_refinement_helpers.hpp index 68367f15f..127fc3ff1 100644 --- a/src/detail/grid_refinement_helpers.hpp +++ b/src/detail/grid_refinement_helpers.hpp @@ -11,13 +11,47 @@ namespace libcloudphxx { private: const bool ref_flag; // true if refinement is actually done - thrust_device::vector arr, arr_ref; // actual data, arr_ref initialized only if refinement is actually done protected: // ctor ref_common(const bool ref_flag): ref_flag(ref_flag) {} + }; + + + // for single values of type T + template + class ref_val : ref_common + { + private: + using parent_t = ref_common; + T val, val_ref; + + public: + using parent_t::parent_t; + + T& val() + { + return val; + } + + T& val_ref() + { + //return ref_flag ? val_ref : val; + return val_ref; + } + }; + + template + class ref_arr_common : ref_common + { + private: + using parent_t = ref_common; + thrust_device::vector arr, arr_ref; // actual data, arr_ref initialized only if refinement is actually done + + protected: + using parent_t::parent_t; auto begin() { @@ -35,27 +69,27 @@ namespace libcloudphxx { return ref_flag ? arr_ref.end() : arr.end(); } - void resize(thrust_size_t n, thrust_size_t n_ref) + void resize(thrust_size_t size, thrust_size_t size_ref) { - arr.resize(n); - if(ref_flag) arr_ref.resize(n_ref); + arr.resize(size); + if(ref_flag) arr_ref.resize(size_ref); } - void reserve(thrust_size_t n, thrust_size_t n_ref) + void reserve(thrust_size_t size, thrust_size_t size_ref) { - arr.reserve(n); - if(ref_flag) arr_ref.reserve(n_ref); + arr.reserve(size); + if(ref_flag) arr_ref.reserve(size_ref); } - } + }; // for arrays of the size of the grid template - class grid_ref : ref_common + class ref_grid : ref_arr_common { - using parent_t = ref_common; + using parent_t = ref_arr_common; public: // ctor - grid_ref(const int &n_cell, const int &n_cell_ref): + ref_grid(const int &n_cell, const int &n_cell_ref): parent_t(n_cell != n_cell_ref) { parent_t::resize(n_cell, n_cell_ref); @@ -64,13 +98,13 @@ namespace libcloudphxx // for arrays of the size of the number of particles template - class part_ref + class ref_part : ref_arr_common { - using parent_t = ref_common; + using parent_t = ref_arr_common; public: // ctor - part_ref(const int &n_ref): + ref_part(const int &n_ref): parent_t(n_ref > 1) {} diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 47f048582..3f0349c74 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -45,7 +45,7 @@ namespace libcloudphxx opts_init_t opts_init; // a copy const int n_dims; const int nx_ref, ny_ref, nz_ref; // number of refined cells in each direction - const thrust_size_t n_cell, n_cell_ref; // total number of normal and refined cells + const ref_val n_cell; // total number of normal and refined cells thrust_size_t n_part, // total number of SDs n_part_old, // total number of SDs before source n_part_to_init; // number of SDs to be initialized by source @@ -108,7 +108,7 @@ namespace libcloudphxx thrust_device::vector i, j, k, // Eulerian grid cell indices (always zero for 0D) sorted_id, sorted_ijk; - part_ref ijk; // ijk in the normal and in the refined grid + ref_part ijk; // ijk in the normal and in the refined grid // helpers that reference refined arrays if refinement is done and non-refined otherwise // used to conserve memory if n_ref==1 (no refinement) @@ -125,26 +125,27 @@ namespace libcloudphxx lft, rgt, abv, blw, fre, hnd; // TODO: could be reused after advection! // moment-counting stuff on normal and refined grid - thrust_device::vector - count_ijk, count_ijk_ref; // key-value pair for sorting particles by cell index - thrust_device::vector - count_num, count_num_ref; // number of particles in a given grid cell - thrust_device::vector - count_mom, count_mom_ref; // statistical moment // TODO (perhaps tmp_device_real_cell could be referenced?) - thrust_size_t count_n, count_n_ref; + ref_grid + count_ijk; // key-value pair for sorting particles by cell index + ref_grid + count_num; // number of particles in a given grid cell + ref_grid + count_mom; // statistical moment // TODO (perhaps tmp_device_real_cell could be referenced?) + ref_val count_n; // Eulerian-Lagrangian interface vars - thrust_device::vector + ref_grid rhod, // dry air density - rhod_ref, // same on the refined grid - th_ref, // potential temperature (dry) - rv_ref, // water vapour mixing ratio - sstp_tmp_chem_0, // ditto for trace gases - sstp_tmp_chem_1, // ditto for trace gases - sstp_tmp_chem_2, // ditto for trace gases - sstp_tmp_chem_3, // ditto for trace gases - sstp_tmp_chem_4, // ditto for trace gases - sstp_tmp_chem_5, // ditto for trace gases + th, // potential temperature (dry) + rv; // water vapour mixing ratio + + thrust_device::vector + sstp_tmp_chem_0, // trace gases + sstp_tmp_chem_1, // trace gases + sstp_tmp_chem_2, // trace gases + sstp_tmp_chem_3, // trace gases + sstp_tmp_chem_4, // trace gases + sstp_tmp_chem_5, // trace gases courant_x, courant_y, courant_z; @@ -154,13 +155,11 @@ namespace libcloudphxx // map of the accumulated volume/volume/mass of water/dry/chem that fell out of the domain std::map output_puddle; - // _ref means that the array is on the refined grid (which is the same as the normal grid if n_ref==1), some only need to be on one of the grids, some need to be on both - thrust_device::vector - T_ref, // temperature [K] - p_ref, // pressure [Pa] - RH_ref, // relative humisity - //eta,// dynamic viscosity, commented to show that we dont know how to calculate it! require input of non-refined tht? average from refined cells? - eta_ref, // dynamic viscosity + grid_ref + T, // temperature [K] + p, // pressure [Pa] + RH, // relative humisity + eta, // dynamic viscosity diss_rate; // turbulent kinetic energy dissipation rate thrust_device::vector w_LS; // large-scale subsidence velocity profile @@ -231,11 +230,10 @@ namespace libcloudphxx tmp_device_real_part4, tmp_device_real_part5, tmp_device_real_cell, - tmp_device_real_cell1, - tmp_device_real_cell2, - tmp_device_real_cell_ref, - tmp_device_real_cell_ref1, &u01; // uniform random numbers between 0 and 1 // TODO: use the tmp array as rand argument? + ref_grid + tmp_device_real_cell1, + tmp_device_real_cell2; thrust_device::vector tmp_device_n_part, &un; // uniform natural random numbers between 0 and max value of unsigned int diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index 32a632246..e27e80e35 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -19,8 +19,8 @@ namespace libcloudphxx namespace arg = thrust::placeholders; // needs to be the same as in hskpng_mfp - thrust_device::vector &lambda_D(opts_init.n_ref > 1 ? tmp_device_real_cell_ref : tmp_device_real_cell1); - thrust_device::vector &lambda_K(opts_init.n_ref > 1 ? tmp_device_real_cell_ref1 : tmp_device_real_cell2); + auto &lambda_D(tmp_device_real_cell1); + auto &lambda_K(tmp_device_real_cell2); // --- calc liquid water content before cond --- hskpng_sort(); @@ -41,15 +41,15 @@ namespace libcloudphxx ); 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_ref.begin(), ijk.begin_ref()), - thrust::make_permutation_iterator(eta.begin(), ijk.begin()), + thrust::make_permutation_iterator(rhod.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(rv.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(T.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(eta.begin_ref(), ijk.begin_ref()), rd3.begin(), kpa.begin(), vt.begin(), - thrust::make_permutation_iterator(lambda_D.begin(), ijk.begin()), - thrust::make_permutation_iterator(lambda_K.begin(), ijk.begin()) + thrust::make_permutation_iterator(lambda_D.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(lambda_K.begin_ref(), ijk.begin_ref()) )); // calculating drop growth in a timestep using backward Euler @@ -59,7 +59,7 @@ namespace libcloudphxx thrust_device::vector &RH_plus_ssp(tmp_device_real_part2); thrust::transform( ssp.begin(), ssp.end(), - thrust::make_permutation_iterator(RH_ref.begin(), ijk.begin_ref()), + thrust::make_permutation_iterator(RH.begin_ref(), ijk.begin_ref()), RH_plus_ssp.begin(), arg::_1 + arg::_2 ); @@ -69,7 +69,7 @@ namespace libcloudphxx thrust::make_zip_iterator( // input - 2nd arg thrust::make_tuple( hlpr_zip_iter, - thrust::make_permutation_iterator(p.begin(), ijk.begin()), + thrust::make_permutation_iterator(p.begin_ref(), ijk.begin_ref()), RH_plus_ssp.begin() ) ), @@ -83,8 +83,8 @@ namespace libcloudphxx thrust::make_zip_iterator( // input - 2nd arg thrust::make_tuple( hlpr_zip_iter, - thrust::make_permutation_iterator(p.begin(), ijk.begin()), - thrust::make_permutation_iterator(RH_ref.begin(), ijk.begin_ref()) + thrust::make_permutation_iterator(p.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(RH.begin_ref(), ijk.begin_ref()) ) ), rw2.begin(), // output diff --git a/src/impl/particles_impl_hskpng_mfp.ipp b/src/impl/particles_impl_hskpng_mfp.ipp index b3586682b..e875436a1 100644 --- a/src/impl/particles_impl_hskpng_mfp.ipp +++ b/src/impl/particles_impl_hskpng_mfp.ipp @@ -42,11 +42,11 @@ namespace libcloudphxx void particles_t::impl::hskpng_mfp() { // tmp_device_real_cell cant be used here, it is used in cond.ipp - thrust_device::vector &lambda_D(opts_init.n_ref > 1 ? tmp_device_real_cell_ref : tmp_device_real_cell1); - thrust_device::vector &lambda_K(opts_init.n_ref > 1 ? tmp_device_real_cell_ref1 : tmp_device_real_cell2); + auto &lambda_D(tmp_device_real_cell1); + auto &lambda_K(tmp_device_real_cell2); - thrust::transform(T_ref.begin(), T_ref.end(), lambda_D.begin(), detail::common__mean_free_path__lambda_D()); - thrust::transform(T_ref.begin(), T_ref.end(), p.begin(), lambda_K.begin(), detail::common__mean_free_path__lambda_K()); + thrust::transform(T.begin_ref(), T.end_ref(), lambda_D.begin_ref(), detail::common__mean_free_path__lambda_D()); + thrust::transform(T.begin_ref(), T.end_ref(), p.begin_ref(), lambda_K.begin_ref(), detail::common__mean_free_path__lambda_K()); } }; }; diff --git a/src/impl/particles_impl_init_ijk.ipp b/src/impl/particles_impl_init_ijk.ipp index ee2282a33..1493f88b1 100644 --- a/src/impl/particles_impl_init_ijk.ipp +++ b/src/impl/particles_impl_init_ijk.ipp @@ -46,7 +46,7 @@ namespace libcloudphxx thrust::make_zip_iterator(thrust::make_tuple( count_num.end(), ptr.end(), zero + n_cell )), - detail::arbitrary_sequence(&(ijk[n_part_old])) + detail::arbitrary_sequence(&(ijk.begin() + n_part_old)) ); } }; diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index ceef8aa6b..5bc56cacf 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -130,13 +130,13 @@ namespace libcloudphxx pimpl->hskpng_Tpr(); thrust::copy( - pimpl->T_ref.begin(), - pimpl->T_ref.end(), - pimpl->count_mom_ref_hlpr.begin() + pimpl->T.begin_ref(), + pimpl->T.end_ref(), + pimpl->count_mom.begin_ref() ); // T defined in all cells - pimpl->count_n_ref = pimpl->n_cell_ref; + pimpl->count_n.val_ref() = pimpl->n_cell_ref.val_ref(); thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } @@ -147,13 +147,13 @@ namespace libcloudphxx pimpl->hskpng_Tpr(); thrust::copy( - pimpl->RH_ref.begin(), - pimpl->RH_ref.end(), - pimpl->count_mom_ref_hlpr.begin() + pimpl->RH.begin_ref(), + pimpl->RH.end_ref(), + pimpl->count_mom.begin_ref() ); // RH defined in all cells - pimpl->count_n_ref = pimpl->n_cell_ref; + pimpl->count_n.val_ref() = pimpl->n_cell_ref.val_ref(); thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } @@ -256,11 +256,11 @@ namespace libcloudphxx thrust::make_zip_iterator(make_tuple( pimpl->kpa.begin(), thrust::make_permutation_iterator( - pimpl->T_ref.begin(), + pimpl->T.begin_ref(), pimpl->ijk.begin_ref() ), thrust::make_permutation_iterator( - pimpl->RH_ref.begin(), + pimpl->RH.begin_ref(), pimpl->ijk.begin_ref() ) )), // input - 2nd arg From b16cd329db71f8b41831d252b078a1faf21d5269 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 18 Jul 2022 15:40:39 +0200 Subject: [PATCH 16/38] n_cell ref ctor --- src/detail/grid_refinement_helpers.hpp | 20 ++++++++++++++++---- src/impl/particles_impl.ipp | 16 ++++++---------- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/src/detail/grid_refinement_helpers.hpp b/src/detail/grid_refinement_helpers.hpp index 127fc3ff1..a76d2ad8b 100644 --- a/src/detail/grid_refinement_helpers.hpp +++ b/src/detail/grid_refinement_helpers.hpp @@ -29,14 +29,26 @@ namespace libcloudphxx T val, val_ref; public: - using parent_t::parent_t; + ref_val(T val, T val_ref): + val(val), + val_ref(val_ref), + parent_t(true) + {} - T& val() + T& get() { return val; } - - T& val_ref() + const T& get() + { + return val; + } + T& get_ref() + { + //return ref_flag ? val_ref : val; + return val_ref; + } + const T& get_ref() { //return ref_flag ? val_ref : val; return val_ref; diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 3f0349c74..68f618a38 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -312,11 +312,11 @@ namespace libcloudphxx _opts_init.ny/m1(_opts_init.ny) + _opts_init.nz/m1(_opts_init.nz) ), - n_cell( - m1(_opts_init.nx) * - m1(_opts_init.ny) * - m1(_opts_init.nz) - ), + nx_ref(n_dims >= 1 ? (opts_init.nx - 1) * opts_init.n_ref + 1 : 0), // TODO: with MPI this is probably not true + ny_ref(n_dims == 3 ? (opts_init.ny - 1) * opts_init.n_ref + 1 : 0), + nz_ref(n_dims >= 2 ? (opts_init.nz - 1) * opts_init.n_ref + 1 : 0), + n_cell(m1(_opts_init.nx) * m1(_opts_init.ny) * m1(_opts_init.nz), + m1(nx_ref) * m1(ny_ref) * m1(nz_ref)), // NOTE: needs to be equal to n_cell for n_ref == 1 zero(0), n_part(0), sorted(false), @@ -347,10 +347,6 @@ namespace libcloudphxx n_dims == 2 ? halo_size * (_opts_init.nz + 1): // 2D halo_size * (_opts_init.nz + 1) * _opts_init.ny // 3D ), - nx_ref(n_dims >= 1 ? (opts_init.nx - 1) * opts_init.n_ref + 1 : 0), // TODO: with MPI this is probably not true - ny_ref(n_dims == 3 ? (opts_init.ny - 1) * opts_init.n_ref + 1 : 0), - nz_ref(n_dims >= 2 ? (opts_init.nz - 1) * opts_init.n_ref + 1 : 0), - n_cell_ref(m1(nx_ref) * m1(ny_ref) * m1(nz_ref)), // NOTE: needs to be equal to n_cell for n_ref == 1 w_LS(_opts_init.w_LS), SGS_mix_len(_opts_init.SGS_mix_len), adve_scheme(_opts_init.adve_scheme), @@ -398,7 +394,7 @@ namespace libcloudphxx break; default: assert(false); } - if (n_dims != 0) assert(n_grid > n_cell); + if (n_dims != 0) assert(n_grid > n_cell.get()); tmp_host_real_grid.resize(n_grid); } From fdb5b8cb4651a6348a605bea4e0b7eeee59de4cb Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 18 Jul 2022 15:45:00 +0200 Subject: [PATCH 17/38] n_cell ref in diag() --- src/particles_diag.ipp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index 5bc56cacf..ef056f2f8 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -119,7 +119,7 @@ namespace libcloudphxx ); // p defined in all cells - pimpl->count_n = pimpl->n_cell; + pimpl->count_n = pimpl->n_cell.get(); thrust::sequence(pimpl->count_ijk.begin(), pimpl->count_ijk.end()); } @@ -136,7 +136,7 @@ namespace libcloudphxx ); // T defined in all cells - pimpl->count_n.val_ref() = pimpl->n_cell_ref.val_ref(); + pimpl->count_n.val_ref() = pimpl->n_cell.get_ref(); thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } @@ -153,7 +153,7 @@ namespace libcloudphxx ); // RH defined in all cells - pimpl->count_n.val_ref() = pimpl->n_cell_ref.val_ref(); + pimpl->count_n.val_ref() = pimpl->n_cell.get_ref(); thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } @@ -379,7 +379,7 @@ namespace libcloudphxx case 3: thrust::transform( pimpl->count_mom.begin(), //arg1 - pimpl->count_mom.begin() + pimpl->n_cell, + pimpl->count_mom.begin() + pimpl->n_cell.get(), thrust::make_zip_iterator(thrust::make_tuple( thrust::make_permutation_iterator(pimpl->courant_y.begin(), pi(pimpl->fre.begin(), thrust::make_counting_iterator(pimpl->halo_x))), // fre counts from the start of halo, but here we need only real cells thrust::make_permutation_iterator(pimpl->courant_y.begin(), pi(pimpl->hnd.begin(), thrust::make_counting_iterator(pimpl->halo_x))) @@ -390,7 +390,7 @@ namespace libcloudphxx case 2: thrust::transform( pimpl->count_mom.begin(), //arg1 - pimpl->count_mom.begin() + pimpl->n_cell, + pimpl->count_mom.begin() + pimpl->n_cell.get(), thrust::make_zip_iterator(thrust::make_tuple( thrust::make_permutation_iterator(pimpl->courant_z.begin(), pi(pimpl->blw.begin(), thrust::make_counting_iterator(pimpl->halo_x))), thrust::make_permutation_iterator(pimpl->courant_z.begin(), pi(pimpl->abv.begin(), thrust::make_counting_iterator(pimpl->halo_x))) @@ -401,7 +401,7 @@ namespace libcloudphxx case 1: thrust::transform( pimpl->count_mom.begin(), //arg1 - pimpl->count_mom.begin() + pimpl->n_cell, + pimpl->count_mom.begin() + pimpl->n_cell.get(), thrust::make_zip_iterator(thrust::make_tuple( thrust::make_permutation_iterator(pimpl->courant_x.begin(), pi(pimpl->lft.begin(), thrust::make_counting_iterator(pimpl->halo_x))), thrust::make_permutation_iterator(pimpl->courant_x.begin(), pi(pimpl->rgt.begin(), thrust::make_counting_iterator(pimpl->halo_x))) @@ -411,7 +411,7 @@ namespace libcloudphxx ); } // divergence defined in all cells - pimpl->count_n = pimpl->n_cell; + pimpl->count_n = pimpl->n_cell.get(); thrust::sequence(pimpl->count_ijk.begin(), pimpl->count_ijk.end()); } @@ -477,7 +477,7 @@ namespace libcloudphxx ); pimpl->count_n = n.first - pimpl->count_ijk.begin(); - assert(pimpl->count_n > 0 && pimpl->count_n <= pimpl->n_cell); + assert(pimpl->count_n > 0 && pimpl->count_n <= pimpl->n_cell.get()); } // computes mean chemical properties for the selected particles From f38fc4b648c2729aae017041b319d91d8acd5688 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 18 Jul 2022 15:57:34 +0200 Subject: [PATCH 18/38] n_cell ref updates --- src/impl/particles_impl_chem_henry.ipp | 2 +- src/impl/particles_impl_coal.ipp | 4 +-- src/impl/particles_impl_cond.ipp | 3 ++- src/impl/particles_impl_hskpng_count.ipp | 2 +- src/impl/particles_impl_init_count_num.ipp | 4 +-- src/impl/particles_impl_init_grid.ipp | 26 +++++++++---------- src/impl/particles_impl_init_hskpng_ncell.ipp | 26 +++++++++---------- src/impl/particles_impl_init_ijk.ipp | 2 +- src/impl/particles_impl_init_sync.ipp | 6 ++--- src/impl/particles_impl_mass_dens.ipp | 2 +- ...articles_impl_src_dry_distros_matching.ipp | 4 +-- src/impl/particles_impl_sstp_chem.ipp | 12 ++++----- src/impl/particles_impl_xchng_courants.ipp | 2 +- 13 files changed, 48 insertions(+), 47 deletions(-) diff --git a/src/impl/particles_impl_chem_henry.ipp b/src/impl/particles_impl_chem_henry.ipp index 65ee5f4bb..84bf1dc45 100644 --- a/src/impl/particles_impl_chem_henry.ipp +++ b/src/impl/particles_impl_chem_henry.ipp @@ -370,7 +370,7 @@ namespace libcloudphxx mass_new.begin() ); count_n = it_pair.first - count_ijk.begin(); - assert(count_n > 0 && count_n <= n_cell); + assert(count_n > 0 && count_n <= n_cell.get()); // apply the change to the mixing ratios of trace gases thrust::transform( diff --git a/src/impl/particles_impl_coal.ipp b/src/impl/particles_impl_coal.ipp index 205ccacb2..50e3f32cd 100644 --- a/src/impl/particles_impl_coal.ipp +++ b/src/impl/particles_impl_coal.ipp @@ -281,7 +281,7 @@ namespace libcloudphxx // 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(count_n != n_cell.get()) thrust::fill(scl.begin(), scl.end(), real_t(0.)); thrust::copy( count_mom.begin(), // input - begin @@ -295,7 +295,7 @@ namespace libcloudphxx // cumulative sum of count_num -> (i - cumsum(ijk(i))) gives droplet index in a given cell // fill with 0s if not all cells will be updated in the following copy - if(count_n!=n_cell) thrust::fill(off.begin(), off.end(), real_t(0.)); + if(count_n != n_cell.get()) thrust::fill(off.begin(), off.end(), real_t(0.)); thrust::copy( count_num.begin(), count_num.begin() + count_n, diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index e27e80e35..c3bf4c9b0 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -33,7 +33,8 @@ namespace libcloudphxx // permute-copying the result to -dm_3 // fill with 0s if not all cells will be updated in the following transform - if(count_n!=n_cell) thrust::fill(drv.begin(), drv.end(), real_t(0.)); + if(count_n != n_cell.get()) thrust::fill(drv.begin(), drv.end(), real_t(0.)); + thrust::transform( count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output diff --git a/src/impl/particles_impl_hskpng_count.ipp b/src/impl/particles_impl_hskpng_count.ipp index 2619094c1..dfe2b5188 100644 --- a/src/impl/particles_impl_hskpng_count.ipp +++ b/src/impl/particles_impl_hskpng_count.ipp @@ -28,7 +28,7 @@ namespace libcloudphxx count_num.begin() // output - values ); count_n = it_pair.first - count_ijk.begin(); - assert(count_n <= n_cell); + assert(count_n <= n_cell.get()); } }; }; diff --git a/src/impl/particles_impl_init_count_num.ipp b/src/impl/particles_impl_init_count_num.ipp index ad3dcb414..a1ad339bc 100644 --- a/src/impl/particles_impl_init_count_num.ipp +++ b/src/impl/particles_impl_init_count_num.ipp @@ -129,7 +129,7 @@ namespace libcloudphxx case 2: thrust::transform( zero, - zero + n_cell, + zero + n_cell.get(), count_num.begin(), real_t(number) // no of SDs to create * ((arg::_1 % opts_init.nz) < k1) @@ -141,7 +141,7 @@ namespace libcloudphxx case 3: thrust::transform( zero, - zero + n_cell, + zero + n_cell.get(), count_num.begin(), real_t(number) // no of SDs to create * ((arg::_1 % opts_init.nz) < k1) diff --git a/src/impl/particles_impl_init_grid.ipp b/src/impl/particles_impl_init_grid.ipp index bfdcf9319..2a9759999 100644 --- a/src/impl/particles_impl_init_grid.ipp +++ b/src/impl/particles_impl_init_grid.ipp @@ -59,7 +59,7 @@ namespace libcloudphxx namespace arg = thrust::placeholders; // filling in sample volume data - dv.resize(n_cell); + dv.resize(n_cell.get()); int n_cell_halo(2 * halo_x); // number of halo cells @@ -67,27 +67,27 @@ namespace libcloudphxx { // in parcel set-up hskpng_Tpr takes care of keeping dv up-to-date with rho (dealing with 1kg of dry air) thrust::transform( - zero, zero + n_cell, // input - 1st arg + zero, zero + n_cell.get(), // input - 1st arg dv.begin(), // output detail::dv_eval(opts_init) ); // memory allocation - lft.resize(n_cell + n_cell_halo); - rgt.resize(n_cell + n_cell_halo); + lft.resize(n_cell.get() + n_cell_halo); + rgt.resize(n_cell.get() + n_cell_halo); } if (n_dims > 1) { // memory allocation - abv.resize(n_cell + n_cell_halo); - blw.resize(n_cell + n_cell_halo); + abv.resize(n_cell.get() + n_cell_halo); + blw.resize(n_cell.get() + n_cell_halo); } if (n_dims == 3) { // memory allocation - fre.resize(n_cell + n_cell_halo); - hnd.resize(n_cell + n_cell_halo); + fre.resize(n_cell.get() + n_cell_halo); + hnd.resize(n_cell.get() + n_cell_halo); } // filling in neighbour info data @@ -96,7 +96,7 @@ namespace libcloudphxx case 0: break; case 2: thrust::transform( - zero, zero + n_cell + n_cell_halo, // input - 1st arg + zero, zero + n_cell.get() + n_cell_halo, // input - 1st arg blw.begin(), // output arg::_1 + (arg::_1 / opts_init.nz) ); @@ -108,7 +108,7 @@ namespace libcloudphxx // intentionally no break! case 1: thrust::transform( - zero, zero + n_cell + n_cell_halo, // input - 1st arg + zero, zero + n_cell.get() + n_cell_halo, // input - 1st arg lft.begin(), // output arg::_1 ); @@ -120,7 +120,7 @@ namespace libcloudphxx break; case 3: thrust::transform( - zero, zero + n_cell + n_cell_halo, // input - 1st arg + zero, zero + n_cell.get() + n_cell_halo, // input - 1st arg lft.begin(), // output arg::_1 ); @@ -130,7 +130,7 @@ namespace libcloudphxx arg::_1 + opts_init.nz * opts_init.ny ); thrust::transform( - zero, zero + n_cell + n_cell_halo, // input - 1st arg + zero, zero + n_cell.get() + n_cell_halo, // input - 1st arg blw.begin(), // output arg::_1 + opts_init.ny * (arg::_1 / (opts_init.nz * opts_init.ny)) @@ -142,7 +142,7 @@ namespace libcloudphxx arg::_1 + 1 ); thrust::transform( - zero, zero + n_cell + n_cell_halo, // input - 1st arg + zero, zero + n_cell.get() + n_cell_halo, // input - 1st arg fre.begin(), // output arg::_1 + (arg::_1 / (opts_init.nz * opts_init.ny)) * opts_init.nz ); diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index 25e84c732..fb3fff0c6 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -13,25 +13,25 @@ namespace libcloudphxx void particles_t::impl::init_hskpng_ncell() { // memory allocation, p already resized in init_sync() - eta.resize(n_cell); + eta.resize(n_cell.get()); - count_ijk.resize(n_cell); - count_num.resize(n_cell); - count_mom.resize(n_cell); + count_ijk.resize(n_cell.get()); + count_num.resize(n_cell.get()); + count_mom.resize(n_cell.get()); count_n = 0; // initialising device temporary arrays - tmp_device_real_cell.resize(n_cell); - tmp_device_real_cell1.resize(n_cell); - tmp_device_real_cell2.resize(n_cell); - tmp_device_size_cell.resize(n_cell); - tmp_host_size_cell.resize(n_cell); - tmp_host_real_cell.resize(n_cell); + tmp_device_real_cell.resize(n_cell.get()); + tmp_device_real_cell1.resize(n_cell.get()); + tmp_device_real_cell2.resize(n_cell.get()); + tmp_device_size_cell.resize(n_cell.get()); + tmp_host_size_cell.resize(n_cell.get()); + tmp_host_real_cell.resize(n_cell.get()); if(allow_sstp_cond && !opts_init.exact_sstp_cond) { - sstp_tmp_rv.resize(n_cell); - sstp_tmp_th.resize(n_cell); - sstp_tmp_rh.resize(n_cell); + sstp_tmp_rv.resize(n_cell.get()); + sstp_tmp_th.resize(n_cell.get()); + sstp_tmp_rh.resize(n_cell.get()); } // arrays on the refined grid diff --git a/src/impl/particles_impl_init_ijk.ipp b/src/impl/particles_impl_init_ijk.ipp index 1493f88b1..37bde74ad 100644 --- a/src/impl/particles_impl_init_ijk.ipp +++ b/src/impl/particles_impl_init_ijk.ipp @@ -44,7 +44,7 @@ namespace libcloudphxx count_num.begin(), ptr.begin(), zero )), thrust::make_zip_iterator(thrust::make_tuple( - count_num.end(), ptr.end(), zero + n_cell + count_num.end(), ptr.end(), zero + n_cell.get() )), detail::arbitrary_sequence(&(ijk.begin() + n_part_old)) ); diff --git a/src/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index 3c7d4c405..f53bd54fa 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/src/impl/particles_impl_init_sync.ipp @@ -13,16 +13,16 @@ namespace libcloudphxx void particles_t::impl::init_sync() { // memory allocation for scalar fields - rhod.resize(n_cell); + rhod.resize(n_cell.get()); rhod_ref.resize(n_cell_ref); p_ref.resize(n_cell_ref); th_ref.resize(n_cell_ref); rv_ref.resize(n_cell_ref); if(opts_init.chem_switch) for (int i = 0; i < chem_gas_n; ++i) - ambient_chem[(chem_species_t)i].resize(n_cell); + ambient_chem[(chem_species_t)i].resize(n_cell.get()); if(opts_init.turb_cond_switch || opts_init.turb_adve_switch || opts_init.turb_coal_switch) - diss_rate.resize(n_cell); + diss_rate.resize(n_cell.get()); // memory allocation for vector fields (Arakawa-C grid) // halo in the x dimensions (which could be distmem boundary) diff --git a/src/impl/particles_impl_mass_dens.ipp b/src/impl/particles_impl_mass_dens.ipp index 7975ce2ac..8340d83e1 100644 --- a/src/impl/particles_impl_mass_dens.ipp +++ b/src/impl/particles_impl_mass_dens.ipp @@ -81,7 +81,7 @@ namespace libcloudphxx ); count_n = it_pair.first - count_ijk.begin(); - assert(count_n > 0 && count_n <= n_cell); + assert(count_n > 0 && count_n <= n_cell.get()); #if !defined(__NVCC__) using std::sqrt; diff --git a/src/impl/particles_impl_src_dry_distros_matching.ipp b/src/impl/particles_impl_src_dry_distros_matching.ipp index f0c201fd9..c51e23cb2 100644 --- a/src/impl/particles_impl_src_dry_distros_matching.ipp +++ b/src/impl/particles_impl_src_dry_distros_matching.ipp @@ -257,9 +257,9 @@ namespace libcloudphxx } // tmp vector to hold number of particles in a given size bin in a given cell - thrust_device::vector bin_cell_count(n_part_tot_in_src + n_cell + 1); // needs space for out_of_bins + thrust_device::vector bin_cell_count(n_part_tot_in_src + n_cell.get() + 1); // needs space for out_of_bins // tmp vector for number of particles in bins up to this one - thrust_device::vector bin_cell_count_ptr(n_part_tot_in_src + n_cell + 1); + thrust_device::vector bin_cell_count_ptr(n_part_tot_in_src + n_cell.get() + 1); thrust_size_t count_bins; { diff --git a/src/impl/particles_impl_sstp_chem.ipp b/src/impl/particles_impl_sstp_chem.ipp index cfc10d2e2..b347d3cf2 100644 --- a/src/impl/particles_impl_sstp_chem.ipp +++ b/src/impl/particles_impl_sstp_chem.ipp @@ -38,12 +38,12 @@ namespace libcloudphxx if (!allow_sstp_chem) return; // memory allocation - sstp_tmp_chem_0.resize(n_cell); - sstp_tmp_chem_1.resize(n_cell); - sstp_tmp_chem_2.resize(n_cell); - sstp_tmp_chem_3.resize(n_cell); - sstp_tmp_chem_4.resize(n_cell); - sstp_tmp_chem_5.resize(n_cell); + sstp_tmp_chem_0.resize(n_cell.get()); + sstp_tmp_chem_1.resize(n_cell.get()); + sstp_tmp_chem_2.resize(n_cell.get()); + sstp_tmp_chem_3.resize(n_cell.get()); + sstp_tmp_chem_4.resize(n_cell.get()); + sstp_tmp_chem_5.resize(n_cell.get()); // initialise _old values sstp_save_chem(); diff --git a/src/impl/particles_impl_xchng_courants.ipp b/src/impl/particles_impl_xchng_courants.ipp index bf9ac8b5d..5a06f8f0a 100644 --- a/src/impl/particles_impl_xchng_courants.ipp +++ b/src/impl/particles_impl_xchng_courants.ipp @@ -29,7 +29,7 @@ namespace libcloudphxx n_dims == 2 ? (halo_size + 1) * opts_init.nz: // 2D (halo_size + 1) * opts_init.nz * opts_init.ny // 3D ); - const int cx_rgt_internal_idx = n_cell; + const int cx_rgt_internal_idx = n_cell.get(); const int cz_lft_internal_idx = halo_z; const int cz_rgt_internal_idx( From 73fa21d7bb3995a79f20b4c6347e2cd3cea41d96 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 18 Jul 2022 16:06:37 +0200 Subject: [PATCH 19/38] count_n ref update --- src/impl/particles_impl_ante_adding_SD.ipp | 2 +- src/impl/particles_impl_chem_henry.ipp | 6 ++--- src/impl/particles_impl_coal.ipp | 12 +++++----- src/impl/particles_impl_cond.ipp | 10 ++++----- src/impl/particles_impl_fill_outbuf.ipp | 2 +- src/impl/particles_impl_hskpng_count.ipp | 4 ++-- src/impl/particles_impl_hskpng_turb_ss.ipp | 4 ++-- src/impl/particles_impl_init_hskpng_ncell.ipp | 4 ++-- src/impl/particles_impl_mass_dens.ipp | 6 ++--- src/impl/particles_impl_moms.ipp | 22 +++++++++---------- src/impl/particles_impl_post_adding_SD.ipp | 2 +- src/impl/particles_impl_rlx_dry_distros.ipp | 6 ++--- src/impl/particles_impl_update_th_rv.ipp | 4 ++-- src/particles_diag.ipp | 14 ++++++------ 14 files changed, 49 insertions(+), 49 deletions(-) diff --git a/src/impl/particles_impl_ante_adding_SD.ipp b/src/impl/particles_impl_ante_adding_SD.ipp index b77502a11..e7b496a6e 100644 --- a/src/impl/particles_impl_ante_adding_SD.ipp +++ b/src/impl/particles_impl_ante_adding_SD.ipp @@ -22,7 +22,7 @@ namespace libcloudphxx // drv = - tot_vol_bfr thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg + count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output thrust::negate() ); diff --git a/src/impl/particles_impl_chem_henry.ipp b/src/impl/particles_impl_chem_henry.ipp index 84bf1dc45..df92f5a02 100644 --- a/src/impl/particles_impl_chem_henry.ipp +++ b/src/impl/particles_impl_chem_henry.ipp @@ -369,12 +369,12 @@ namespace libcloudphxx count_ijk.begin(), mass_new.begin() ); - count_n = it_pair.first - count_ijk.begin(); - assert(count_n > 0 && count_n <= n_cell.get()); + count_n.get() = it_pair.first - count_ijk.begin(); + assert(count_n.get() > 0 && count_n.get() <= n_cell.get()); // apply the change to the mixing ratios of trace gases thrust::transform( - mass_new.begin(), mass_new.begin() + count_n, // input - 1st arg + mass_new.begin(), mass_new.begin() + count_n.get(), // input - 1st arg thrust::make_zip_iterator(thrust::make_tuple( // input - 2nd arg mass_old.begin(), thrust::make_permutation_iterator(rhod.begin(), count_ijk.begin()), diff --git a/src/impl/particles_impl_coal.ipp b/src/impl/particles_impl_coal.ipp index 50e3f32cd..67e106b03 100644 --- a/src/impl/particles_impl_coal.ipp +++ b/src/impl/particles_impl_coal.ipp @@ -265,11 +265,11 @@ namespace libcloudphxx // placing scale_factors in count_mom (of size count_n!) thrust::transform( - count_num.begin(), count_num.begin() + count_n, // input - 1st arg + count_num.begin(), count_num.begin() + count_n.get(), // input - 1st arg count_mom.begin(), // output detail::scale_factor() ); - nancheck_range(count_mom.begin(), count_mom.begin() + count_n, "count_mom storing scale_factors"); + nancheck_range(count_mom.begin(), count_mom.begin() + count_n.get(), "count_mom storing scale_factors"); // references to tmp data thrust_device::vector @@ -281,11 +281,11 @@ namespace libcloudphxx // 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.get()) thrust::fill(scl.begin(), scl.end(), real_t(0.)); + if(count_n.get() != n_cell.get()) thrust::fill(scl.begin(), scl.end(), real_t(0.)); thrust::copy( count_mom.begin(), // input - begin - count_mom.begin() + count_n, // input - end + count_mom.begin() + count_n.get(), // input - end thrust::make_permutation_iterator( // output scl.begin(), // data count_ijk.begin() // permutation @@ -295,10 +295,10 @@ namespace libcloudphxx // cumulative sum of count_num -> (i - cumsum(ijk(i))) gives droplet index in a given cell // fill with 0s if not all cells will be updated in the following copy - if(count_n != n_cell.get()) thrust::fill(off.begin(), off.end(), real_t(0.)); + if(count_n.get() != n_cell.get()) thrust::fill(off.begin(), off.end(), real_t(0.)); thrust::copy( count_num.begin(), - count_num.begin() + count_n, + count_num.begin() + count_n.get(), thrust::make_permutation_iterator( // output off.begin(), // data count_ijk.begin() // permutation diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index c3bf4c9b0..59ad35308 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -29,14 +29,14 @@ namespace libcloudphxx // calculating the 3rd wet moment before condensation moms_all(); moms_calc(rw2.begin(), real_t(3./2.)); - nancheck_range(count_mom.begin(), count_mom.begin() + count_n, "count_mom (3rd wet moment) before condensation"); + nancheck_range(count_mom.begin(), count_mom.begin() + count_n.get(), "count_mom (3rd wet moment) before condensation"); // permute-copying the result to -dm_3 // fill with 0s if not all cells will be updated in the following transform - if(count_n != n_cell.get()) thrust::fill(drv.begin(), drv.end(), real_t(0.)); + if(count_n.get() != n_cell.get()) thrust::fill(drv.begin(), drv.end(), real_t(0.)); thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg + count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output thrust::negate() ); @@ -95,11 +95,11 @@ namespace libcloudphxx // calculating the 3rd wet moment after condensation moms_calc(rw2.begin(), real_t(3./2.)); - nancheck_range(count_mom.begin(), count_mom.begin() + count_n, "count_mom (3rd wet moment) after condensation"); + nancheck_range(count_mom.begin(), count_mom.begin() + count_n.get(), "count_mom (3rd wet moment) after condensation"); // adding the third moment after condensation to dm_3 thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg + count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // input - 2nd arg thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output thrust::plus() diff --git a/src/impl/particles_impl_fill_outbuf.ipp b/src/impl/particles_impl_fill_outbuf.ipp index 2855fe550..a38e2af08 100644 --- a/src/impl/particles_impl_fill_outbuf.ipp +++ b/src/impl/particles_impl_fill_outbuf.ipp @@ -29,7 +29,7 @@ namespace libcloudphxx thrust::copy( count_mom.begin(), // input - begin - count_mom.begin() + count_n, // input - end + count_mom.begin() + count_n.get(), // input - end thrust::make_permutation_iterator( // output tmp_host_real_cell.begin(), // data pi.begin() // permutation diff --git a/src/impl/particles_impl_hskpng_count.ipp b/src/impl/particles_impl_hskpng_count.ipp index dfe2b5188..c9733059e 100644 --- a/src/impl/particles_impl_hskpng_count.ipp +++ b/src/impl/particles_impl_hskpng_count.ipp @@ -27,8 +27,8 @@ namespace libcloudphxx count_ijk.begin(), // output - keys count_num.begin() // output - values ); - count_n = it_pair.first - count_ijk.begin(); - assert(count_n <= n_cell.get()); + count_n.get() = it_pair.first - count_ijk.begin(); + assert(count_n.get() <= n_cell.get()); } }; }; diff --git a/src/impl/particles_impl_hskpng_turb_ss.ipp b/src/impl/particles_impl_hskpng_turb_ss.ipp index f14d8ea74..735dc91b8 100644 --- a/src/impl/particles_impl_hskpng_turb_ss.ipp +++ b/src/impl/particles_impl_hskpng_turb_ss.ipp @@ -55,7 +55,7 @@ namespace libcloudphxx moms_calc(rw2.begin(), real_t(1./2), false); thrust::transform( count_mom.begin(), - count_mom.begin() + count_n, + count_mom.begin() + count_n.get(), thrust::make_permutation_iterator( dv.begin(), count_ijk.begin() @@ -67,7 +67,7 @@ namespace libcloudphxx // copy the calculated relaxation times to the array of length n_cell thrust::copy( count_mom.begin(), - count_mom.begin() + count_n, + count_mom.begin() + count_n.get(), thrust::make_permutation_iterator( tau_rlx.begin(), count_ijk.begin() diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index fb3fff0c6..f5e1f8f79 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -18,7 +18,8 @@ namespace libcloudphxx count_ijk.resize(n_cell.get()); count_num.resize(n_cell.get()); count_mom.resize(n_cell.get()); - count_n = 0; + count_n.get() = 0; + count_n.get_ref() = 0; // initialising device temporary arrays tmp_device_real_cell.resize(n_cell.get()); @@ -47,7 +48,6 @@ namespace libcloudphxx count_ijk_ref.resize(n_cell_ref); count_num_ref.resize(n_cell_ref); count_mom_ref.resize(n_cell_ref); - count_n_ref = 0; } } }; diff --git a/src/impl/particles_impl_mass_dens.ipp b/src/impl/particles_impl_mass_dens.ipp index 8340d83e1..bf3736fff 100644 --- a/src/impl/particles_impl_mass_dens.ipp +++ b/src/impl/particles_impl_mass_dens.ipp @@ -80,8 +80,8 @@ namespace libcloudphxx count_mom.begin() ); - count_n = it_pair.first - count_ijk.begin(); - assert(count_n > 0 && count_n <= n_cell.get()); + count_n.get() = it_pair.first - count_ijk.begin(); + assert(count_n.get() > 0 && count_n.get() <= n_cell.get()); #if !defined(__NVCC__) using std::sqrt; @@ -97,7 +97,7 @@ namespace libcloudphxx //multiply by prefactor and divide by dv thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - first arg + count_mom.begin(), count_mom.begin() + count_n.get(), // input - first arg thrust::make_permutation_iterator( // input - second arg dv.begin(), count_ijk.begin() diff --git a/src/impl/particles_impl_moms.ipp b/src/impl/particles_impl_moms.ipp index 4baac9cd8..aad203795 100644 --- a/src/impl/particles_impl_moms.ipp +++ b/src/impl/particles_impl_moms.ipp @@ -226,22 +226,22 @@ namespace libcloudphxx count_mom.begin() ); - count_n = it_pair.first - count_ijk.begin(); + count_n.get() = it_pair.first - count_ijk.begin(); #if !defined(NDEBUG) { - int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n, isnaninf(), 0, thrust::plus()); + int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n.get(), isnaninf(), 0, thrust::plus()); if(nan_count>0) { std::cout << nan_count << " nan/inf numbers detected in count_mom after reduce_by_key " << std::endl; } } #endif - assert(count_n <= n_cell); + assert(count_n.get() <= n_cell.get()); if(specific) { // dividing by dv thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - first arg + count_mom.begin(), count_mom.begin() + count_n.get(), // input - first arg thrust::make_permutation_iterator( // input - second arg dv.begin(), count_ijk.begin() @@ -251,7 +251,7 @@ namespace libcloudphxx ); #if !defined(NDEBUG) { - int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n, isnaninf(), 0, thrust::plus()); + int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n.get(), isnaninf(), 0, thrust::plus()); if(nan_count>0) { std::cout << nan_count << " nan/inf numbers detected in count_mom after dividing by dv " << std::endl; @@ -261,7 +261,7 @@ namespace libcloudphxx // dividing by rhod to get specific moments // (for compatibility with blk_1m and blk_2m reporting mixing ratios) thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - first arg + count_mom.begin(), count_mom.begin() + count_n.get(), // input - first arg thrust::make_permutation_iterator( // input - second arg rhod.begin(), count_ijk.begin() @@ -271,7 +271,7 @@ namespace libcloudphxx ); #if !defined(NDEBUG) { - int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n, isnaninf(), 0, thrust::plus()); + int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n.get(), isnaninf(), 0, thrust::plus()); if(nan_count>0) { std::cout << nan_count << " nan/inf numbers detected in count_mom after dividing by rhod " << std::endl; @@ -281,15 +281,15 @@ namespace libcloudphxx } #if !defined(NDEBUG) { - int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n, isnaninf(), 0, thrust::plus()); + int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n.get(), isnaninf(), 0, thrust::plus()); if(nan_count>0) { std::cout << nan_count << " nan/inf numbers detected in count_mom " << std::endl; - std::cout << "count_n:" << count_n << std::endl; + std::cout << "count_n:" << count_n.get() << std::endl; std::cout << "count_mom:" << std::endl; - debug::print(count_mom.begin(), count_mom.begin() + count_n); + debug::print(count_mom.begin(), count_mom.begin() + count_n.get()); std::cout << "count_ijk:" << std::endl; - debug::print(count_ijk.begin(), count_ijk.begin() + count_n); + debug::print(count_ijk.begin(), count_ijk.begin() + count_n.get()); std::cout << "n_filtered:" << std::endl; debug::print(n_filtered); std::cout << "sorted_ijk:" << std::endl; diff --git a/src/impl/particles_impl_post_adding_SD.ipp b/src/impl/particles_impl_post_adding_SD.ipp index 9bebe5352..675c58485 100644 --- a/src/impl/particles_impl_post_adding_SD.ipp +++ b/src/impl/particles_impl_post_adding_SD.ipp @@ -24,7 +24,7 @@ namespace libcloudphxx // drv = tot_vol_after -tot_vol_bfr + dry_vol_bfr thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg + count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // 2nd arg thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output thrust::plus() diff --git a/src/impl/particles_impl_rlx_dry_distros.ipp b/src/impl/particles_impl_rlx_dry_distros.ipp index 528b9eca7..191770879 100644 --- a/src/impl/particles_impl_rlx_dry_distros.ipp +++ b/src/impl/particles_impl_rlx_dry_distros.ipp @@ -166,10 +166,10 @@ namespace libcloudphxx // horizontal sum of this moment thrust::fill(hor_sum.begin(), hor_sum.end(), 0); thrust_device::vector &count_k(tmp_device_size_cell); // NOTE: tmp_device_size_cell is also used in some other inits, careful not to overwrite it! - thrust::transform(count_ijk.begin(), count_ijk.begin() + count_n, count_k.begin(), arg::_1 % opts_init.nz); - thrust::sort_by_key(count_k.begin(), count_k.begin() + count_n, count_mom.begin()); + thrust::transform(count_ijk.begin(), count_ijk.begin() + count_n.get(), count_k.begin(), arg::_1 % opts_init.nz); + thrust::sort_by_key(count_k.begin(), count_k.begin() + count_n.get(), count_mom.begin()); - auto new_end = thrust::reduce_by_key(count_k.begin(), count_k.begin() + count_n, count_mom.begin(), hor_sum_k.begin(), hor_sum_count.begin()); + auto new_end = thrust::reduce_by_key(count_k.begin(), count_k.begin() + count_n.get(), count_mom.begin(), hor_sum_k.begin(), hor_sum_count.begin()); int number_of_levels_with_droplets = new_end.first - hor_sum_k.begin(); // number of levels with any SD, not with SD in this size and kappa range diff --git a/src/impl/particles_impl_update_th_rv.ipp b/src/impl/particles_impl_update_th_rv.ipp index af575e038..7c65227ce 100644 --- a/src/impl/particles_impl_update_th_rv.ipp +++ b/src/impl/particles_impl_update_th_rv.ipp @@ -113,11 +113,11 @@ namespace libcloudphxx count_ijk.begin(), count_mom.begin() ); - count_n = it_pair.first - count_ijk.begin(); + count_n.get() = it_pair.first - count_ijk.begin(); // add this sum to dstate thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg + count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg thrust::make_permutation_iterator(dstate.begin(), count_ijk.begin()), // 2nd arg thrust::make_permutation_iterator(dstate.begin(), count_ijk.begin()), // output thrust::plus() diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index ef056f2f8..7905506bc 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -119,7 +119,7 @@ namespace libcloudphxx ); // p defined in all cells - pimpl->count_n = pimpl->n_cell.get(); + pimpl->count_n.get() = pimpl->n_cell.get(); thrust::sequence(pimpl->count_ijk.begin(), pimpl->count_ijk.end()); } @@ -136,7 +136,7 @@ namespace libcloudphxx ); // T defined in all cells - pimpl->count_n.val_ref() = pimpl->n_cell.get_ref(); + pimpl->count_n.get_ref() = pimpl->n_cell.get_ref(); thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } @@ -153,7 +153,7 @@ namespace libcloudphxx ); // RH defined in all cells - pimpl->count_n.val_ref() = pimpl->n_cell.get_ref(); + pimpl->count_n.get_ref() = pimpl->n_cell.get_ref(); thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } @@ -179,7 +179,7 @@ namespace libcloudphxx pimpl->count_ijk.begin(), // output - keys pimpl->count_mom.begin() // output - values ); - pimpl->count_n = n.first - pimpl->count_ijk.begin(); + pimpl->count_n.get() = n.first - pimpl->count_ijk.begin(); } // selected all particles @@ -411,7 +411,7 @@ namespace libcloudphxx ); } // divergence defined in all cells - pimpl->count_n = pimpl->n_cell.get(); + pimpl->count_n.get() = pimpl->n_cell.get(); thrust::sequence(pimpl->count_ijk.begin(), pimpl->count_ijk.end()); } @@ -476,8 +476,8 @@ namespace libcloudphxx thrust::maximum() ); - pimpl->count_n = n.first - pimpl->count_ijk.begin(); - assert(pimpl->count_n > 0 && pimpl->count_n <= pimpl->n_cell.get()); + pimpl->count_n.get() = n.first - pimpl->count_ijk.begin(); + assert(pimpl->count_n.get() > 0 && pimpl->count_n.get() <= pimpl->n_cell.get()); } // computes mean chemical properties for the selected particles From 7a752c3f5905d9cc07ba61b398e6314f337d38eb Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 20 Jul 2022 19:01:50 +0200 Subject: [PATCH 20/38] _ref WIP --- src/detail/grid_refinement_helpers.hpp | 53 ++++++++++++++++--- src/impl/particles_impl.ipp | 12 ++--- src/impl/particles_impl_chem_dissoc.ipp | 4 +- src/impl/particles_impl_chem_henry.ipp | 2 +- src/impl/particles_impl_chem_react.ipp | 2 +- src/impl/particles_impl_chem_strength.ipp | 4 +- src/impl/particles_impl_hskpng_Tpr.ipp | 28 +++++----- src/impl/particles_impl_hskpng_ijk.ipp | 2 +- src/impl/particles_impl_hskpng_remove.ipp | 2 +- src/impl/particles_impl_hskpng_vterm.ipp | 16 +++--- src/impl/particles_impl_init_hskpng_ncell.ipp | 33 ++++-------- src/impl/particles_impl_init_sync.ipp | 9 ++-- src/impl/particles_impl_init_wet.ipp | 4 +- .../particles_impl_update_incloud_time.ipp | 2 +- src/impl/particles_impl_update_th_rv.ipp | 2 +- src/particles_init.ipp | 17 +++--- 16 files changed, 109 insertions(+), 83 deletions(-) diff --git a/src/detail/grid_refinement_helpers.hpp b/src/detail/grid_refinement_helpers.hpp index a76d2ad8b..7d03797a3 100644 --- a/src/detail/grid_refinement_helpers.hpp +++ b/src/detail/grid_refinement_helpers.hpp @@ -39,20 +39,24 @@ namespace libcloudphxx { return val; } + /* const T& get() { return val; } + */ T& get_ref() { //return ref_flag ? val_ref : val; return val_ref; } + /* const T& get_ref() { //return ref_flag ? val_ref : val; return val_ref; } + */ }; template @@ -60,7 +64,8 @@ namespace libcloudphxx { private: using parent_t = ref_common; - thrust_device::vector arr, arr_ref; // actual data, arr_ref initialized only if refinement is actually done + using vec_t = thrust_device::vector; + vec_t arr, arr_ref; // actual data, arr_ref initialized only if refinement is actually done protected: using parent_t::parent_t; @@ -71,7 +76,7 @@ namespace libcloudphxx } auto begin_ref() { - return ref_flag ? arr_ref.begin() : arr.begin(); + return this->ref_flag ? arr_ref.begin() : arr.begin(); } auto end() { @@ -79,23 +84,55 @@ namespace libcloudphxx } auto end_ref() { - return ref_flag ? arr_ref.end() : arr.end(); + return this->ref_flag ? arr_ref.end() : arr.end(); } void resize(thrust_size_t size, thrust_size_t size_ref) { arr.resize(size); - if(ref_flag) arr_ref.resize(size_ref); + if(this->ref_flag) arr_ref.resize(size_ref); + } + void resize(ref_val size) + { + arr.resize(size.get()); + if(this->ref_flag) arr_ref.resize(size.get_ref()); } void reserve(thrust_size_t size, thrust_size_t size_ref) { arr.reserve(size); - if(ref_flag) arr_ref.reserve(size_ref); + if(this->ref_flag) arr_ref.reserve(size_ref); + } + void reserve(ref_val size) + { + arr.reserve(size.get()); + if(this->ref_flag) arr_ref.reserve(size.get_ref()); + } + + public: + + auto ptr() + { + return &arr; + } + + auto ptr_ref() + { + return &arr_ref; + } + + vec_t& get() + { + return arr; + } + + vec_t& get_ref() + { + return arr_ref; } }; // for arrays of the size of the grid template - class ref_grid : ref_arr_common + class ref_grid : public ref_arr_common { using parent_t = ref_arr_common; @@ -106,6 +143,10 @@ namespace libcloudphxx { parent_t::resize(n_cell, n_cell_ref); } + + ref_grid(ref_val n): + ref_grid(n.get(), n.get_ref()) + {} }; // for arrays of the size of the number of particles diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 68f618a38..ac1e73c59 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -45,7 +45,7 @@ namespace libcloudphxx opts_init_t opts_init; // a copy const int n_dims; const int nx_ref, ny_ref, nz_ref; // number of refined cells in each direction - const ref_val n_cell; // total number of normal and refined cells + ref_val n_cell; // total number of normal and refined cells thrust_size_t n_part, // total number of SDs n_part_old, // total number of SDs before source n_part_to_init; // number of SDs to be initialized by source @@ -85,9 +85,6 @@ namespace libcloudphxx wp, // turbulent perturbation of velocity ssp, // turbulent perturbation of supersaturation dot_ssp, // time derivative of the turbulent perturbation of supersaturation - sstp_tmp_rv, // either rv_old or advection-caused change in water vapour mixing ratio - sstp_tmp_th, // ditto for theta - sstp_tmp_rh, // ditto for rho sstp_tmp_p, // ditto for pressure incloud_time; // time this SD has been within a cloud @@ -137,7 +134,10 @@ namespace libcloudphxx ref_grid rhod, // dry air density th, // potential temperature (dry) - rv; // water vapour mixing ratio + rv, // water vapour mixing ratio + 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 thrust_device::vector sstp_tmp_chem_0, // trace gases @@ -155,7 +155,7 @@ namespace libcloudphxx // map of the accumulated volume/volume/mass of water/dry/chem that fell out of the domain std::map output_puddle; - grid_ref + ref_grid T, // temperature [K] p, // pressure [Pa] RH, // relative humisity diff --git a/src/impl/particles_impl_chem_dissoc.ipp b/src/impl/particles_impl_chem_dissoc.ipp index 74ebbf7d5..cfa0a691b 100644 --- a/src/impl/particles_impl_chem_dissoc.ipp +++ b/src/impl/particles_impl_chem_dissoc.ipp @@ -172,11 +172,11 @@ namespace libcloudphxx zip_it_t(thrust::make_tuple( chem_bgn[SO2], chem_bgn[CO2], chem_bgn[HNO3], chem_bgn[NH3], chem_bgn[S_VI], V.begin(), - thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()))), // input - begin + thrust::make_permutation_iterator(T.begin(), ijk.begin()))), // input - begin zip_it_t(thrust::make_tuple( chem_end[SO2], chem_end[CO2], chem_end[HNO3], chem_end[NH3], chem_end[S_VI], V.end(), - thrust::make_permutation_iterator(T_ref.end(), ijk.end_ref()))), // input - end + thrust::make_permutation_iterator(T.end(), ijk.end()))), // input - end chem_flag.begin(), // stencil chem_bgn[H], // output detail::chem_electroneutral(), // op diff --git a/src/impl/particles_impl_chem_henry.ipp b/src/impl/particles_impl_chem_henry.ipp index df92f5a02..d9c109606 100644 --- a/src/impl/particles_impl_chem_henry.ipp +++ b/src/impl/particles_impl_chem_henry.ipp @@ -335,7 +335,7 @@ namespace libcloudphxx V.begin(), V.end(), // input - 1st arg thrust::make_zip_iterator(thrust::make_tuple( // input - 2nd arg thrust::make_permutation_iterator(p.begin(), ijk.begin()), - thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()), + thrust::make_permutation_iterator(T.begin(), ijk.begin()), thrust::make_permutation_iterator(ambient_chem[(chem_species_t)i].begin(), ijk.begin()), chem_bgn[i], rw2.begin(), diff --git a/src/impl/particles_impl_chem_react.ipp b/src/impl/particles_impl_chem_react.ipp index f79f7fc76..b95e750fe 100644 --- a/src/impl/particles_impl_chem_react.ipp +++ b/src/impl/particles_impl_chem_react.ipp @@ -278,7 +278,7 @@ namespace libcloudphxx detail::chem_rhs( dt, V, - thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()), + thrust::make_permutation_iterator(T.begin(), ijk.begin()), chem_bgn[H], chem_flag ), // TODO: make it an impl member field diff --git a/src/impl/particles_impl_chem_strength.ipp b/src/impl/particles_impl_chem_strength.ipp index e485c90c0..43338ea9f 100644 --- a/src/impl/particles_impl_chem_strength.ipp +++ b/src/impl/particles_impl_chem_strength.ipp @@ -97,10 +97,10 @@ namespace libcloudphxx thrust::transform( zip_it_t(thrust::make_tuple( chem_bgn[SO2], chem_bgn[CO2], chem_bgn[HNO3], chem_bgn[NH3], chem_bgn[S_VI], chem_bgn[H], - V.begin(), thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()), rw2.begin())), // input - begin + V.begin(), thrust::make_permutation_iterator(T.begin(), ijk.begin()), rw2.begin())), // input - begin zip_it_t(thrust::make_tuple( chem_end[SO2], chem_end[CO2], chem_end[HNO3], chem_end[NH3], chem_end[S_VI], chem_end[H], - V.end(), thrust::make_permutation_iterator(T_ref.end(), ijk.end_ref()), rw2.end())), // input - end + V.end(), thrust::make_permutation_iterator(T.end(), ijk.end()), rw2.end())), // input - end chem_flag.begin(), // output detail::set_chem_flag() // op ); diff --git a/src/impl/particles_impl_hskpng_Tpr.ipp b/src/impl/particles_impl_hskpng_Tpr.ipp index 6b8a49af8..dd882fad6 100644 --- a/src/impl/particles_impl_hskpng_Tpr.ipp +++ b/src/impl/particles_impl_hskpng_Tpr.ipp @@ -177,9 +177,9 @@ namespace libcloudphxx { // T = common::theta_dry::T(th, rhod); thrust::transform( - th_ref.begin(), th_ref.end(), // input - first arg - rhod_ref.begin(), // input - second arg - T_ref.begin(), // output + th.begin_ref(), th.end_ref(), // input - first arg + rhod.begin_ref(), // input - second arg + T.begin_ref(), // output detail::common__theta_dry__T_rhod() ); } @@ -187,9 +187,9 @@ namespace libcloudphxx { // T = th * exner(p_tot) thrust::transform( - th_ref.begin(), th_ref.end(), // input - first arg - thrust::make_zip_iterator(thrust::make_tuple(rv_ref.begin(), p_ref.begin())), // input - second and third args - T_ref.begin(), // output + th.begin_ref(), th.end_ref(), // input - first arg + thrust::make_zip_iterator(thrust::make_tuple(rv.begin_ref(), p.begin_ref())), // input - second and third args + T.begin_ref(), // output detail::common__theta_dry__T_p() ); } @@ -207,17 +207,17 @@ namespace libcloudphxx { // p = common::theta_dry::p(rhod, r, T); thrust::transform( - zip_it_t(thrust::make_tuple(rhod_ref.begin(), rv_ref.begin(), T_ref.begin())), // input - begin - zip_it_t(thrust::make_tuple(rhod_ref.end(), rv_ref.end(), T_ref.end() )), // input - end - p_ref.begin(), // output + zip_it_t(thrust::make_tuple(rhod.begin_ref(), rv.begin_ref(), T.begin_ref())), // input - begin + zip_it_t(thrust::make_tuple(rhod.end_ref(), rv.end_ref(), T.end_ref() )), // input - end + p.begin_ref(), // output detail::common__theta_dry__p() ); } thrust::transform( - zip_it_t(thrust::make_tuple(p_ref.begin(), rv_ref.begin(), T_ref.begin())), // input - begin - zip_it_t(thrust::make_tuple(p_ref.end(), rv_ref.end(), T_ref.end() )), // input - end - RH_ref.begin(), // output + zip_it_t(thrust::make_tuple(p.begin_ref(), rv.begin_ref(), T.begin_ref())), // input - begin + zip_it_t(thrust::make_tuple(p.end_ref(), rv.end_ref(), T.end_ref() )), // input - end + RH.begin_ref(), // output detail::RH(opts_init.RH_formula) ); } @@ -226,8 +226,8 @@ namespace libcloudphxx { // on refined grid thrust::transform( - T_ref.begin(), T_ref.end(), // 1st arg - eta_ref.begin(), // output + T.begin_ref(), T.end_ref(), // 1st arg + eta.begin_ref(), // output detail::common__vterm__visc() ); // on normal grid: how if we dont know T, because we dont know th? diff --git a/src/impl/particles_impl_hskpng_ijk.ipp b/src/impl/particles_impl_hskpng_ijk.ipp index 081d4ec2d..113e79676 100644 --- a/src/impl/particles_impl_hskpng_ijk.ipp +++ b/src/impl/particles_impl_hskpng_ijk.ipp @@ -59,7 +59,7 @@ namespace libcloudphxx arg::_2 * opts_init.nz ); thrust::transform( - k.begin()+begin_shift, k.end() + k.begin()+begin_shift, k.end(), ijk_begin+begin_shift, ijk_begin+begin_shift, // in-place! arg::_1 + arg::_2 diff --git a/src/impl/particles_impl_hskpng_remove.ipp b/src/impl/particles_impl_hskpng_remove.ipp index 6ef943dc5..1d4e343e0 100644 --- a/src/impl/particles_impl_hskpng_remove.ipp +++ b/src/impl/particles_impl_hskpng_remove.ipp @@ -25,7 +25,7 @@ namespace libcloudphxx if (opts_init.ny != 0) n_t_vctrs.insert(&j); if (opts_init.nz != 0) n_t_vctrs.insert(&k); - std::set*> n_t_vctrs_ref; + std::set*> n_t_vctrs_ref; n_t_vctrs_ref.insert(&ijk); namespace arg = thrust::placeholders; diff --git a/src/impl/particles_impl_hskpng_vterm.ipp b/src/impl/particles_impl_hskpng_vterm.ipp index cf9c1ff22..024581daa 100644 --- a/src/impl/particles_impl_hskpng_vterm.ipp +++ b/src/impl/particles_impl_hskpng_vterm.ipp @@ -172,10 +172,10 @@ namespace libcloudphxx thrust::transform_if( rw2.begin(), rw2.end(), // input - 1st arg zip_it_t(thrust::make_tuple( - thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()), - thrust::make_permutation_iterator(p.begin(), ijk.begin()), - thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), - thrust::make_permutation_iterator(eta.begin(), ijk.begin()) + thrust::make_permutation_iterator(T.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(p.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(rhod.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(eta.begin_ref(), ijk.begin_ref()) )), // input - 2nd arg vt.begin(), // condition argument vt.begin(), // output @@ -219,10 +219,10 @@ namespace libcloudphxx thrust::transform( rw2.begin(), rw2.end(), // input - 1st arg zip_it_t(thrust::make_tuple( - thrust::make_permutation_iterator(T_ref.begin(), ijk.begin_ref()), - thrust::make_permutation_iterator(p.begin(), ijk.begin()), - thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), - thrust::make_permutation_iterator(eta.begin(), ijk.begin()) + thrust::make_permutation_iterator(T.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(p.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(rhod.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(eta.begin_ref(), ijk.begin_ref()) )), // input - 2nd arg vt.begin(), // output detail::common__vterm__vt(opts_init.terminal_velocity) diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index f5e1f8f79..5f267a15e 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -13,41 +13,26 @@ namespace libcloudphxx void particles_t::impl::init_hskpng_ncell() { // memory allocation, p already resized in init_sync() - eta.resize(n_cell.get()); + eta.resize(n_cell); - count_ijk.resize(n_cell.get()); - count_num.resize(n_cell.get()); - count_mom.resize(n_cell.get()); + count_ijk.resize(n_cell); + count_num.resize(n_cell); + count_mom.resize(n_cell); count_n.get() = 0; count_n.get_ref() = 0; // initialising device temporary arrays tmp_device_real_cell.resize(n_cell.get()); - tmp_device_real_cell1.resize(n_cell.get()); - tmp_device_real_cell2.resize(n_cell.get()); + tmp_device_real_cell1.resize(n_cell); + tmp_device_real_cell2.resize(n_cell); tmp_device_size_cell.resize(n_cell.get()); tmp_host_size_cell.resize(n_cell.get()); tmp_host_real_cell.resize(n_cell.get()); if(allow_sstp_cond && !opts_init.exact_sstp_cond) { - sstp_tmp_rv.resize(n_cell.get()); - sstp_tmp_th.resize(n_cell.get()); - sstp_tmp_rh.resize(n_cell.get()); - } - - // arrays on the refined grid - if(opts_init.n_ref > 1) - { - T_ref.resize(n_cell_ref); - RH_ref.resize(n_cell_ref); - eta_ref.resize(n_cell_ref); - - tmp_device_real_cell_ref.resize(n_cell_ref); - tmp_device_real_cell_ref1.resize(n_cell_ref); - - count_ijk_ref.resize(n_cell_ref); - count_num_ref.resize(n_cell_ref); - count_mom_ref.resize(n_cell_ref); + sstp_tmp_rv.resize(0, n_cell.get_ref()); + sstp_tmp_th.resize(0, n_cell.get_ref()); + sstp_tmp_rh.resize(0, n_cell.get_ref()); } } }; diff --git a/src/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index f53bd54fa..65edd3924 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/src/impl/particles_impl_init_sync.ipp @@ -13,11 +13,10 @@ namespace libcloudphxx void particles_t::impl::init_sync() { // memory allocation for scalar fields - rhod.resize(n_cell.get()); - rhod_ref.resize(n_cell_ref); - p_ref.resize(n_cell_ref); - th_ref.resize(n_cell_ref); - rv_ref.resize(n_cell_ref); + rhod.resize(n_cell); + p.resize(n_cell); + th.resize(n_cell); + rv.resize(n_cell); if(opts_init.chem_switch) for (int i = 0; i < chem_gas_n; ++i) ambient_chem[(chem_species_t)i].resize(n_cell.get()); diff --git a/src/impl/particles_impl_init_wet.ipp b/src/impl/particles_impl_init_wet.ipp index 046c606aa..a0b2e8a61 100644 --- a/src/impl/particles_impl_init_wet.ipp +++ b/src/impl/particles_impl_init_wet.ipp @@ -63,8 +63,8 @@ namespace libcloudphxx zip_it_t zip_it(thrust::make_tuple( rd3.begin() + n_part_old, kpa.begin() + n_part_old, - pi_t(RH_ref.begin(), ijk.begin_ref() + n_part_old), - pi_t(T_ref.begin(), ijk.begin_ref() + n_part_old) + pi_t(RH.begin_ref(), ijk.begin_ref() + n_part_old), + pi_t(T.begin_ref(), ijk.begin_ref() + n_part_old) )); thrust::transform( diff --git a/src/impl/particles_impl_update_incloud_time.ipp b/src/impl/particles_impl_update_incloud_time.ipp index 9db16eb35..b85937a46 100644 --- a/src/impl/particles_impl_update_incloud_time.ipp +++ b/src/impl/particles_impl_update_incloud_time.ipp @@ -44,7 +44,7 @@ namespace libcloudphxx thrust::make_zip_iterator(make_tuple( kpa.begin(), thrust::make_permutation_iterator( - T_ref.begin(), + T.begin_ref(), ijk.begin_ref() ) )), // input - 2nd arg diff --git a/src/impl/particles_impl_update_th_rv.ipp b/src/impl/particles_impl_update_th_rv.ipp index 7c65227ce..375cbadf1 100644 --- a/src/impl/particles_impl_update_th_rv.ipp +++ b/src/impl/particles_impl_update_th_rv.ipp @@ -77,7 +77,7 @@ namespace libcloudphxx thrust::make_transform_iterator( zip_it_t(thrust::make_tuple( drv.begin(), // - T_ref.begin(), // dth = drv * d_th_d_rv(T, th) + T.begin(), // dth = drv * d_th_d_rv(T, th) th.begin() // )), detail::dth() diff --git a/src/particles_init.ipp b/src/particles_init.ipp index 18ab2e830..fc3d2e0b1 100644 --- a/src/particles_init.ipp +++ b/src/particles_init.ipp @@ -39,11 +39,12 @@ namespace libcloudphxx // initialising Eulerian-Lagrangian coupling pimpl->init_sync(); // also, init of ambient_chem vectors - pimpl->init_e2l(th, &pimpl->th, true); - pimpl->init_e2l(rv, &pimpl->rv, true); - pimpl->init_e2l(rhod, &pimpl->rhod); + pimpl->init_e2l(th, pimpl->th.ptr_ref(), true); + pimpl->init_e2l(rv, pimpl->rv.ptr_ref(), true); + //pimpl->init_e2l(rhod, &pimpl->rhod); + pimpl->init_e2l(rhod, pimpl->rhod.ptr_ref(), true); if(pimpl->const_p) - pimpl->init_e2l(p, &pimpl->p); + pimpl->init_e2l(p, pimpl->p.ptr_ref(), true); #if !defined(__NVCC__) using std::max; @@ -57,10 +58,10 @@ namespace libcloudphxx pimpl->init_e2l(ambient_chem.at((chem_species_t)i), &pimpl->ambient_chem[(chem_species_t)i]); // feeding in Eulerian fields - pimpl->sync(th, pimpl->th); - pimpl->sync(rv, pimpl->rv); - pimpl->sync(rhod, pimpl->rhod); - pimpl->sync(p, pimpl->p); + pimpl->sync(th, pimpl->th.get_ref()); + pimpl->sync(rv, pimpl->rv.get_ref()); + pimpl->sync(rhod, pimpl->rhod.get_ref()); + pimpl->sync(p, pimpl->p.get_ref()); pimpl->sync(courant_x, pimpl->courant_x); pimpl->sync(courant_y, pimpl->courant_y); From da1bbedcd3825211106e32c80026354c2d7d7c3f Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 20 Jul 2022 20:18:45 +0200 Subject: [PATCH 21/38] _ref WIP; TODO: compilation errors --- src/particles_step.ipp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 608b85748..5419bdced 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -86,8 +86,8 @@ namespace libcloudphxx if (!courant_z.is_null()) pimpl->init_e2l(courant_z, &pimpl->courant_z, false, 0, 0, 1, pimpl->n_x_bfr * std::max(1, pimpl->opts_init.ny) - pimpl->halo_z); } - if (pimpl->l2e[&pimpl->diss_rate].size() == 0) - if (!diss_rate.is_null()) pimpl->init_e2l(diss_rate, &pimpl->diss_rate); + if (pimpl->l2e[pimpl->diss_rate.ptr()].size() == 0) + if (!diss_rate.is_null()) pimpl->init_e2l(diss_rate, pimpl->diss_rate.ptr()); // did rhod change pimpl->var_rho = !rhod.is_null(); From 12e11eec907d2038fa9087d58d10772809548ee7 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 21 Sep 2022 11:44:19 +0200 Subject: [PATCH 22/38] ref grid: init in ctor --- src/detail/grid_refinement_helpers.hpp | 77 ++++++++++++++----- src/impl/particles_impl.ipp | 45 +++++++++-- src/impl/particles_impl_init_hskpng_ncell.ipp | 18 +++-- src/impl/particles_impl_init_sync.ipp | 8 +- src/impl/particles_impl_sstp.ipp | 2 +- src/particles_diag.ipp | 2 +- src/particles_step.ipp | 4 +- 7 files changed, 113 insertions(+), 43 deletions(-) diff --git a/src/detail/grid_refinement_helpers.hpp b/src/detail/grid_refinement_helpers.hpp index 7d03797a3..e71dffc1d 100644 --- a/src/detail/grid_refinement_helpers.hpp +++ b/src/detail/grid_refinement_helpers.hpp @@ -17,6 +17,8 @@ namespace libcloudphxx ref_common(const bool ref_flag): ref_flag(ref_flag) {} + + ref_common() = delete; }; @@ -34,6 +36,7 @@ namespace libcloudphxx val_ref(val_ref), parent_t(true) {} + ref_val() = delete; T& get() { @@ -65,27 +68,11 @@ namespace libcloudphxx private: using parent_t = ref_common; using vec_t = thrust_device::vector; - vec_t arr, arr_ref; // actual data, arr_ref initialized only if refinement is actually done protected: using parent_t::parent_t; + vec_t arr, arr_ref; // actual data, arr_ref initialized only if refinement is actually done - auto begin() - { - return arr.begin(); - } - auto begin_ref() - { - return this->ref_flag ? arr_ref.begin() : arr.begin(); - } - auto end() - { - return arr.end(); - } - auto end_ref() - { - return this->ref_flag ? arr_ref.end() : arr.end(); - } void resize(thrust_size_t size, thrust_size_t size_ref) { arr.resize(size); @@ -109,6 +96,26 @@ namespace libcloudphxx public: + auto begin() + { + return arr.begin(); + } + + auto begin_ref() + { + return this->ref_flag ? arr_ref.begin() : arr.begin(); + } + + auto end() + { + return arr.end(); + } + + auto end_ref() + { + return this->ref_flag ? arr_ref.end() : arr.end(); + } + auto ptr() { return &arr; @@ -130,7 +137,8 @@ namespace libcloudphxx } }; - // for arrays of the size of the grid + // for arrays of the size of the grid that need values + // in the normal grid and in the refined grid (or in only one of these) template class ref_grid : public ref_arr_common { @@ -147,11 +155,36 @@ namespace libcloudphxx ref_grid(ref_val n): ref_grid(n.get(), n.get_ref()) {} + + ref_grid() = delete; + + // grid size cant be 0, ensure that correct arrays are pointed to + auto begin() + { + assert(this->arr.size() > 0); + return parent_t::begin(); + } + auto begin_ref() + { + assert((this->ref_flag && this->arr_ref.size() > 0) || (!this->ref_flag && this->arr.size() > 0)); + return parent_t::begin_ref(); + } + auto end() + { + assert(this->arr.size() > 0); + return parent_t::end(); + } + auto end_ref() + { + assert((this->ref_flag && this->arr_ref.size() > 0) || (!this->ref_flag && this->arr.size() > 0)); + return parent_t::end_ref(); + } }; - // for arrays of the size of the number of particles + // for arrays of the size of the number of particles that need + // two independent values: for the normal and for the refined grid (e.g. ijk) template - class ref_part : ref_arr_common + class ref_part : public ref_arr_common { using parent_t = ref_arr_common; @@ -159,7 +192,9 @@ namespace libcloudphxx // ctor ref_part(const int &n_ref): parent_t(n_ref > 1) - {} + {assert(n_ref > 0);} + + ref_part() = delete; void resize(int n) { diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index ac1e73c59..9efe701c3 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -85,7 +85,10 @@ namespace libcloudphxx wp, // turbulent perturbation of velocity ssp, // turbulent perturbation of supersaturation dot_ssp, // time derivative of the turbulent perturbation of supersaturation - sstp_tmp_p, // ditto for pressure + sstp_pp_tmp_rv, // helper for per-particle cond substepping + sstp_pp_tmp_rh, // ditto for density + sstp_pp_tmp_th, // ditto for temperature + sstp_pp_tmp_p, // ditto for pressure incloud_time; // time this SD has been within a cloud // dry radii distribution characteristics @@ -135,9 +138,9 @@ namespace libcloudphxx rhod, // dry air density th, // potential temperature (dry) rv, // water vapour mixing ratio - sstp_tmp_rv, // either rv_old or advection-caused change in water vapour mixing ratio - sstp_tmp_th, // ditto for theta - sstp_tmp_rh; // ditto for rho + sstp_pc_tmp_rv, // either rv_old or advection-caused change in water vapour mixing ratio, per-cell substepping + sstp_pc_tmp_th, // ditto for theta + sstp_pc_tmp_rh; // ditto for rho thrust_device::vector sstp_tmp_chem_0, // trace gases @@ -154,13 +157,15 @@ namespace libcloudphxx // map of the accumulated volume/volume/mass of water/dry/chem that fell out of the domain std::map output_puddle; + + thrust_device::vector + diss_rate; // turbulent kinetic energy dissipation rate ref_grid T, // temperature [K] p, // pressure [Pa] RH, // relative humisity - eta, // dynamic viscosity - diss_rate; // turbulent kinetic energy dissipation rate + eta; // dynamic viscosity thrust_device::vector w_LS; // large-scale subsidence velocity profile thrust_device::vector SGS_mix_len; // SGS mixing length profile @@ -183,6 +188,14 @@ namespace libcloudphxx // is it allowed to do substepping, if not, some memory can be saved bool allow_sstp_cond, allow_sstp_chem; + + // aliases to substepping arrays, they point to per-cell or per-particle arrays + // depending on runtime options + thrust_device::vector + &sstp_tmp_rv, + &sstp_tmp_rh, + &sstp_tmp_th, + &sstp_tmp_p; // timestep counter n_t src_stp_ctr, rlx_stp_ctr; @@ -317,6 +330,19 @@ namespace libcloudphxx nz_ref(n_dims >= 2 ? (opts_init.nz - 1) * opts_init.n_ref + 1 : 0), n_cell(m1(_opts_init.nx) * m1(_opts_init.ny) * m1(_opts_init.nz), m1(nx_ref) * m1(ny_ref) * m1(nz_ref)), // NOTE: needs to be equal to n_cell for n_ref == 1 + rhod(n_cell), + p(0, n_cell.get_ref()), + th(0, n_cell.get_ref()), + rv(0, n_cell.get_ref()), + count_ijk(n_cell), + count_num(n_cell), + count_mom(n_cell), + count_n(0,0), + tmp_device_real_cell1(n_cell), + tmp_device_real_cell2(n_cell), + T (0, n_cell.get_ref()), + RH(0, n_cell.get_ref()), + eta(n_cell), zero(0), n_part(0), sorted(false), @@ -352,6 +378,13 @@ namespace libcloudphxx adve_scheme(_opts_init.adve_scheme), allow_sstp_cond(_opts_init.sstp_cond > 1 || _opts_init.variable_dt_switch), allow_sstp_chem(_opts_init.sstp_chem > 1 || _opts_init.variable_dt_switch), + sstp_pc_tmp_rv(0, n_cell.get_ref()), + sstp_pc_tmp_th(0, n_cell.get_ref()), + sstp_pc_tmp_rh(0, n_cell.get_ref()), + sstp_tmp_rv(_opts_init.exact_sstp_cond ? &sstp_pp_tmp_rv : sstp_pc_tmp_rv.ptr()), + sstp_tmp_th(_opts_init.exact_sstp_cond ? &sstp_pp_tmp_th : sstp_pc_tmp_th.ptr()), + sstp_tmp_rh(_opts_init.exact_sstp_cond ? &sstp_pp_tmp_rh : sstp_pc_tmp_rh.ptr()), + sstp_tmp_p(sstp_pp_tmp_p), // not used in per-cell substepping pure_const_multi (((_opts_init.sd_conc) == 0) && (_opts_init.sd_const_multi > 0 || _opts_init.dry_sizes.size() > 0)), // coal prob can be greater than one only in sd_conc simulations ijk(opts_init.n_ref) { diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index 5f267a15e..353a4ff95 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -13,18 +13,20 @@ namespace libcloudphxx void particles_t::impl::init_hskpng_ncell() { // memory allocation, p already resized in init_sync() - eta.resize(n_cell); + // eta.resize(n_cell); + // T.resize(n_cell); + // RH.resize(n_cell); - count_ijk.resize(n_cell); - count_num.resize(n_cell); - count_mom.resize(n_cell); - count_n.get() = 0; - count_n.get_ref() = 0; +// count_ijk.resize(n_cell); +// count_num.resize(n_cell); +// count_mom.resize(n_cell); +// count_n.get() = 0; +// count_n.get_ref() = 0; // initialising device temporary arrays tmp_device_real_cell.resize(n_cell.get()); - tmp_device_real_cell1.resize(n_cell); - tmp_device_real_cell2.resize(n_cell); +// tmp_device_real_cell1.resize(n_cell); +// tmp_device_real_cell2.resize(n_cell); tmp_device_size_cell.resize(n_cell.get()); tmp_host_size_cell.resize(n_cell.get()); tmp_host_real_cell.resize(n_cell.get()); diff --git a/src/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index 65edd3924..954c94dd3 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/src/impl/particles_impl_init_sync.ipp @@ -13,10 +13,10 @@ namespace libcloudphxx void particles_t::impl::init_sync() { // memory allocation for scalar fields - rhod.resize(n_cell); - p.resize(n_cell); - th.resize(n_cell); - rv.resize(n_cell); +// rhod.resize(n_cell); +// p.resize(n_cell); +// th.resize(n_cell); +// rv.resize(n_cell); if(opts_init.chem_switch) for (int i = 0; i < chem_gas_n; ++i) ambient_chem[(chem_species_t)i].resize(n_cell.get()); diff --git a/src/impl/particles_impl_sstp.ipp b/src/impl/particles_impl_sstp.ipp index 28ce6191e..50b013240 100644 --- a/src/impl/particles_impl_sstp.ipp +++ b/src/impl/particles_impl_sstp.ipp @@ -108,7 +108,7 @@ namespace libcloudphxx void particles_t::impl::sstp_step_exact( const int &step ) - { + { if (sstp_cond == 1) return; namespace arg = thrust::placeholders; diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index 7905506bc..5ebdc707c 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -120,7 +120,7 @@ namespace libcloudphxx // p defined in all cells pimpl->count_n.get() = pimpl->n_cell.get(); - thrust::sequence(pimpl->count_ijk.begin(), pimpl->count_ijk.end()); + thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } // records temperature diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 5419bdced..608b85748 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -86,8 +86,8 @@ namespace libcloudphxx if (!courant_z.is_null()) pimpl->init_e2l(courant_z, &pimpl->courant_z, false, 0, 0, 1, pimpl->n_x_bfr * std::max(1, pimpl->opts_init.ny) - pimpl->halo_z); } - if (pimpl->l2e[pimpl->diss_rate.ptr()].size() == 0) - if (!diss_rate.is_null()) pimpl->init_e2l(diss_rate, pimpl->diss_rate.ptr()); + if (pimpl->l2e[&pimpl->diss_rate].size() == 0) + if (!diss_rate.is_null()) pimpl->init_e2l(diss_rate, &pimpl->diss_rate); // did rhod change pimpl->var_rho = !rhod.is_null(); From f91c9ea63611843d740934d5368e428fd7af71eb Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 21 Sep 2022 15:53:38 +0200 Subject: [PATCH 23/38] wip on ref, compilation errors' --- include/libcloudph++/lgrngn/opts_init.hpp | 8 ++-- src/detail/grid_refinement_helpers.hpp | 11 +++-- src/impl/particles_impl.ipp | 44 ++++++++++--------- src/impl/particles_impl_cond.ipp | 2 +- src/impl/particles_impl_fill_outbuf.ipp | 17 ++++--- src/impl/particles_impl_hskpng_Tpr.ipp | 21 ++++++--- src/impl/particles_impl_init_hskpng_ncell.ipp | 2 + src/impl/particles_impl_init_ijk.ipp | 2 +- src/impl/particles_impl_init_sync.ipp | 2 +- src/impl/particles_impl_sstp.ipp | 4 +- src/impl/particles_impl_update_th_rv.ipp | 13 ++++++ src/particles_diag.ipp | 8 ++-- src/particles_init.ipp | 6 +-- src/particles_step.ipp | 8 ++-- 14 files changed, 91 insertions(+), 57 deletions(-) diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index 121c3648b..65e1a7991 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -49,9 +49,11 @@ namespace libcloudphxx unsigned int n_ref; // number of grid refinements per 'normal' cell; refined cell size is (dx, dy, dz) / n_ref; // condensation (and sedimentation?) is resolved on refined cells, rest (coalescence, advection, etc.) on normal cells // input arrays on the refined grid: pressure p (if supplied), rv, tht - // input arrays that require two versions, on the normal and on the refined grid: density (rhod, rhod_ref) - // internal arrays on the refined grid: RH, T - // internal arrays that require two versions, on the normal and on the refined grid: dynamic viscosity eta (eta, eta_ref), cell volume dv (not sure...) + // internal arrays on the refined grid: RH, T, eta + // TODO: + // dynamic viscosity eta on both grids? + // what about cell volume dv? + // what about density rhod? // no. of substeps int sstp_cond, sstp_coal; diff --git a/src/detail/grid_refinement_helpers.hpp b/src/detail/grid_refinement_helpers.hpp index e71dffc1d..09083ea3a 100644 --- a/src/detail/grid_refinement_helpers.hpp +++ b/src/detail/grid_refinement_helpers.hpp @@ -9,10 +9,9 @@ namespace libcloudphxx template class ref_common { - private: + protected: const bool ref_flag; // true if refinement is actually done - protected: // ctor ref_common(const bool ref_flag): ref_flag(ref_flag) @@ -161,22 +160,22 @@ namespace libcloudphxx // grid size cant be 0, ensure that correct arrays are pointed to auto begin() { - assert(this->arr.size() > 0); + assert(this->arr.size() > 0); return parent_t::begin(); } auto begin_ref() { - assert((this->ref_flag && this->arr_ref.size() > 0) || (!this->ref_flag && this->arr.size() > 0)); + assert((this->ref_flag && this->arr_ref.size() > 0) || (!this->ref_flag && this->arr.size() > 0)); return parent_t::begin_ref(); } auto end() { - assert(this->arr.size() > 0); + assert(this->arr.size() > 0); return parent_t::end(); } auto end_ref() { - assert((this->ref_flag && this->arr_ref.size() > 0) || (!this->ref_flag && this->arr.size() > 0)); + assert((this->ref_flag && this->arr_ref.size() > 0) || (!this->ref_flag && this->arr.size() > 0)); return parent_t::end_ref(); } }; diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 9efe701c3..0d6f4fe7b 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -109,6 +109,7 @@ namespace libcloudphxx i, j, k, // Eulerian grid cell indices (always zero for 0D) sorted_id, sorted_ijk; ref_part ijk; // ijk in the normal and in the refined grid + ref_grid ijk_ref2ijk; // maps refined cell index to the index of the normal cell containing this refined cell // helpers that reference refined arrays if refinement is done and non-refined otherwise // used to conserve memory if n_ref==1 (no refinement) @@ -135,7 +136,6 @@ namespace libcloudphxx // Eulerian-Lagrangian interface vars ref_grid - rhod, // dry air density th, // potential temperature (dry) rv, // water vapour mixing ratio sstp_pc_tmp_rv, // either rv_old or advection-caused change in water vapour mixing ratio, per-cell substepping @@ -159,6 +159,7 @@ namespace libcloudphxx std::map output_puddle; thrust_device::vector + rhod, // dry air density diss_rate; // turbulent kinetic energy dissipation rate ref_grid @@ -330,19 +331,19 @@ namespace libcloudphxx nz_ref(n_dims >= 2 ? (opts_init.nz - 1) * opts_init.n_ref + 1 : 0), n_cell(m1(_opts_init.nx) * m1(_opts_init.ny) * m1(_opts_init.nz), m1(nx_ref) * m1(ny_ref) * m1(nz_ref)), // NOTE: needs to be equal to n_cell for n_ref == 1 - rhod(n_cell), - p(0, n_cell.get_ref()), - th(0, n_cell.get_ref()), - rv(0, n_cell.get_ref()), - count_ijk(n_cell), - count_num(n_cell), - count_mom(n_cell), - count_n(0,0), - tmp_device_real_cell1(n_cell), - tmp_device_real_cell2(n_cell), - T (0, n_cell.get_ref()), - RH(0, n_cell.get_ref()), - eta(n_cell), + p(0, n_cell.get_ref()), + th(0, n_cell.get_ref()), + rv(0, n_cell.get_ref()), + ijk_ref2ijk(0, n_cell.get_ref()), + count_ijk(n_cell), + count_num(n_cell), + count_mom(n_cell), + count_n(0,0), + tmp_device_real_cell1(n_cell), + tmp_device_real_cell2(n_cell), + T (0, n_cell.get_ref()), + RH(0, n_cell.get_ref()), + eta(0, n_cell.get_ref()), zero(0), n_part(0), sorted(false), @@ -378,13 +379,13 @@ namespace libcloudphxx adve_scheme(_opts_init.adve_scheme), allow_sstp_cond(_opts_init.sstp_cond > 1 || _opts_init.variable_dt_switch), allow_sstp_chem(_opts_init.sstp_chem > 1 || _opts_init.variable_dt_switch), - sstp_pc_tmp_rv(0, n_cell.get_ref()), - sstp_pc_tmp_th(0, n_cell.get_ref()), - sstp_pc_tmp_rh(0, n_cell.get_ref()), - sstp_tmp_rv(_opts_init.exact_sstp_cond ? &sstp_pp_tmp_rv : sstp_pc_tmp_rv.ptr()), - sstp_tmp_th(_opts_init.exact_sstp_cond ? &sstp_pp_tmp_th : sstp_pc_tmp_th.ptr()), - sstp_tmp_rh(_opts_init.exact_sstp_cond ? &sstp_pp_tmp_rh : sstp_pc_tmp_rh.ptr()), - sstp_tmp_p(sstp_pp_tmp_p), // not used in per-cell substepping + sstp_pc_tmp_rv(0, n_cell.get_ref()), + sstp_pc_tmp_th(0, n_cell.get_ref()), + sstp_pc_tmp_rh(0, n_cell.get_ref()), + sstp_tmp_rv(_opts_init.exact_sstp_cond ? sstp_pp_tmp_rv : sstp_pc_tmp_rv.get_ref()), + sstp_tmp_th(_opts_init.exact_sstp_cond ? sstp_pp_tmp_th : sstp_pc_tmp_th.get_ref()), + sstp_tmp_rh(_opts_init.exact_sstp_cond ? sstp_pp_tmp_rh : sstp_pc_tmp_rh.get_ref()), + sstp_tmp_p(sstp_pp_tmp_p), // not used in per-cell substepping pure_const_multi (((_opts_init.sd_conc) == 0) && (_opts_init.sd_const_multi > 0 || _opts_init.dry_sizes.size() > 0)), // coal prob can be greater than one only in sd_conc simulations ijk(opts_init.n_ref) { @@ -611,6 +612,7 @@ namespace libcloudphxx 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_state(ref_grid &, thrust_device::vector &); void update_pstate(thrust_device::vector &, thrust_device::vector &); void update_incloud_time(const real_t &dt); diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index 59ad35308..d2ea5dfad 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -42,7 +42,7 @@ namespace libcloudphxx ); auto hlpr_zip_iter = thrust::make_zip_iterator(thrust::make_tuple( - thrust::make_permutation_iterator(rhod.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(rv.begin_ref(), ijk.begin_ref()), thrust::make_permutation_iterator(T.begin_ref(), ijk.begin_ref()), thrust::make_permutation_iterator(eta.begin_ref(), ijk.begin_ref()), diff --git a/src/impl/particles_impl_fill_outbuf.ipp b/src/impl/particles_impl_fill_outbuf.ipp index a38e2af08..4f2813ab8 100644 --- a/src/impl/particles_impl_fill_outbuf.ipp +++ b/src/impl/particles_impl_fill_outbuf.ipp @@ -28,13 +28,18 @@ namespace libcloudphxx #endif thrust::copy( - count_mom.begin(), // input - begin - count_mom.begin() + count_n.get(), // input - end - thrust::make_permutation_iterator( // output - tmp_host_real_cell.begin(), // data - pi.begin() // permutation - ) + count_mom.begin(), // input - begin + count_mom.begin() + count_n.get(), // input - end + thrust::make_permutation_iterator( // output + tmp_host_real_cell.begin(), // data + pi.begin() // permutation + ) ); } + + template + void particles_t::impl::fill_outbuf_ref() + { + } }; }; diff --git a/src/impl/particles_impl_hskpng_Tpr.ipp b/src/impl/particles_impl_hskpng_Tpr.ipp index dd882fad6..a80f11a69 100644 --- a/src/impl/particles_impl_hskpng_Tpr.ipp +++ b/src/impl/particles_impl_hskpng_Tpr.ipp @@ -178,7 +178,7 @@ namespace libcloudphxx // T = common::theta_dry::T(th, rhod); thrust::transform( th.begin_ref(), th.end_ref(), // input - first arg - rhod.begin_ref(), // input - second arg + thrust::make_permutation_iterator(rhod.begin(), ijk_ref2ijk.begin()), // input - second arg T.begin_ref(), // output detail::common__theta_dry__T_rhod() ); @@ -206,10 +206,15 @@ namespace libcloudphxx if(!const_p) { // p = common::theta_dry::p(rhod, r, T); + zip_it_t it(thrust::make_tuple( + thrust::make_permutation_iterator(rhod.begin(), ijk_ref2ijk.begin()), + rv.begin_ref(), + T.begin_ref() + )); thrust::transform( - zip_it_t(thrust::make_tuple(rhod.begin_ref(), rv.begin_ref(), T.begin_ref())), // input - begin - zip_it_t(thrust::make_tuple(rhod.end_ref(), rv.end_ref(), T.end_ref() )), // input - end - p.begin_ref(), // output + it, // input - begin + it + n_cell.get_ref(), // input - end + p.begin_ref(), // output detail::common__theta_dry__p() ); } @@ -231,9 +236,15 @@ namespace libcloudphxx detail::common__vterm__visc() ); // on normal grid: how if we dont know T, because we dont know th? - // TODO: use averaging from refined to normal? + // TODO: require input of normal th? + // or use averaging from refined to normal? + // cumbersome: + // 0. make copies of eta and ijk_ref2ijk + // 1. sort copies of eta by ijk_ref2ijk + // 2. reduce_by_key and add the result to normal eta } + // adjusting dv if using a parcel set-up (1kg of dry air) if (n_dims == 0) { diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index 353a4ff95..f81b1aea9 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -36,6 +36,8 @@ namespace libcloudphxx sstp_tmp_th.resize(0, n_cell.get_ref()); sstp_tmp_rh.resize(0, n_cell.get_ref()); } + + //TODO: init ijk_ref2ijk. it's size is already initialized } }; }; diff --git a/src/impl/particles_impl_init_ijk.ipp b/src/impl/particles_impl_init_ijk.ipp index 37bde74ad..c80c97ab7 100644 --- a/src/impl/particles_impl_init_ijk.ipp +++ b/src/impl/particles_impl_init_ijk.ipp @@ -46,7 +46,7 @@ namespace libcloudphxx thrust::make_zip_iterator(thrust::make_tuple( count_num.end(), ptr.end(), zero + n_cell.get() )), - detail::arbitrary_sequence(&(ijk.begin() + n_part_old)) + detail::arbitrary_sequence(&(*(ijk.begin() + n_part_old))) ); } }; diff --git a/src/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index 954c94dd3..616ac01e0 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/src/impl/particles_impl_init_sync.ipp @@ -13,7 +13,7 @@ namespace libcloudphxx void particles_t::impl::init_sync() { // memory allocation for scalar fields -// rhod.resize(n_cell); + rhod.resize(n_cell.get()); // p.resize(n_cell); // th.resize(n_cell); // rv.resize(n_cell); diff --git a/src/impl/particles_impl_sstp.ipp b/src/impl/particles_impl_sstp.ipp index 50b013240..7bbc340f3 100644 --- a/src/impl/particles_impl_sstp.ipp +++ b/src/impl/particles_impl_sstp.ipp @@ -45,8 +45,8 @@ namespace libcloudphxx const int n = 4; thrust_device::vector - *fr[n] = { &rv, &th, &rhod, &p }, - *to[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh, &sstp_tmp_p }; + *fr[n] = { rv.ptr_ref(), th.ptr_ref(), &rhod, p.ptr_ref() }, + *to[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh, &sstp_tmp_p }; for (int ix = 0; ix < ( (const_p) ? n : n-1); ++ix) // TODO: var_rho { thrust::copy( diff --git a/src/impl/particles_impl_update_th_rv.ipp b/src/impl/particles_impl_update_th_rv.ipp index 375cbadf1..3b5241a96 100644 --- a/src/impl/particles_impl_update_th_rv.ipp +++ b/src/impl/particles_impl_update_th_rv.ipp @@ -134,6 +134,7 @@ namespace libcloudphxx // update cell state based on particle-specific cell state // particles have to be sorted + /* template void particles_t::impl::update_state( thrust_device::vector &state, // cell state @@ -145,5 +146,17 @@ namespace libcloudphxx thrust::make_permutation_iterator(state.begin(), ijk.begin()) ); } + */ + template + void particles_t::impl::update_state( + ref_grid &state, // cell state + thrust_device::vector &pstate // particle-specific cell state (same for all particles in one cell) + ) + { + thrust::copy( + pstate.begin(), pstate.end(), + thrust::make_permutation_iterator(state.begin_ref(), ijk.begin_ref()) + ); + } }; }; diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index 5ebdc707c..35ad31216 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -113,9 +113,9 @@ namespace libcloudphxx pimpl->hskpng_Tpr(); thrust::copy( - pimpl->p.begin(), - pimpl->p.end(), - pimpl->count_mom.begin() + pimpl->p.begin_ref(), + pimpl->p.end_ref(), + pimpl->count_mom.begin_ref() ); // p defined in all cells @@ -285,7 +285,7 @@ namespace libcloudphxx thrust::make_zip_iterator(make_tuple( pimpl->kpa.begin(), thrust::make_permutation_iterator( - pimpl->T_ref.begin(), + pimpl->T.begin_ref(), pimpl->ijk.begin_ref() ) )), // input - 2nd arg diff --git a/src/particles_init.ipp b/src/particles_init.ipp index fc3d2e0b1..a9c8d5460 100644 --- a/src/particles_init.ipp +++ b/src/particles_init.ipp @@ -41,8 +41,8 @@ namespace libcloudphxx pimpl->init_sync(); // also, init of ambient_chem vectors pimpl->init_e2l(th, pimpl->th.ptr_ref(), true); pimpl->init_e2l(rv, pimpl->rv.ptr_ref(), true); - //pimpl->init_e2l(rhod, &pimpl->rhod); - pimpl->init_e2l(rhod, pimpl->rhod.ptr_ref(), true); + pimpl->init_e2l(rhod, &pimpl->rhod); + //pimpl->init_e2l(rhod, pimpl->rhod.ptr_ref(), true); if(pimpl->const_p) pimpl->init_e2l(p, pimpl->p.ptr_ref(), true); @@ -60,7 +60,7 @@ namespace libcloudphxx // feeding in Eulerian fields pimpl->sync(th, pimpl->th.get_ref()); pimpl->sync(rv, pimpl->rv.get_ref()); - pimpl->sync(rhod, pimpl->rhod.get_ref()); + pimpl->sync(rhod, pimpl->rhod); pimpl->sync(p, pimpl->p.get_ref()); pimpl->sync(courant_x, pimpl->courant_x); diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 608b85748..50de407e7 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -93,8 +93,8 @@ namespace libcloudphxx pimpl->var_rho = !rhod.is_null(); // syncing in Eulerian fields (if not null) - pimpl->sync(th, pimpl->th); - pimpl->sync(rv, pimpl->rv); + pimpl->sync(th, pimpl->th.get_ref()); + pimpl->sync(rv, pimpl->rv.get_ref()); pimpl->sync(diss_rate, pimpl->diss_rate); pimpl->sync(rhod, pimpl->rhod); pimpl->sync(courant_x, pimpl->courant_x); @@ -259,8 +259,8 @@ namespace libcloudphxx if(opts.cond) // || (opts.src && pimpl->src_stp_ctr % pimpl->opts_init.supstp_src == 0) || (opts.rlx && pimpl->rlx_stp_ctr % pimpl->opts_init.supstp_rlx == 0)) { // syncing out // TODO: this is not necesarry in off-line mode (see coupling with DALES) - pimpl->sync(pimpl->th, th); - pimpl->sync(pimpl->rv, rv); + pimpl->sync(pimpl->th.get_ref(), th); + pimpl->sync(pimpl->rv.get_ref(), rv); } if (opts.chem_dsl == true) From 139aa4be982bb69c0d8489b7daed7abe1dc978e4 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 22 Sep 2022 12:54:23 +0200 Subject: [PATCH 24/38] ref grid: init in ctor --- CMakeLists.txt | 2 +- src/detail/grid_refinement_helpers.hpp | 22 ++++++++------- src/impl/particles_impl.ipp | 12 ++++++--- src/impl/particles_impl_fill_outbuf.ipp | 27 +++++++++++++++++-- src/impl/particles_impl_hskpng_ijk.ipp | 2 +- src/impl/particles_impl_init_hskpng_ncell.ipp | 4 +-- src/impl/particles_impl_init_n.ipp | 2 +- 7 files changed, 50 insertions(+), 21 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 00d9fa58a..b1f7cefd1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -116,7 +116,7 @@ target_include_directories(cloudphxx_lgrngn # enabling c++14 target_compile_features(cloudphxx_lgrngn PUBLIC cxx_std_14) # config-specific flags -target_compile_options(cloudphxx_lgrngn PRIVATE $<$,$>: -g -Og -DTHRUST_DEBUG>) +target_compile_options(cloudphxx_lgrngn PRIVATE $<$,$>: -g -Og -DTHRUST_DEBUG -Wfatal-errors>) target_compile_options(cloudphxx_lgrngn PRIVATE $<$,$>: -Ofast -march=native -Winline -DNDEBUG >) target_compile_options(cloudphxx_lgrngn PRIVATE $<$,$>: -Ofast -march=native>) target_compile_options(cloudphxx_lgrngn PRIVATE $<$,$>: -Ofast>) diff --git a/src/detail/grid_refinement_helpers.hpp b/src/detail/grid_refinement_helpers.hpp index 09083ea3a..de0c71a2e 100644 --- a/src/detail/grid_refinement_helpers.hpp +++ b/src/detail/grid_refinement_helpers.hpp @@ -1,4 +1,5 @@ // helper classes to store arrays associated with grid refinement +#include #pragma once @@ -9,9 +10,10 @@ namespace libcloudphxx template class ref_common { - protected: + public: const bool ref_flag; // true if refinement is actually done + protected: // ctor ref_common(const bool ref_flag): ref_flag(ref_flag) @@ -61,12 +63,12 @@ namespace libcloudphxx */ }; - template - class ref_arr_common : ref_common + template + class ref_arr_common : public ref_common { private: using parent_t = ref_common; - using vec_t = thrust_device::vector; + using vec_t = typename std::conditional, thrust_device::vector>::type; protected: using parent_t::parent_t; @@ -138,10 +140,10 @@ namespace libcloudphxx // for arrays of the size of the grid that need values // in the normal grid and in the refined grid (or in only one of these) - template - class ref_grid : public ref_arr_common + template // value type and container type + class ref_grid : public ref_arr_common { - using parent_t = ref_arr_common; + using parent_t = ref_arr_common; public: // ctor @@ -182,10 +184,10 @@ namespace libcloudphxx // for arrays of the size of the number of particles that need // two independent values: for the normal and for the refined grid (e.g. ijk) - template - class ref_part : public ref_arr_common + template // value type and container type + class ref_part : public ref_arr_common { - using parent_t = ref_arr_common; + using parent_t = ref_arr_common; public: // ctor diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 0d6f4fe7b..6980174ec 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -232,10 +232,7 @@ namespace libcloudphxx // temporary data thrust::host_vector - tmp_host_real_grid, - tmp_host_real_cell; - thrust::host_vector - tmp_host_size_cell; + tmp_host_real_grid; thrust_device::vector tmp_device_real_part, tmp_device_real_part1, @@ -248,6 +245,10 @@ namespace libcloudphxx ref_grid tmp_device_real_cell1, tmp_device_real_cell2; + ref_grid + tmp_host_real_cell; + ref_grid + tmp_host_size_cell; thrust_device::vector tmp_device_n_part, &un; // uniform natural random numbers between 0 and max value of unsigned int @@ -341,6 +342,8 @@ namespace libcloudphxx count_n(0,0), tmp_device_real_cell1(n_cell), tmp_device_real_cell2(n_cell), + tmp_host_real_cell(n_cell), + tmp_host_size_cell(n_cell), T (0, n_cell.get_ref()), RH(0, n_cell.get_ref()), eta(0, n_cell.get_ref()), @@ -529,6 +532,7 @@ namespace libcloudphxx void init_vterm(); void fill_outbuf(); + void fill_outbuf_ref(); void mpi_exchange(); // rename hskpng_ -> step_? diff --git a/src/impl/particles_impl_fill_outbuf.ipp b/src/impl/particles_impl_fill_outbuf.ipp index 4f2813ab8..ab78cfd4f 100644 --- a/src/impl/particles_impl_fill_outbuf.ipp +++ b/src/impl/particles_impl_fill_outbuf.ipp @@ -22,9 +22,9 @@ namespace libcloudphxx #endif #if !defined(__NVCC__) - thrust_device::vector &pi(count_ijk); + auto &pi(count_ijk); #else - thrust::host_vector &pi(tmp_host_size_cell); + auto &pi(tmp_host_size_cell); #endif thrust::copy( @@ -40,6 +40,29 @@ namespace libcloudphxx template void particles_t::impl::fill_outbuf_ref() { + thrust::fill(tmp_host_real_cell.begin_ref(), tmp_host_real_cell.end_ref(), 0); + +#if defined(__NVCC__) + thrust::copy( + count_ijk.begin_ref(), count_ijk.end_ref(), // from + tmp_host_size_cell.begin_ref() + ); +#endif + +#if !defined(__NVCC__) + auto &pi(count_ijk); +#else + auto &pi(tmp_host_size_cell); +#endif + + thrust::copy( + count_mom.begin_ref(), // input - begin + count_mom.begin_ref() + count_n.get_ref(), // input - end + thrust::make_permutation_iterator( // output + tmp_host_real_cell.begin_ref(), // data + pi.begin_ref() // permutation + ) + ); } }; }; diff --git a/src/impl/particles_impl_hskpng_ijk.ipp b/src/impl/particles_impl_hskpng_ijk.ipp index 113e79676..a598c1125 100644 --- a/src/impl/particles_impl_hskpng_ijk.ipp +++ b/src/impl/particles_impl_hskpng_ijk.ipp @@ -32,7 +32,7 @@ namespace libcloudphxx template void particles_t::impl::ravel_ijk(const thrust_size_t begin_shift, const bool refined) { - auto &ijk_begin = refined ? ijk.begin_ref() : ijk.begin(); + auto ijk_begin = refined ? ijk.begin_ref() : ijk.begin(); switch (n_dims) { case 0: diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index f81b1aea9..31312507e 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -28,8 +28,8 @@ namespace libcloudphxx // tmp_device_real_cell1.resize(n_cell); // tmp_device_real_cell2.resize(n_cell); tmp_device_size_cell.resize(n_cell.get()); - tmp_host_size_cell.resize(n_cell.get()); - tmp_host_real_cell.resize(n_cell.get()); +// tmp_host_size_cell.resize(n_cell.get()); +// tmp_host_real_cell.resize(n_cell.get()); if(allow_sstp_cond && !opts_init.exact_sstp_cond) { sstp_tmp_rv.resize(0, n_cell.get_ref()); diff --git a/src/impl/particles_impl_init_n.ipp b/src/impl/particles_impl_init_n.ipp index 50c4880a7..74ae9e91f 100644 --- a/src/impl/particles_impl_init_n.ipp +++ b/src/impl/particles_impl_init_n.ipp @@ -52,7 +52,7 @@ namespace libcloudphxx // temporary space on the host thrust::host_vector tmp_real(n_part_to_init); thrust::host_vector tmp_ijk(n_part_to_init); - thrust::host_vector &tmp_rhod(tmp_host_real_cell); + thrust::host_vector &tmp_rhod(tmp_host_real_cell.get()); thrust::copy( rhod.begin(), rhod.end(), // from From a28b5c93ccd3c35c8511b42c2ea25c7eb97eb764 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 22 Sep 2022 13:14:02 +0200 Subject: [PATCH 25/38] ref grid: sstp tuff --- src/impl/particles_impl_chem_henry.ipp | 2 +- src/impl/particles_impl_hskpng_vterm.ipp | 4 ++-- src/impl/particles_impl_sstp.ipp | 27 ++++++++++++++---------- 3 files changed, 19 insertions(+), 14 deletions(-) diff --git a/src/impl/particles_impl_chem_henry.ipp b/src/impl/particles_impl_chem_henry.ipp index d9c109606..297a634d8 100644 --- a/src/impl/particles_impl_chem_henry.ipp +++ b/src/impl/particles_impl_chem_henry.ipp @@ -308,7 +308,7 @@ namespace libcloudphxx // temporarily needed to store old mass per cell thrust_device::vector &mass_old(tmp_device_real_cell); - thrust_device::vector &mass_new(tmp_device_real_cell1); + thrust_device::vector &mass_new(tmp_device_real_cell1.get()); for (int i = 0; i < chem_gas_n; ++i) { diff --git a/src/impl/particles_impl_hskpng_vterm.ipp b/src/impl/particles_impl_hskpng_vterm.ipp index 024581daa..d449e913e 100644 --- a/src/impl/particles_impl_hskpng_vterm.ipp +++ b/src/impl/particles_impl_hskpng_vterm.ipp @@ -174,7 +174,7 @@ namespace libcloudphxx zip_it_t(thrust::make_tuple( thrust::make_permutation_iterator(T.begin_ref(), ijk.begin_ref()), thrust::make_permutation_iterator(p.begin_ref(), ijk.begin_ref()), - thrust::make_permutation_iterator(rhod.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin_ref(), ijk.begin_ref()) )), // input - 2nd arg vt.begin(), // condition argument @@ -221,7 +221,7 @@ namespace libcloudphxx zip_it_t(thrust::make_tuple( thrust::make_permutation_iterator(T.begin_ref(), ijk.begin_ref()), thrust::make_permutation_iterator(p.begin_ref(), ijk.begin_ref()), - thrust::make_permutation_iterator(rhod.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin_ref(), ijk.begin_ref()) )), // input - 2nd arg vt.begin(), // output diff --git a/src/impl/particles_impl_sstp.ipp b/src/impl/particles_impl_sstp.ipp index 7bbc340f3..5970229f4 100644 --- a/src/impl/particles_impl_sstp.ipp +++ b/src/impl/particles_impl_sstp.ipp @@ -18,8 +18,11 @@ namespace libcloudphxx const int n = 4; thrust_device::vector - *fr[n] = { &rv, &th, &rhod, &p }, - *to[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh, &sstp_tmp_p }; + *fr[n] = { rv.ptr_ref(), th.ptr_ref(), &rhod, p.ptr_ref() }, + *to[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh, &sstp_tmp_p }; + thrust_device::vector + *_ijk[n] = { ijk.ptr_ref(), ijk.ptr_ref(), ijk.ptr(), ijk.ptr_ref() }; + for (int ix = 0; ix < ( (const_p && opts_init.exact_sstp_cond) ? n : n-1); ++ix) // TODO: var_rho { if(opts_init.exact_sstp_cond) // per-particle version @@ -68,8 +71,8 @@ namespace libcloudphxx const int n = 3; thrust_device::vector - *scl[n] = { &rv, &th, &rhod }, - *tmp[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh }; + *scl[n] = { rv.ptr_ref(), th.ptr_ref(), &rhod }, + *tmp[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh }; for (int ix = 0; ix < (var_rho ? n : n-1); ++ix) { @@ -115,9 +118,11 @@ namespace libcloudphxx const int n = 4; thrust_device::vector - *scl[n] = { &rv, &th, &rhod, &p }, - *tmp[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh, &sstp_tmp_p }, - *dlt[n] = { &tmp_device_real_part, &tmp_device_real_part1, &tmp_device_real_part2, &tmp_device_real_part5 }; + *scl[n] = { rv.ptr_ref(), th.ptr_ref(), &rhod, p.ptr_ref() }, + *tmp[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh, &sstp_tmp_p }, + *dlt[n] = { &tmp_device_real_part, &tmp_device_real_part1, &tmp_device_real_part2, &tmp_device_real_part5 }; + thrust_device::vector + *_ijk[n] = { ijk.ptr_ref(), ijk.ptr_ref(), ijk.ptr(), ijk.ptr_ref() }; for (int ix = 0; ix < (const_p ? n : n-1); ++ix) { @@ -126,16 +131,16 @@ namespace libcloudphxx { // sstp_tmp_scl = dscl_adv (i.e. delta, i.e. new - old) thrust::transform( - thrust::make_permutation_iterator(scl[ix]->begin(), ijk.begin()), // 1st arg: rv_new - thrust::make_permutation_iterator(scl[ix]->begin(), ijk.end()), // 1st arg: rv_new + thrust::make_permutation_iterator(scl[ix]->begin(), _ijk[ix]->begin()), // 1st arg: rv_new + thrust::make_permutation_iterator(scl[ix]->begin(), _ijk[ix]->end()), // 1st arg: rv_new tmp[ix]->begin(), // 2nd arg: rv_old dlt[ix]->begin(), // output (arg::_1 - arg::_2) / sstp // op: dscl_adv = (scl_new - scl_old) / sstp ); // scl -= (sstp - 1) * dscl_adv / sstp thrust::transform( - thrust::make_permutation_iterator(scl[ix]->begin(), ijk.begin()), // 1st arg: rv_new - thrust::make_permutation_iterator(scl[ix]->begin(), ijk.end()), // 1st arg: rv_new + thrust::make_permutation_iterator(scl[ix]->begin(), _ijk[ix]->begin()), // 1st arg: rv_new + thrust::make_permutation_iterator(scl[ix]->begin(), _ijk[ix]->end()), // 1st arg: rv_new dlt[ix]->begin(), // 2nd arg tmp[ix]->begin(), // output arg::_1 - (sstp - 1) * arg::_2 // op: rv = rv - (sstp - 1) * dscl_adv From cc4733475a2e6c0ea49aacc233709974c8ddd229 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 22 Sep 2022 13:20:15 +0200 Subject: [PATCH 26/38] ref grid: update th rv --- src/impl/particles_impl.ipp | 2 +- src/impl/particles_impl_ante_adding_SD.ipp | 6 +++--- src/impl/particles_impl_post_adding_SD.ipp | 6 +++--- src/impl/particles_impl_update_th_rv.ipp | 20 ++++++++++---------- 4 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 6980174ec..86e1b3392 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -614,7 +614,7 @@ namespace libcloudphxx 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_th_rv(ref_grid &); void update_state(thrust_device::vector &, thrust_device::vector &); void update_state(ref_grid &, thrust_device::vector &); void update_pstate(thrust_device::vector &, thrust_device::vector &); diff --git a/src/impl/particles_impl_ante_adding_SD.ipp b/src/impl/particles_impl_ante_adding_SD.ipp index e7b496a6e..9e5a603dd 100644 --- a/src/impl/particles_impl_ante_adding_SD.ipp +++ b/src/impl/particles_impl_ante_adding_SD.ipp @@ -14,8 +14,8 @@ namespace libcloudphxx { // --- calc liquid water content before src --- hskpng_sort(); - thrust_device::vector &drv(tmp_device_real_cell1); // NOTE: this can't be changed by any function called before a call to after_adding_SD... - thrust::fill(drv.begin(), drv.end(), real_t(0.)); + auto &drv(tmp_device_real_cell1); // NOTE: this can't be changed by any function called before a call to after_adding_SD... + thrust::fill(drv.begin_ref(), drv.end_ref(), real_t(0.)); moms_all(); moms_calc(rw2.begin(), real_t(3./2.)); @@ -23,7 +23,7 @@ namespace libcloudphxx // drv = - tot_vol_bfr thrust::transform( count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output + thrust::make_permutation_iterator(drv.begin_ref(), count_ijk.begin_ref()), // output thrust::negate() ); diff --git a/src/impl/particles_impl_post_adding_SD.ipp b/src/impl/particles_impl_post_adding_SD.ipp index 675c58485..e07c8825c 100644 --- a/src/impl/particles_impl_post_adding_SD.ipp +++ b/src/impl/particles_impl_post_adding_SD.ipp @@ -12,7 +12,7 @@ namespace libcloudphxx template void particles_t::impl::post_adding_SD() { - thrust_device::vector &drv(tmp_device_real_cell1); // NOTE: this can't be changed by any function called since ante_adding_SD() was called..... + auto &drv(tmp_device_real_cell1); // NOTE: this can't be changed by any function called since ante_adding_SD() was called..... // --- after source particles are no longer sorted --- sorted = false; @@ -25,8 +25,8 @@ namespace libcloudphxx // drv = tot_vol_after -tot_vol_bfr + dry_vol_bfr thrust::transform( count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // 2nd arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output + thrust::make_permutation_iterator(drv.begin_ref(), count_ijk.begin_ref()), // 2nd arg + thrust::make_permutation_iterator(drv.begin_ref(), count_ijk.begin_ref()), // output thrust::plus() ); diff --git a/src/impl/particles_impl_update_th_rv.ipp b/src/impl/particles_impl_update_th_rv.ipp index 3b5241a96..f9f963215 100644 --- a/src/impl/particles_impl_update_th_rv.ipp +++ b/src/impl/particles_impl_update_th_rv.ipp @@ -35,29 +35,29 @@ namespace libcloudphxx // particles have to be sorted template void particles_t::impl::update_th_rv( - thrust_device::vector &drv // change in water vapor mixing ratio + ref_grid &drv // change in water vapor mixing ratio ) { if(!sorted) throw std::runtime_error("update_th_rv called on an unsorted set"); - nancheck(drv, "update_th_rv: input drv"); + nancheck(drv.get_ref(), "update_th_rv: input drv"); // multiplying specific 3rd moms diff by -rho_w*4/3*pi thrust::transform( - drv.begin(), drv.end(), // input - 1st arg + drv.begin_ref(), drv.end_ref(), // input - 1st arg thrust::make_constant_iterator( // input - 2nd arg - common::moist_air::rho_w() / si::kilograms * si::cubic_metres * real_t(4./3) * pi() ), - drv.begin(), // output + drv.begin_ref(), // output thrust::multiplies() ); // updating rv assert(*thrust::min_element(rv.begin(), rv.end()) >= 0); thrust::transform( - rv.begin(), rv.end(), // input - 1st arg - drv.begin(), // input - 2nd arg - rv.begin(), // output + rv.begin_ref(), rv.end_ref(), // input - 1st arg + drv.begin_ref(), // input - 2nd arg + rv.begin_ref(), // output thrust::plus() ); assert(*thrust::min_element(rv.begin(), rv.end()) >= 0); @@ -76,9 +76,9 @@ namespace libcloudphxx th.begin(), th.end(), // input - 1st arg thrust::make_transform_iterator( zip_it_t(thrust::make_tuple( - drv.begin(), // - T.begin(), // dth = drv * d_th_d_rv(T, th) - th.begin() // + drv.begin_ref(), // + T.begin_ref(), // dth = drv * d_th_d_rv(T, th) + th.begin_ref() // )), detail::dth() ), From 25574a1d977235695107182b5bd9267230cd97f6 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 22 Sep 2022 13:56:25 +0200 Subject: [PATCH 27/38] ref grid: wip --- src/impl/particles_impl.ipp | 12 +++++++----- src/impl/particles_impl_chem_henry.ipp | 2 +- src/impl/particles_impl_coal.ipp | 2 +- src/impl/particles_impl_cond.ipp | 10 +++++----- src/impl/particles_impl_hskpng_resize.ipp | 2 +- src/impl/particles_impl_hskpng_turb_ss.ipp | 3 ++- src/impl/particles_impl_hskpng_turb_vel.ipp | 2 +- src/impl/particles_impl_init_count_num.ipp | 2 +- src/impl/particles_impl_init_hskpng_ncell.ipp | 2 +- src/impl/particles_impl_init_n.ipp | 2 +- src/impl/particles_impl_mass_dens.ipp | 2 +- src/impl/particles_impl_update_th_rv.ipp | 10 +++++----- 12 files changed, 27 insertions(+), 24 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 86e1b3392..0d22e8a50 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -106,9 +106,8 @@ namespace libcloudphxx // housekeeping data (per particle) thrust_device::vector - i, j, k, // Eulerian grid cell indices (always zero for 0D) - sorted_id, sorted_ijk; - ref_part ijk; // ijk in the normal and in the refined grid + i, j, k; // Eulerian grid cell indices (always zero for 0D) + ref_part ijk, sorted_id, sorted_ijk; // ijk in the normal and in the refined grid ref_grid ijk_ref2ijk; // maps refined cell index to the index of the normal cell containing this refined cell // helpers that reference refined arrays if refinement is done and non-refined otherwise @@ -240,9 +239,9 @@ namespace libcloudphxx tmp_device_real_part3, tmp_device_real_part4, tmp_device_real_part5, - tmp_device_real_cell, &u01; // uniform random numbers between 0 and 1 // TODO: use the tmp array as rand argument? ref_grid + tmp_device_real_cell, tmp_device_real_cell1, tmp_device_real_cell2; ref_grid @@ -340,6 +339,7 @@ namespace libcloudphxx count_num(n_cell), count_mom(n_cell), count_n(0,0), + tmp_device_real_cell(n_cell), tmp_device_real_cell1(n_cell), tmp_device_real_cell2(n_cell), tmp_host_real_cell(n_cell), @@ -390,7 +390,9 @@ namespace libcloudphxx sstp_tmp_rh(_opts_init.exact_sstp_cond ? sstp_pp_tmp_rh : sstp_pc_tmp_rh.get_ref()), sstp_tmp_p(sstp_pp_tmp_p), // not used in per-cell substepping pure_const_multi (((_opts_init.sd_conc) == 0) && (_opts_init.sd_const_multi > 0 || _opts_init.dry_sizes.size() > 0)), // coal prob can be greater than one only in sd_conc simulations - ijk(opts_init.n_ref) + ijk(opts_init.n_ref), + sorted_ijk(opts_init.n_ref), + sorted_id(opts_init.n_ref) { // set 0 dev_count to mark that its not a multi_CUDA spawn diff --git a/src/impl/particles_impl_chem_henry.ipp b/src/impl/particles_impl_chem_henry.ipp index 297a634d8..b4f5e8d79 100644 --- a/src/impl/particles_impl_chem_henry.ipp +++ b/src/impl/particles_impl_chem_henry.ipp @@ -307,7 +307,7 @@ namespace libcloudphxx hskpng_sort(); // temporarily needed to store old mass per cell - thrust_device::vector &mass_old(tmp_device_real_cell); + thrust_device::vector &mass_old(tmp_device_real_cell.get()); thrust_device::vector &mass_new(tmp_device_real_cell1.get()); for (int i = 0; i < chem_gas_n; ++i) diff --git a/src/impl/particles_impl_coal.ipp b/src/impl/particles_impl_coal.ipp index 67e106b03..0a2c6c36b 100644 --- a/src/impl/particles_impl_coal.ipp +++ b/src/impl/particles_impl_coal.ipp @@ -273,7 +273,7 @@ namespace libcloudphxx // references to tmp data thrust_device::vector - &scl(tmp_device_real_cell), // scale factor for probablility + &scl(tmp_device_real_cell.get()), // 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 // 1st one of a pair stores number of collisions, 2nd one stores info on which one has greater multiplicity thrust_device::vector diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index d2ea5dfad..d77017b2f 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -24,7 +24,7 @@ namespace libcloudphxx // --- calc liquid water content before cond --- hskpng_sort(); - thrust_device::vector &drv(tmp_device_real_cell); + ref_grid &drv(tmp_device_real_cell); // calculating the 3rd wet moment before condensation moms_all(); @@ -33,11 +33,11 @@ namespace libcloudphxx // permute-copying the result to -dm_3 // fill with 0s if not all cells will be updated in the following transform - if(count_n.get() != n_cell.get()) thrust::fill(drv.begin(), drv.end(), real_t(0.)); + if(count_n.get() != n_cell.get()) thrust::fill(drv.begin_ref(), drv.end_ref(), real_t(0.)); thrust::transform( count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output + thrust::make_permutation_iterator(drv.begin_ref(), count_ijk.begin_ref()), // output thrust::negate() ); @@ -100,8 +100,8 @@ namespace libcloudphxx // adding the third moment after condensation to dm_3 thrust::transform( count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // input - 2nd arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output + thrust::make_permutation_iterator(drv.begin_ref(), count_ijk.begin_ref()), // input - 2nd arg + thrust::make_permutation_iterator(drv.begin_ref(), count_ijk.begin_ref()), // output thrust::plus() ); diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index 757b58791..121e3810f 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -15,7 +15,7 @@ namespace libcloudphxx } } { - thrust_device::vector *vec[] = {&sorted_id, &sorted_ijk}; + ref_part *vec[] = {&sorted_id, &sorted_ijk}; for(int i=0; i<3; ++i) { vec[i]->resize(n_part); diff --git a/src/impl/particles_impl_hskpng_turb_ss.ipp b/src/impl/particles_impl_hskpng_turb_ss.ipp index 735dc91b8..1a44b4e91 100644 --- a/src/impl/particles_impl_hskpng_turb_ss.ipp +++ b/src/impl/particles_impl_hskpng_turb_ss.ipp @@ -41,10 +41,11 @@ namespace libcloudphxx }; // calc the SGS turbulent supersaturation + // TODO: do it on refined grid template void particles_t::impl::hskpng_turb_dot_ss() { - thrust_device::vector &tau_rlx(tmp_device_real_cell); // tau_rlx needs to have length n_cell + thrust_device::vector &tau_rlx(tmp_device_real_cell.get()); // 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); diff --git a/src/impl/particles_impl_hskpng_turb_vel.ipp b/src/impl/particles_impl_hskpng_turb_vel.ipp index e20221a9e..32ee95b24 100644 --- a/src/impl/particles_impl_hskpng_turb_vel.ipp +++ b/src/impl/particles_impl_hskpng_turb_vel.ipp @@ -53,7 +53,7 @@ namespace libcloudphxx { namespace arg = thrust::placeholders; - thrust_device::vector &tau(tmp_device_real_cell); + thrust_device::vector &tau(tmp_device_real_cell.get()); thrust_device::vector &tke(diss_rate); // should be called after hskpng_tke, which replaces diss_rate with tke thrust::transform(tke.begin(), tke.end(), diff --git a/src/impl/particles_impl_init_count_num.ipp b/src/impl/particles_impl_init_count_num.ipp index a1ad339bc..c135ed43e 100644 --- a/src/impl/particles_impl_init_count_num.ipp +++ b/src/impl/particles_impl_init_count_num.ipp @@ -65,7 +65,7 @@ namespace libcloudphxx template void particles_t::impl::init_count_num_hlpr(const real_t &conc, const thrust_size_t &const_multi) { - thrust_device::vector &concentration(tmp_device_real_cell); + thrust_device::vector &concentration(tmp_device_real_cell.get()); thrust::fill(concentration.begin(), concentration.end(), conc); conc_to_number(concentration); diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index 31312507e..8c515316d 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -24,7 +24,7 @@ namespace libcloudphxx // count_n.get_ref() = 0; // initialising device temporary arrays - tmp_device_real_cell.resize(n_cell.get()); +// tmp_device_real_cell.resize(n_cell.get()); // tmp_device_real_cell1.resize(n_cell); // tmp_device_real_cell2.resize(n_cell); tmp_device_size_cell.resize(n_cell.get()); diff --git a/src/impl/particles_impl_init_n.ipp b/src/impl/particles_impl_init_n.ipp index 74ae9e91f..2aa6f18e8 100644 --- a/src/impl/particles_impl_init_n.ipp +++ b/src/impl/particles_impl_init_n.ipp @@ -135,7 +135,7 @@ namespace libcloudphxx void particles_t::impl::init_n_dry_sizes(const real_t &conc, const thrust_size_t &sd_conc) { namespace arg = thrust::placeholders; - thrust_device::vector &concentration(tmp_device_real_cell); + thrust_device::vector &concentration(tmp_device_real_cell.get()); thrust::fill(concentration.begin(), concentration.end(), conc); conc_to_number(concentration); thrust::transform( diff --git a/src/impl/particles_impl_mass_dens.ipp b/src/impl/particles_impl_mass_dens.ipp index bf3736fff..a57247a21 100644 --- a/src/impl/particles_impl_mass_dens.ipp +++ b/src/impl/particles_impl_mass_dens.ipp @@ -39,7 +39,7 @@ namespace libcloudphxx thrust_device::vector &n_filtered(tmp_device_real_part); // number of SD in each cell casted to real_t - thrust_device::vector &count_num_real_t(tmp_device_real_cell); + thrust_device::vector &count_num_real_t(tmp_device_real_cell.get()); // get number of SD in each cell hskpng_count(); diff --git a/src/impl/particles_impl_update_th_rv.ipp b/src/impl/particles_impl_update_th_rv.ipp index f9f963215..773e16fdc 100644 --- a/src/impl/particles_impl_update_th_rv.ipp +++ b/src/impl/particles_impl_update_th_rv.ipp @@ -100,9 +100,9 @@ namespace libcloudphxx if(!sorted) throw std::runtime_error("update_uh_rv called on an unsorted set"); // cell-wise change in state - thrust_device::vector &dstate(tmp_device_real_cell); + ref_grid &dstate(tmp_device_real_cell); // init dstate with 0s - thrust::fill(dstate.begin(), dstate.end(), real_t(0)); + thrust::fill(dstate.begin_ref(), dstate.end_ref(), real_t(0)); // calc sum of pdstate in each cell thrust::pair< thrust_device::vector::iterator, @@ -118,15 +118,15 @@ namespace libcloudphxx // add this sum to dstate thrust::transform( count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg - thrust::make_permutation_iterator(dstate.begin(), count_ijk.begin()), // 2nd arg - thrust::make_permutation_iterator(dstate.begin(), count_ijk.begin()), // output + thrust::make_permutation_iterator(dstate.begin_ref(), count_ijk.begin_ref()), // 2nd arg + thrust::make_permutation_iterator(dstate.begin_ref(), count_ijk.begin_ref()), // output thrust::plus() ); // add dstate to pstate thrust::transform( pstate.begin(), pstate.end(), - thrust::make_permutation_iterator(dstate.begin(), ijk.begin()), + thrust::make_permutation_iterator(dstate.begin_ref(), ijk.begin_ref()), pstate.begin(), thrust::plus() ); From 315c6af59aae010403683871be91355948bcf0f2 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 22 Sep 2022 14:53:45 +0200 Subject: [PATCH 28/38] ref wip --- src/detail/debug.hpp | 35 ++++++++++++++++++++++++++++++++ src/detail/thrust.hpp | 19 ----------------- src/impl/particles_impl_rcyc.ipp | 6 +++--- src/particles.tpp | 1 + 4 files changed, 39 insertions(+), 22 deletions(-) create mode 100644 src/detail/debug.hpp diff --git a/src/detail/debug.hpp b/src/detail/debug.hpp new file mode 100644 index 000000000..aad726189 --- /dev/null +++ b/src/detail/debug.hpp @@ -0,0 +1,35 @@ +#pragma once + +namespace libcloudphxx +{ + namespace lgrngn + { +//#if !defined(NDEBUG) // TODO (CMake defaults) + namespace debug + { + template + void print(const vec_t &v) + { + thrust::copy(v.begin(), v.end(), std::ostream_iterator(std::cerr, "\t")); + std::cerr << std::endl; + } + + template + void print(const vec_bgn_t &bgn, const vec_end_t &end) + { + thrust::copy(bgn, end, std::ostream_iterator(std::cerr, "\t")); + std::cerr << std::endl; + } + + template + void print(ref_grid v) + { + std::cerr << "normal grid values:" << std::endl; + print(v.begin(), v.end()); + std::cerr << "refined grid values:" << std::endl; + print(v.get_ref()); + } + }; +//#endif + }; +}; diff --git a/src/detail/thrust.hpp b/src/detail/thrust.hpp index 15e9468dc..d5c5c89c4 100644 --- a/src/detail/thrust.hpp +++ b/src/detail/thrust.hpp @@ -9,24 +9,5 @@ namespace libcloudphxx namespace lgrngn { typedef thrust_device::vector::size_type thrust_size_t; - -//#if !defined(NDEBUG) // TODO (CMake defaults) - namespace debug - { - template - void print(const vec_t &v) - { - thrust::copy(v.begin(), v.end(), std::ostream_iterator(std::cerr, "\t")); - std::cerr << std::endl; - } - - template - void print(const vec_bgn_t &bgn, const vec_end_t &end) - { - thrust::copy(bgn, end, std::ostream_iterator(std::cerr, "\t")); - std::cerr << std::endl; - } - }; -//#endif }; }; diff --git a/src/impl/particles_impl_rcyc.ipp b/src/impl/particles_impl_rcyc.ipp index 053cf559a..4a5a05de5 100644 --- a/src/impl/particles_impl_rcyc.ipp +++ b/src/impl/particles_impl_rcyc.ipp @@ -71,7 +71,7 @@ namespace libcloudphxx #else static_assert(sizeof(thrust_size_t) == sizeof(n_t), ""); #endif - thrust_device::vector &tmp(sorted_ijk); + thrust_device::vector &tmp(sorted_ijk.get()); thrust::copy(n.begin(), n.end(), tmp.begin()); thrust::sort_by_key( @@ -98,12 +98,12 @@ namespace libcloudphxx // for each property... for(auto vec: distmem_real_vctrs) - detail::copy_prop(vec->begin(), sorted_id, n_flagged); + detail::copy_prop(vec->begin(), sorted_id.get(), n_flagged); // ... chemical properties only if chem enabled if (opts_init.chem_switch){ for (int i = 0; i < chem_all; ++i) - detail::copy_prop(chem_bgn[i], sorted_id, n_flagged); + detail::copy_prop(chem_bgn[i], sorted_id.get(), n_flagged); } { diff --git a/src/particles.tpp b/src/particles.tpp index ad1327f9f..2e048dc71 100644 --- a/src/particles.tpp +++ b/src/particles.tpp @@ -29,6 +29,7 @@ #include "detail/functors_host.hpp" #include "detail/ran_with_mpi.hpp" #include "detail/grid_refinement_helpers.hpp" +#include "detail/debug.hpp" //kernel definitions #include "detail/kernel_definitions/hall_efficiencies.hpp" From 94fad7892c5ebf46892e6a63e5d982b7b3222445 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 22 Sep 2022 15:07:10 +0200 Subject: [PATCH 29/38] ref grid: debug output --- src/detail/debug.hpp | 13 +++++++++++-- src/detail/grid_refinement_helpers.hpp | 20 ++++++++++++++++++++ 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/src/detail/debug.hpp b/src/detail/debug.hpp index aad726189..567cca741 100644 --- a/src/detail/debug.hpp +++ b/src/detail/debug.hpp @@ -25,9 +25,18 @@ namespace libcloudphxx void print(ref_grid v) { std::cerr << "normal grid values:" << std::endl; - print(v.begin(), v.end()); + print(v.cbegin(), v.cend()); std::cerr << "refined grid values:" << std::endl; - print(v.get_ref()); + print(v.cbegin_ref(), v.cend_ref()); + } + + template + void print(ref_part v) + { + std::cerr << "particle parameters on normal grid:" << std::endl; + print(v.cbegin(), v.cend()); + std::cerr << "particle parameters on refined grid:" << std::endl; + print(v.cbegin_ref(), v.cend_ref()); } }; //#endif diff --git a/src/detail/grid_refinement_helpers.hpp b/src/detail/grid_refinement_helpers.hpp index de0c71a2e..b71a8eb6a 100644 --- a/src/detail/grid_refinement_helpers.hpp +++ b/src/detail/grid_refinement_helpers.hpp @@ -102,21 +102,41 @@ namespace libcloudphxx return arr.begin(); } + const auto cbegin() + { + return arr.begin(); + } + auto begin_ref() { return this->ref_flag ? arr_ref.begin() : arr.begin(); } + auto cbegin_ref() + { + return this->ref_flag ? arr_ref.begin() : arr.begin(); + } + auto end() { return arr.end(); } + const auto cend() + { + return arr.end(); + } + auto end_ref() { return this->ref_flag ? arr_ref.end() : arr.end(); } + auto cend_ref() + { + return this->ref_flag ? arr_ref.end() : arr.end(); + } + auto ptr() { return &arr; From 319e59a96d573edff3846462b539e881d0f50fb2 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 22 Sep 2022 15:20:34 +0200 Subject: [PATCH 30/38] ref grid: fixed comp errors --- src/impl/particles_impl_cond_sstp.ipp | 8 ++++---- src/impl/particles_impl_hskpng_Tpr.ipp | 18 +++++++++--------- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/impl/particles_impl_cond_sstp.ipp b/src/impl/particles_impl_cond_sstp.ipp index 01c145724..fdfcef427 100644 --- a/src/impl/particles_impl_cond_sstp.ipp +++ b/src/impl/particles_impl_cond_sstp.ipp @@ -38,8 +38,8 @@ namespace libcloudphxx const RH_iter &rhi ) { - thrust_device::vector &lambda_D(tmp_device_real_cell1); // real_cell used in cond.ipp - thrust_device::vector &lambda_K(tmp_device_real_cell2); // real_cell used in cond.ipp + auto &lambda_D(tmp_device_real_cell1); // real_cell used in cond.ipp + auto &lambda_K(tmp_device_real_cell2); // real_cell used in cond.ipp auto hlpr_zip_iter = thrust::make_zip_iterator(thrust::make_tuple( sstp_tmp_rh.begin(), @@ -53,8 +53,8 @@ namespace libcloudphxx rd3.begin(), kpa.begin(), vt.begin(), - thrust::make_permutation_iterator(lambda_D.begin(), ijk.begin()), - thrust::make_permutation_iterator(lambda_K.begin(), ijk.begin()) + thrust::make_permutation_iterator(lambda_D.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(lambda_K.begin_ref(), ijk.begin_ref()) )); // calculating drop growth in a timestep using backward Euler diff --git a/src/impl/particles_impl_hskpng_Tpr.ipp b/src/impl/particles_impl_hskpng_Tpr.ipp index a80f11a69..f97286b08 100644 --- a/src/impl/particles_impl_hskpng_Tpr.ipp +++ b/src/impl/particles_impl_hskpng_Tpr.ipp @@ -195,18 +195,10 @@ namespace libcloudphxx } { - typedef thrust::zip_iterator< - thrust::tuple< - typename thrust_device::vector::iterator, - typename thrust_device::vector::iterator, - typename thrust_device::vector::iterator - > - > zip_it_t; - if(!const_p) { // p = common::theta_dry::p(rhod, r, T); - zip_it_t it(thrust::make_tuple( + auto it = thrust::make_zip_iterator(thrust::make_tuple( thrust::make_permutation_iterator(rhod.begin(), ijk_ref2ijk.begin()), rv.begin_ref(), T.begin_ref() @@ -219,6 +211,14 @@ namespace libcloudphxx ); } + typedef thrust::zip_iterator< + thrust::tuple< + typename thrust_device::vector::iterator, + typename thrust_device::vector::iterator, + typename thrust_device::vector::iterator + > + > zip_it_t; + thrust::transform( zip_it_t(thrust::make_tuple(p.begin_ref(), rv.begin_ref(), T.begin_ref())), // input - begin zip_it_t(thrust::make_tuple(p.end_ref(), rv.end_ref(), T.end_ref() )), // input - end From 8056761515de69934b333309dc9164b61e239ee7 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 26 Sep 2022 12:55:21 +0200 Subject: [PATCH 31/38] api_lgrngn debug --- src/impl/particles_impl_hskpng_Tpr.ipp | 4 ++-- src/impl/particles_impl_hskpng_resize.ipp | 3 +-- src/particles_step.ipp | 12 ++++++------ 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/src/impl/particles_impl_hskpng_Tpr.ipp b/src/impl/particles_impl_hskpng_Tpr.ipp index f97286b08..3e8818dff 100644 --- a/src/impl/particles_impl_hskpng_Tpr.ipp +++ b/src/impl/particles_impl_hskpng_Tpr.ipp @@ -178,7 +178,7 @@ namespace libcloudphxx // T = common::theta_dry::T(th, rhod); thrust::transform( th.begin_ref(), th.end_ref(), // input - first arg - thrust::make_permutation_iterator(rhod.begin(), ijk_ref2ijk.begin()), // input - second arg + thrust::make_permutation_iterator(rhod.begin(), ijk_ref2ijk.begin_ref()), // input - second arg T.begin_ref(), // output detail::common__theta_dry__T_rhod() ); @@ -199,7 +199,7 @@ namespace libcloudphxx { // p = common::theta_dry::p(rhod, r, T); auto it = thrust::make_zip_iterator(thrust::make_tuple( - thrust::make_permutation_iterator(rhod.begin(), ijk_ref2ijk.begin()), + thrust::make_permutation_iterator(rhod.begin(), ijk_ref2ijk.begin_ref()), rv.begin_ref(), T.begin_ref() )); diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index 121e3810f..914bce5b3 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -15,7 +15,7 @@ namespace libcloudphxx } } { - ref_part *vec[] = {&sorted_id, &sorted_ijk}; + ref_part *vec[] = {&ijk, &sorted_id, &sorted_ijk}; for(int i=0; i<3; ++i) { vec[i]->resize(n_part); @@ -24,7 +24,6 @@ namespace libcloudphxx n.resize(n_part); tmp_device_n_part.resize(n_part); tmp_device_size_part.resize(n_part); - ijk.resize(n_part); vt.resize(n_part, detail::invalid); diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 50de407e7..055ad7a1b 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -104,8 +104,8 @@ namespace libcloudphxx // fill in mpi courant halos pimpl->xchng_courants(); - nancheck(pimpl->th, " th after sync-in"); - nancheck(pimpl->rv, " rv after sync-in"); + nancheck(pimpl->th.get_ref(), " th after sync-in"); + nancheck(pimpl->rv.get_ref(), " rv after sync-in"); 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"); @@ -114,8 +114,8 @@ namespace libcloudphxx 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->rv.begin_ref(), pimpl->rv.end_ref()) >= 0); + assert(*thrust::min_element(pimpl->th.begin_ref(), pimpl->th.end_ref()) >= 0); assert(*thrust::min_element(pimpl->rhod.begin(), pimpl->rhod.end()) >= 0); if(pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_switch || pimpl->opts_init.turb_coal_switch) {assert(*thrust::min_element(pimpl->diss_rate.begin(), pimpl->diss_rate.end()) >= 0);} @@ -206,8 +206,8 @@ namespace libcloudphxx } } - nancheck(pimpl->th, " th after cond"); - nancheck(pimpl->rv, " rv after cond"); + nancheck(pimpl->th.get_ref(), " th after cond"); + nancheck(pimpl->rv.get_ref(), " rv after cond"); // saving rv to be used as rv_old pimpl->sstp_save(); From 707062403f4fa9c48c169a94c54a475f92f12f52 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 26 Sep 2022 13:25:36 +0200 Subject: [PATCH 32/38] ref wip --- src/impl/particles_impl_update_th_rv.ipp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/impl/particles_impl_update_th_rv.ipp b/src/impl/particles_impl_update_th_rv.ipp index 773e16fdc..9e16d5425 100644 --- a/src/impl/particles_impl_update_th_rv.ipp +++ b/src/impl/particles_impl_update_th_rv.ipp @@ -89,7 +89,7 @@ namespace libcloudphxx nancheck(th, "update_th_rv: th after update"); } - // update particle-specific cell state + // update particle-specific cell state; done on refined grid // particles have to be sorted template void particles_t::impl::update_pstate( @@ -97,7 +97,7 @@ namespace libcloudphxx thrust_device::vector &pdstate // change in cell characteristic ) { - if(!sorted) throw std::runtime_error("update_uh_rv called on an unsorted set"); + if(!sorted) throw std::runtime_error("update_th_rv called on an unsorted set"); // cell-wise change in state ref_grid &dstate(tmp_device_real_cell); @@ -108,16 +108,16 @@ namespace libcloudphxx thrust_device::vector::iterator, typename thrust_device::vector::iterator > it_pair = thrust::reduce_by_key( - sorted_ijk.begin(), sorted_ijk.end(), - thrust::make_permutation_iterator(pdstate.begin(), sorted_id.begin()), - count_ijk.begin(), - count_mom.begin() + sorted_ijk.begin_ref(), sorted_ijk.end_ref(), + thrust::make_permutation_iterator(pdstate.begin(), sorted_id.begin_ref()), + count_ijk.begin_ref(), + count_mom.begin_ref() ); - count_n.get() = it_pair.first - count_ijk.begin(); + count_n.get_ref() = it_pair.first - count_ijk.begin_ref(); // add this sum to dstate thrust::transform( - count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg + count_mom.begin_ref(), count_mom.begin_ref() + count_n.get_ref(), // input - 1st arg thrust::make_permutation_iterator(dstate.begin_ref(), count_ijk.begin_ref()), // 2nd arg thrust::make_permutation_iterator(dstate.begin_ref(), count_ijk.begin_ref()), // output thrust::plus() From e8ce178d7bc6a3ed5b4d9a382aba6b56e5180f19 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 26 Sep 2022 13:25:55 +0200 Subject: [PATCH 33/38] ref wip --- src/impl/particles_impl_hskpng_sort.ipp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/impl/particles_impl_hskpng_sort.ipp b/src/impl/particles_impl_hskpng_sort.ipp index 4f4abcd55..7b4f252ab 100644 --- a/src/impl/particles_impl_hskpng_sort.ipp +++ b/src/impl/particles_impl_hskpng_sort.ipp @@ -19,11 +19,11 @@ namespace libcloudphxx if (!shuffle) { - // making a copy of ijk - thrust::copy( - ijk.begin(), ijk.end(), // from - sorted_ijk.begin() // to - ); + // making a copy of ijk + thrust::copy( + ijk.begin(), ijk.end(), // from + sorted_ijk.begin() // to + ); } else { @@ -46,8 +46,8 @@ namespace libcloudphxx // sorting sorted_ijk and sorted_id thrust::sort_by_key( - sorted_ijk.begin(), sorted_ijk.end(), // keys - sorted_id.begin() // values + sorted_ijk.begin(), sorted_ijk.end(), // keys + sorted_id.begin() // values ); // flagging that particles are now sorted From 2c23a9d81e8a4b10bccfa98474f7d95f2c3abf46 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 26 Sep 2022 14:01:51 +0200 Subject: [PATCH 34/38] ref grid: fixed comp errors --- src/impl/particles_impl.ipp | 7 +++---- src/impl/particles_impl_hskpng_resize.ipp | 5 +++-- src/impl/particles_impl_hskpng_sort.ipp | 24 +++++++++++++++++------ src/impl/particles_impl_rcyc.ipp | 4 ++-- src/impl/particles_impl_update_th_rv.ipp | 2 +- 5 files changed, 27 insertions(+), 15 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 0d22e8a50..620fa83bd 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -106,8 +106,8 @@ namespace libcloudphxx // housekeeping data (per particle) thrust_device::vector - i, j, k; // Eulerian grid cell indices (always zero for 0D) - ref_part ijk, sorted_id, sorted_ijk; // ijk in the normal and in the refined grid + i, j, k, sorted_id; // Eulerian grid cell indices (always zero for 0D), one sorted_id needed, because its the same in normal and refined grids as we always sort them together + ref_part ijk, sorted_ijk; // ijk in the normal and in the refined grid ref_grid ijk_ref2ijk; // maps refined cell index to the index of the normal cell containing this refined cell // helpers that reference refined arrays if refinement is done and non-refined otherwise @@ -391,8 +391,7 @@ namespace libcloudphxx sstp_tmp_p(sstp_pp_tmp_p), // not used in per-cell substepping pure_const_multi (((_opts_init.sd_conc) == 0) && (_opts_init.sd_const_multi > 0 || _opts_init.dry_sizes.size() > 0)), // coal prob can be greater than one only in sd_conc simulations ijk(opts_init.n_ref), - sorted_ijk(opts_init.n_ref), - sorted_id(opts_init.n_ref) + sorted_ijk(opts_init.n_ref) { // set 0 dev_count to mark that its not a multi_CUDA spawn diff --git a/src/impl/particles_impl_hskpng_resize.ipp b/src/impl/particles_impl_hskpng_resize.ipp index 914bce5b3..7311688ed 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -15,8 +15,8 @@ namespace libcloudphxx } } { - ref_part *vec[] = {&ijk, &sorted_id, &sorted_ijk}; - for(int i=0; i<3; ++i) + ref_part *vec[] = {&ijk, &sorted_ijk}; + for(int i=0; i<2; ++i) { vec[i]->resize(n_part); } @@ -24,6 +24,7 @@ namespace libcloudphxx n.resize(n_part); tmp_device_n_part.resize(n_part); tmp_device_size_part.resize(n_part); + sorted_id.resize(n_part); vt.resize(n_part, detail::invalid); diff --git a/src/impl/particles_impl_hskpng_sort.ipp b/src/impl/particles_impl_hskpng_sort.ipp index 7b4f252ab..0a0e15f15 100644 --- a/src/impl/particles_impl_hskpng_sort.ipp +++ b/src/impl/particles_impl_hskpng_sort.ipp @@ -20,10 +20,15 @@ namespace libcloudphxx if (!shuffle) { // making a copy of ijk +// thrust_device::vector *ijk[] = {&ijk.get(), &ijk.get_ref()}; thrust::copy( ijk.begin(), ijk.end(), // from sorted_ijk.begin() // to ); + thrust::copy( + ijk.begin_ref(), ijk.end_ref(), // from + sorted_ijk.begin_ref() // to + ); } else { @@ -36,18 +41,25 @@ namespace libcloudphxx sorted_id.begin() ); - // permuting sorted_ijk accordingly + // permuting sorted_ijk ref accordingly thrust::copy( - thrust::make_permutation_iterator(ijk.begin(), sorted_id.begin()), // input - begin - thrust::make_permutation_iterator(ijk.end(), sorted_id.end() ), // input - end - sorted_ijk.begin() // output + thrust::make_permutation_iterator(ijk.begin_ref(), sorted_id.begin()), // input - begin + thrust::make_permutation_iterator(ijk.end_ref(), sorted_id.end() ), // input - end + sorted_ijk.begin_ref() // output ); } // sorting sorted_ijk and sorted_id thrust::sort_by_key( - sorted_ijk.begin(), sorted_ijk.end(), // keys - sorted_id.begin() // values + sorted_ijk.begin_ref(), sorted_ijk.end_ref(), // keys + sorted_id.begin() // values + ); + + // set sroted_ijk for normal grid + thrust::copy( + thrust::make_permutation_iterator(ijk.begin(), sorted_id.begin()), + thrust::make_permutation_iterator(ijk.begin(), sorted_id.begin()) + n_part, + sorted_ijk.begin() ); // flagging that particles are now sorted diff --git a/src/impl/particles_impl_rcyc.ipp b/src/impl/particles_impl_rcyc.ipp index 4a5a05de5..b2d9e23a3 100644 --- a/src/impl/particles_impl_rcyc.ipp +++ b/src/impl/particles_impl_rcyc.ipp @@ -98,12 +98,12 @@ namespace libcloudphxx // for each property... for(auto vec: distmem_real_vctrs) - detail::copy_prop(vec->begin(), sorted_id.get(), n_flagged); + detail::copy_prop(vec->begin(), sorted_id, n_flagged); // ... chemical properties only if chem enabled if (opts_init.chem_switch){ for (int i = 0; i < chem_all; ++i) - detail::copy_prop(chem_bgn[i], sorted_id.get(), n_flagged); + detail::copy_prop(chem_bgn[i], sorted_id, n_flagged); } { diff --git a/src/impl/particles_impl_update_th_rv.ipp b/src/impl/particles_impl_update_th_rv.ipp index 9e16d5425..2caad0baa 100644 --- a/src/impl/particles_impl_update_th_rv.ipp +++ b/src/impl/particles_impl_update_th_rv.ipp @@ -109,7 +109,7 @@ namespace libcloudphxx typename thrust_device::vector::iterator > it_pair = thrust::reduce_by_key( sorted_ijk.begin_ref(), sorted_ijk.end_ref(), - thrust::make_permutation_iterator(pdstate.begin(), sorted_id.begin_ref()), + thrust::make_permutation_iterator(pdstate.begin(), sorted_id.begin()), count_ijk.begin_ref(), count_mom.begin_ref() ); From aa20e1d480999a35ce07b3931383f68d166fdc5b Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 26 Sep 2022 14:50:50 +0200 Subject: [PATCH 35/38] ref grid: fixed comp errors --- src/impl/particles_impl.ipp | 2 +- src/impl/particles_impl_coal.ipp | 4 ++-- src/impl/particles_impl_hskpng_vterm.ipp | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 620fa83bd..2ff9c1b93 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -356,7 +356,7 @@ namespace libcloudphxx rng(_opts_init.rng_seed), src_stp_ctr(0), rlx_stp_ctr(0), - bcond(bcond), + bcond(bcond), n_x_bfr(0), n_cell_bfr(0), mpi_rank(mpi_rank), diff --git a/src/impl/particles_impl_coal.ipp b/src/impl/particles_impl_coal.ipp index 0a2c6c36b..de5b7b3c7 100644 --- a/src/impl/particles_impl_coal.ipp +++ b/src/impl/particles_impl_coal.ipp @@ -375,7 +375,7 @@ 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_ref(), sorted_ijk.begin_ref()), // TODO: use eta.begin(), but right no we dont have et on normal grid // tke dissipation rate thrust::make_permutation_iterator(thrust::make_constant_iterator(0), sorted_ijk.begin()) ) @@ -388,7 +388,7 @@ 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_ref(), sorted_ijk.begin_ref()), // TODO: use eta.begin(), but right no we dont have et on normal grid // tke dissipation rate thrust::make_permutation_iterator(diss_rate.begin(), sorted_ijk.begin()) ) diff --git a/src/impl/particles_impl_hskpng_vterm.ipp b/src/impl/particles_impl_hskpng_vterm.ipp index d449e913e..af6839ddf 100644 --- a/src/impl/particles_impl_hskpng_vterm.ipp +++ b/src/impl/particles_impl_hskpng_vterm.ipp @@ -159,7 +159,7 @@ namespace libcloudphxx thrust::make_permutation_iterator(vt_0.begin(), vt0_bin.begin()), thrust::make_permutation_iterator(p.begin(), ijk.begin()), thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), - thrust::make_permutation_iterator(eta.begin(), ijk.begin()) + thrust::make_permutation_iterator(eta.begin_ref(), ijk.begin_ref()) // TODO: use eta on normal grid, currently it is not calculated; same in multiple lines below )), // input - 2nd arg vt.begin(), // condition argument vt.begin(), // output @@ -209,7 +209,7 @@ namespace libcloudphxx thrust::make_permutation_iterator(vt_0.begin(), vt0_bin.begin()), thrust::make_permutation_iterator(p.begin(), ijk.begin()), thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), - thrust::make_permutation_iterator(eta.begin(), ijk.begin()) + thrust::make_permutation_iterator(eta.begin_ref(), ijk.begin_ref()) )), // input - 2nd arg vt.begin(), // output detail::common__vterm__vt__cached(opts_init.terminal_velocity) From 8ed63cc1a1e63897ff13f9f4dd0c2f2e9f0d4e53 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 26 Sep 2022 15:10:19 +0200 Subject: [PATCH 36/38] wip on moms calc on ref grid --- src/impl/particles_impl_hskpng_vterm.ipp | 4 +- src/impl/particles_impl_moms.ipp | 74 ++++++++++++++---------- 2 files changed, 46 insertions(+), 32 deletions(-) diff --git a/src/impl/particles_impl_hskpng_vterm.ipp b/src/impl/particles_impl_hskpng_vterm.ipp index af6839ddf..85438080a 100644 --- a/src/impl/particles_impl_hskpng_vterm.ipp +++ b/src/impl/particles_impl_hskpng_vterm.ipp @@ -157,7 +157,7 @@ namespace libcloudphxx rw2.begin(), rw2.end(), // input - 1st arg zip_it_t(thrust::make_tuple( thrust::make_permutation_iterator(vt_0.begin(), vt0_bin.begin()), - thrust::make_permutation_iterator(p.begin(), ijk.begin()), + thrust::make_permutation_iterator(p.begin_ref(), ijk.begin_ref()), thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin_ref(), ijk.begin_ref()) // TODO: use eta on normal grid, currently it is not calculated; same in multiple lines below )), // input - 2nd arg @@ -207,7 +207,7 @@ namespace libcloudphxx rw2.begin(), rw2.end(), // input - 1st arg zip_it_t(thrust::make_tuple( thrust::make_permutation_iterator(vt_0.begin(), vt0_bin.begin()), - thrust::make_permutation_iterator(p.begin(), ijk.begin()), + thrust::make_permutation_iterator(p.begin_ref(), ijk.begin_ref()), thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin_ref(), ijk.begin_ref()) )), // input - 2nd arg diff --git a/src/impl/particles_impl_moms.ipp b/src/impl/particles_impl_moms.ipp index aad203795..6f73b5149 100644 --- a/src/impl/particles_impl_moms.ipp +++ b/src/impl/particles_impl_moms.ipp @@ -192,7 +192,8 @@ namespace libcloudphxx const typename thrust_device::vector::iterator &vec_bgn, const thrust_size_t npart, const real_t power, - const bool specific + const bool specific, + const bool refined ) { assert(selected_before_counting); @@ -200,6 +201,18 @@ namespace libcloudphxx // same as above thrust_device::vector &n_filtered(tmp_device_real_part); + thrust_device::vector + &_sorted_ijk ( refined ? sorted_ijk.get_ref() : sorted_ijk.get()), + &_sorted_id ( refined ? sorted_id.get_ref() : sorted_id.get()), + &_count_ijk ( refined ? count_ijk.get_ref() : count_ijk.get()); + thrust_device::vector + &_count_mom ( refined ? count_mom.get_ref() : count_mom.get()), + &_rhod ( refined ? rhod.get_ref() : rhod.get()), + &_dv ( refined ? dv.get_ref() : dv.get()); + auto + &_count_n ( refined ? count_n.get_ref() : count_n.get()), + &_n_cell ( refined ? n_cell.get_ref() : n_cell.get()); + typedef thrust::permutation_iterator< typename thrust_device::vector::const_iterator, typename thrust_device::vector::iterator @@ -211,47 +224,47 @@ namespace libcloudphxx typename thrust_device::vector::iterator > it_pair = thrust::reduce_by_key( // input - keys - sorted_ijk.begin(), sorted_ijk.begin()+npart, + _sorted_ijk.begin(), _sorted_ijk.begin()+npart, // input - values thrust::make_transform_iterator( zip_it_t(thrust::make_tuple( - pi_t(n_filtered.begin(), sorted_id.begin()), - pi_t(vec_bgn, sorted_id.begin()) + pi_t(n_filtered.begin(), _sorted_id.begin()), + pi_t(vec_bgn, _sorted_id.begin()) )), detail::moment_counter(power) ), // output - keys - count_ijk.begin(), + _count_ijk.begin(), // output - values - count_mom.begin() + _count_mom.begin() ); - count_n.get() = it_pair.first - count_ijk.begin(); + _count_n = it_pair.first - _count_ijk.begin(); #if !defined(NDEBUG) { - int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n.get(), isnaninf(), 0, thrust::plus()); + int nan_count = thrust::transform_reduce(_count_mom.begin(), _count_mom.begin() + _count_n, isnaninf(), 0, thrust::plus()); if(nan_count>0) { std::cout << nan_count << " nan/inf numbers detected in count_mom after reduce_by_key " << std::endl; } } #endif - assert(count_n.get() <= n_cell.get()); + assert(_count_n <= _n_cell); if(specific) { // dividing by dv thrust::transform( - count_mom.begin(), count_mom.begin() + count_n.get(), // input - first arg + _count_mom.begin(), _count_mom.begin() + _count_n, // input - first arg thrust::make_permutation_iterator( // input - second arg - dv.begin(), - count_ijk.begin() + _dv.begin(), + _count_ijk.begin() ), - count_mom.begin(), // output (in place) + _count_mom.begin(), // output (in place) thrust::divides() ); #if !defined(NDEBUG) { - int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n.get(), isnaninf(), 0, thrust::plus()); + int nan_count = thrust::transform_reduce(_count_mom.begin(), _count_mom.begin() + _count_n, isnaninf(), 0, thrust::plus()); if(nan_count>0) { std::cout << nan_count << " nan/inf numbers detected in count_mom after dividing by dv " << std::endl; @@ -261,17 +274,17 @@ namespace libcloudphxx // dividing by rhod to get specific moments // (for compatibility with blk_1m and blk_2m reporting mixing ratios) thrust::transform( - count_mom.begin(), count_mom.begin() + count_n.get(), // input - first arg - thrust::make_permutation_iterator( // input - second arg - rhod.begin(), - count_ijk.begin() + _count_mom.begin(), _count_mom.begin() + _count_n, // input - first arg + thrust::make_permutation_iterator( // input - second arg + _rhod.begin(), + _count_ijk.begin() ), - count_mom.begin(), // output (in place) + _count_mom.begin(), // output (in place) thrust::divides() ); #if !defined(NDEBUG) { - int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n.get(), isnaninf(), 0, thrust::plus()); + int nan_count = thrust::transform_reduce(_count_mom.begin(), _count_mom.begin() + _count_n, isnaninf(), 0, thrust::plus()); if(nan_count>0) { std::cout << nan_count << " nan/inf numbers detected in count_mom after dividing by rhod " << std::endl; @@ -281,27 +294,27 @@ namespace libcloudphxx } #if !defined(NDEBUG) { - int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n.get(), isnaninf(), 0, thrust::plus()); + int nan_count = thrust::transform_reduce(_count_mom.begin(), _count_mom.begin() + _count_n, isnaninf(), 0, thrust::plus()); if(nan_count>0) { std::cout << nan_count << " nan/inf numbers detected in count_mom " << std::endl; - std::cout << "count_n:" << count_n.get() << std::endl; + std::cout << "count_n:" << _count_n << std::endl; std::cout << "count_mom:" << std::endl; - debug::print(count_mom.begin(), count_mom.begin() + count_n.get()); + debug::print(_count_mom.begin(), _count_mom.begin() + _count_n); std::cout << "count_ijk:" << std::endl; - debug::print(count_ijk.begin(), count_ijk.begin() + count_n.get()); + debug::print(_count_ijk.begin(), _count_ijk.begin() + _count_n); std::cout << "n_filtered:" << std::endl; debug::print(n_filtered); std::cout << "sorted_ijk:" << std::endl; - debug::print(sorted_ijk); + debug::print(_sorted_ijk); std::cout << "sorted_id:" << std::endl; - debug::print(sorted_id); + debug::print(_sorted_id); std::cout << "vec:" << std::endl; debug::print(vec_bgn, vec_bgn + npart); std::cout << "dv:" << std::endl; - debug::print(dv); + debug::print(_dv); std::cout << "rhod:" << std::endl; - debug::print(rhod); + debug::print(_rhod); } } #endif @@ -311,10 +324,11 @@ namespace libcloudphxx void particles_t::impl::moms_calc( const typename thrust_device::vector::iterator &vec_bgn, const real_t power, - const bool specific + const bool specific, + const bool refined ) { - moms_calc(vec_bgn, n_part, power, specific); + moms_calc(vec_bgn, n_part, power, specific, refined); } }; }; From ad76ec15645b1de065a2b7f684659d208a1666a7 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 26 Sep 2022 18:29:17 +0200 Subject: [PATCH 37/38] wip on moms calc on ref grid --- src/impl/particles_impl.ipp | 34 +++++++------- src/impl/particles_impl_dist_analysis.ipp | 2 +- src/impl/particles_impl_hskpng_Tpr.ipp | 54 ++++++++++++++++++----- src/impl/particles_impl_init_grid.ipp | 3 -- src/impl/particles_impl_init_sync.ipp | 2 +- src/impl/particles_impl_moms.ipp | 7 ++- src/impl/particles_impl_sstp.ipp | 4 +- src/particles_init.ipp | 8 ++-- src/particles_step.ipp | 6 +-- 9 files changed, 76 insertions(+), 44 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 2ff9c1b93..2ed820081 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -101,14 +101,11 @@ namespace libcloudphxx // sea level term velocity according to Beard 1977, compute once thrust_device::vector vt_0; - // grid-cell volumes (per grid cell) - thrust_device::vector dv; - // housekeeping data (per particle) thrust_device::vector i, j, k, sorted_id; // Eulerian grid cell indices (always zero for 0D), one sorted_id needed, because its the same in normal and refined grids as we always sort them together ref_part ijk, sorted_ijk; // ijk in the normal and in the refined grid - ref_grid ijk_ref2ijk; // maps refined cell index to the index of the normal cell containing this refined cell +// ref_grid ijk_ref2ijk; // maps refined cell index to the index of the normal cell containing this refined cell // helpers that reference refined arrays if refinement is done and non-refined otherwise // used to conserve memory if n_ref==1 (no refinement) @@ -139,7 +136,8 @@ namespace libcloudphxx rv, // water vapour mixing ratio sstp_pc_tmp_rv, // either rv_old or advection-caused change in water vapour mixing ratio, per-cell substepping sstp_pc_tmp_th, // ditto for theta - sstp_pc_tmp_rh; // ditto for rho + sstp_pc_tmp_rh, // ditto for rho + dv; // cell volume thrust_device::vector sstp_tmp_chem_0, // trace gases @@ -158,13 +156,13 @@ namespace libcloudphxx std::map output_puddle; thrust_device::vector - rhod, // dry air density diss_rate; // turbulent kinetic energy dissipation rate ref_grid T, // temperature [K] p, // pressure [Pa] RH, // relative humisity + rhod, // dry air density eta; // dynamic viscosity thrust_device::vector w_LS; // large-scale subsidence velocity profile @@ -334,7 +332,7 @@ namespace libcloudphxx p(0, n_cell.get_ref()), th(0, n_cell.get_ref()), rv(0, n_cell.get_ref()), - ijk_ref2ijk(0, n_cell.get_ref()), +// ijk_ref2ijk(0, n_cell.get_ref()), count_ijk(n_cell), count_num(n_cell), count_mom(n_cell), @@ -347,6 +345,8 @@ namespace libcloudphxx T (0, n_cell.get_ref()), RH(0, n_cell.get_ref()), eta(0, n_cell.get_ref()), + rhod(n_cell), + dv(n_cell), zero(0), n_part(0), sorted(false), @@ -382,13 +382,13 @@ namespace libcloudphxx adve_scheme(_opts_init.adve_scheme), allow_sstp_cond(_opts_init.sstp_cond > 1 || _opts_init.variable_dt_switch), allow_sstp_chem(_opts_init.sstp_chem > 1 || _opts_init.variable_dt_switch), - sstp_pc_tmp_rv(0, n_cell.get_ref()), - sstp_pc_tmp_th(0, n_cell.get_ref()), - sstp_pc_tmp_rh(0, n_cell.get_ref()), - sstp_tmp_rv(_opts_init.exact_sstp_cond ? sstp_pp_tmp_rv : sstp_pc_tmp_rv.get_ref()), - sstp_tmp_th(_opts_init.exact_sstp_cond ? sstp_pp_tmp_th : sstp_pc_tmp_th.get_ref()), - sstp_tmp_rh(_opts_init.exact_sstp_cond ? sstp_pp_tmp_rh : sstp_pc_tmp_rh.get_ref()), - sstp_tmp_p(sstp_pp_tmp_p), // not used in per-cell substepping + sstp_pc_tmp_rv(0, _opts_init.exact_sstp_cond ? 0 : n_cell.get_ref()), + sstp_pc_tmp_th(0, _opts_init.exact_sstp_cond ? 0 : n_cell.get_ref()), + sstp_pc_tmp_rh(0, _opts_init.exact_sstp_cond ? 0 : n_cell.get_ref()), + sstp_tmp_rv(_opts_init.exact_sstp_cond ? sstp_pp_tmp_rv : sstp_pc_tmp_rv.get_ref()), + sstp_tmp_th(_opts_init.exact_sstp_cond ? sstp_pp_tmp_th : sstp_pc_tmp_th.get_ref()), + sstp_tmp_rh(_opts_init.exact_sstp_cond ? sstp_pp_tmp_rh : sstp_pc_tmp_rh.get_ref()), + sstp_tmp_p(sstp_pp_tmp_p), // not used in per-cell substepping pure_const_multi (((_opts_init.sd_conc) == 0) && (_opts_init.sd_const_multi > 0 || _opts_init.dry_sizes.size() > 0)), // coal prob can be greater than one only in sd_conc simulations ijk(opts_init.n_ref), sorted_ijk(opts_init.n_ref) @@ -580,12 +580,14 @@ namespace libcloudphxx const typename thrust_device::vector::iterator &vec_bgn, const thrust_size_t npart, const real_t power, - const bool specific = true + const bool specific = true, + const bool refined = false ); void moms_calc( const typename thrust_device::vector::iterator &vec_bgn, const real_t power, - const bool specific = true + const bool specific = true, + const bool refined = false ); void mass_dens_estim( diff --git a/src/impl/particles_impl_dist_analysis.ipp b/src/impl/particles_impl_dist_analysis.ipp index 74cafaafd..fa2fd9bff 100644 --- a/src/impl/particles_impl_dist_analysis.ipp +++ b/src/impl/particles_impl_dist_analysis.ipp @@ -35,7 +35,7 @@ namespace libcloudphxx / sd_conc * dt * (n_dims == 0 - ? dv[0] + ? dv.get()[0] : (opts_init.dx * opts_init.dy * opts_init.dz) ); diff --git a/src/impl/particles_impl_hskpng_Tpr.ipp b/src/impl/particles_impl_hskpng_Tpr.ipp index 3e8818dff..97c9692d7 100644 --- a/src/impl/particles_impl_hskpng_Tpr.ipp +++ b/src/impl/particles_impl_hskpng_Tpr.ipp @@ -178,7 +178,7 @@ namespace libcloudphxx // T = common::theta_dry::T(th, rhod); thrust::transform( th.begin_ref(), th.end_ref(), // input - first arg - thrust::make_permutation_iterator(rhod.begin(), ijk_ref2ijk.begin_ref()), // input - second arg + rhod.begin_ref(), // input - second arg T.begin_ref(), // output detail::common__theta_dry__T_rhod() ); @@ -199,7 +199,7 @@ namespace libcloudphxx { // p = common::theta_dry::p(rhod, r, T); auto it = thrust::make_zip_iterator(thrust::make_tuple( - thrust::make_permutation_iterator(rhod.begin(), ijk_ref2ijk.begin_ref()), + rhod.begin_ref(), rv.begin_ref(), T.begin_ref() )); @@ -229,19 +229,53 @@ namespace libcloudphxx // dynamic viscosity { - // on refined grid - thrust::transform( - T.begin_ref(), T.end_ref(), // 1st arg - eta.begin_ref(), // output - detail::common__vterm__visc() - ); // on normal grid: how if we dont know T, because we dont know th? - // TODO: require input of normal th? - // or use averaging from refined to normal? + // OTHER OPTION: use averaging from refined to normal? // cumbersome: // 0. make copies of eta and ijk_ref2ijk // 1. sort copies of eta by ijk_ref2ijk // 2. reduce_by_key and add the result to normal eta + // + // calc based on th on normal grid: + // T on normal grid: + thrust_device::vector &Tn(tmp_device_real_cell.get()); + + // same code as for T on refined grid, TODO: DRY + if(!const_p) // variable pressure + { + // T = common::theta_dry::T(th, rhod); + thrust::transform( + th.begin(), th.end(), // input - first arg + rhod.begin(), // input - second arg + Tn.begin(), // output + detail::common__theta_dry__T_rhod() + ); + } + else // external pressure profile + { + // T = th * exner(p_tot) + thrust::transform( + th.begin(), th.end(), // input - first arg + thrust::make_tuple(rv.begin(), p.begin()), // input - second and third args + Tn.begin(), // output + detail::common__theta_dry__T_p() + ); + } + + thrust::transform( + Tn.begin(), Tn.end(), // 1st arg + eta.begin(), // output + detail::common__vterm__visc() + ); + + + // on refined grid + if(opts_init.n_ref > 1) + thrust::transform( + T.begin_ref(), T.end_ref(), // 1st arg + eta.begin_ref(), // output + detail::common__vterm__visc() + ); } diff --git a/src/impl/particles_impl_init_grid.ipp b/src/impl/particles_impl_init_grid.ipp index 2a9759999..23deb3f99 100644 --- a/src/impl/particles_impl_init_grid.ipp +++ b/src/impl/particles_impl_init_grid.ipp @@ -58,9 +58,6 @@ namespace libcloudphxx { namespace arg = thrust::placeholders; - // filling in sample volume data - dv.resize(n_cell.get()); - int n_cell_halo(2 * halo_x); // number of halo cells if (n_dims > 0) diff --git a/src/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index 616ac01e0..f1765319d 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/src/impl/particles_impl_init_sync.ipp @@ -13,7 +13,7 @@ namespace libcloudphxx void particles_t::impl::init_sync() { // memory allocation for scalar fields - rhod.resize(n_cell.get()); +// rhod.resize(n_cell.get()); // p.resize(n_cell); // th.resize(n_cell); // rv.resize(n_cell); diff --git a/src/impl/particles_impl_moms.ipp b/src/impl/particles_impl_moms.ipp index 6f73b5149..b145235cc 100644 --- a/src/impl/particles_impl_moms.ipp +++ b/src/impl/particles_impl_moms.ipp @@ -203,7 +203,6 @@ namespace libcloudphxx thrust_device::vector &_sorted_ijk ( refined ? sorted_ijk.get_ref() : sorted_ijk.get()), - &_sorted_id ( refined ? sorted_id.get_ref() : sorted_id.get()), &_count_ijk ( refined ? count_ijk.get_ref() : count_ijk.get()); thrust_device::vector &_count_mom ( refined ? count_mom.get_ref() : count_mom.get()), @@ -228,8 +227,8 @@ namespace libcloudphxx // input - values thrust::make_transform_iterator( zip_it_t(thrust::make_tuple( - pi_t(n_filtered.begin(), _sorted_id.begin()), - pi_t(vec_bgn, _sorted_id.begin()) + pi_t(n_filtered.begin(), sorted_id.begin()), + pi_t(vec_bgn, sorted_id.begin()) )), detail::moment_counter(power) ), @@ -308,7 +307,7 @@ namespace libcloudphxx std::cout << "sorted_ijk:" << std::endl; debug::print(_sorted_ijk); std::cout << "sorted_id:" << std::endl; - debug::print(_sorted_id); + debug::print(sorted_id); std::cout << "vec:" << std::endl; debug::print(vec_bgn, vec_bgn + npart); std::cout << "dv:" << std::endl; diff --git a/src/impl/particles_impl_sstp.ipp b/src/impl/particles_impl_sstp.ipp index 5970229f4..6b7a6d5f3 100644 --- a/src/impl/particles_impl_sstp.ipp +++ b/src/impl/particles_impl_sstp.ipp @@ -48,8 +48,8 @@ namespace libcloudphxx const int n = 4; thrust_device::vector - *fr[n] = { rv.ptr_ref(), th.ptr_ref(), &rhod, p.ptr_ref() }, - *to[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh, &sstp_tmp_p }; + *fr[n] = { rv.ptr_ref(), th.ptr_ref(), rhod.ptr_ref(), p.ptr_ref() }, + *to[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh, &sstp_tmp_p }; for (int ix = 0; ix < ( (const_p) ? n : n-1); ++ix) // TODO: var_rho { thrust::copy( diff --git a/src/particles_init.ipp b/src/particles_init.ipp index a9c8d5460..8c1843c67 100644 --- a/src/particles_init.ipp +++ b/src/particles_init.ipp @@ -41,8 +41,8 @@ namespace libcloudphxx pimpl->init_sync(); // also, init of ambient_chem vectors pimpl->init_e2l(th, pimpl->th.ptr_ref(), true); pimpl->init_e2l(rv, pimpl->rv.ptr_ref(), true); - pimpl->init_e2l(rhod, &pimpl->rhod); - //pimpl->init_e2l(rhod, pimpl->rhod.ptr_ref(), true); +// pimpl->init_e2l(rhod, &pimpl->rhod); + pimpl->init_e2l(rhod, pimpl->rhod.ptr_ref(), true); if(pimpl->const_p) pimpl->init_e2l(p, pimpl->p.ptr_ref(), true); @@ -60,8 +60,8 @@ namespace libcloudphxx // feeding in Eulerian fields pimpl->sync(th, pimpl->th.get_ref()); pimpl->sync(rv, pimpl->rv.get_ref()); - pimpl->sync(rhod, pimpl->rhod); - pimpl->sync(p, pimpl->p.get_ref()); + pimpl->sync(rhod, pimpl->rhod.get_ref()); + pimpl->sync(p, pimpl->p.get_ref()); pimpl->sync(courant_x, pimpl->courant_x); pimpl->sync(courant_y, pimpl->courant_y); diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 055ad7a1b..edf04fd0e 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -96,7 +96,7 @@ namespace libcloudphxx pimpl->sync(th, pimpl->th.get_ref()); pimpl->sync(rv, pimpl->rv.get_ref()); pimpl->sync(diss_rate, pimpl->diss_rate); - pimpl->sync(rhod, pimpl->rhod); + pimpl->sync(rhod, pimpl->rhod.get_ref()); pimpl->sync(courant_x, pimpl->courant_x); pimpl->sync(courant_y, pimpl->courant_y); pimpl->sync(courant_z, pimpl->courant_z); @@ -110,13 +110,13 @@ namespace libcloudphxx 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"); + nancheck(pimpl->rhod.get_ref(), " rhod after sync-in"); if(pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_switch || pimpl->opts_init.turb_coal_switch) {nancheck(pimpl->diss_rate, " diss_rate after sync-in");} assert(*thrust::min_element(pimpl->rv.begin_ref(), pimpl->rv.end_ref()) >= 0); assert(*thrust::min_element(pimpl->th.begin_ref(), pimpl->th.end_ref()) >= 0); - assert(*thrust::min_element(pimpl->rhod.begin(), pimpl->rhod.end()) >= 0); + assert(*thrust::min_element(pimpl->rhod.begin_ref(), pimpl->rhod.end_ref()) >= 0); if(pimpl->opts_init.turb_adve_switch || pimpl->opts_init.turb_cond_switch || pimpl->opts_init.turb_coal_switch) {assert(*thrust::min_element(pimpl->diss_rate.begin(), pimpl->diss_rate.end()) >= 0);} From 0cf3e063dad59fd044f038aa85b093c9b41488bf Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Tue, 27 Sep 2022 14:54:03 +0200 Subject: [PATCH 38/38] ref wip --- src/impl/particles_impl.ipp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 2ed820081..188ce887b 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -329,9 +329,12 @@ namespace libcloudphxx nz_ref(n_dims >= 2 ? (opts_init.nz - 1) * opts_init.n_ref + 1 : 0), n_cell(m1(_opts_init.nx) * m1(_opts_init.ny) * m1(_opts_init.nz), m1(nx_ref) * m1(ny_ref) * m1(nz_ref)), // NOTE: needs to be equal to n_cell for n_ref == 1 - p(0, n_cell.get_ref()), - th(0, n_cell.get_ref()), - rv(0, n_cell.get_ref()), + // p(0, n_cell.get_ref()), + // th(0, n_cell.get_ref()), + // rv(0, n_cell.get_ref()), + p (n_cell), + th(n_cell), + rv(n_cell), // ijk_ref2ijk(0, n_cell.get_ref()), count_ijk(n_cell), count_num(n_cell), @@ -344,7 +347,8 @@ namespace libcloudphxx tmp_host_size_cell(n_cell), T (0, n_cell.get_ref()), RH(0, n_cell.get_ref()), - eta(0, n_cell.get_ref()), + //eta(0, n_cell.get_ref()), + eta(n_cell), rhod(n_cell), dv(n_cell), zero(0),