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/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index c012d8823..65e1a7991 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -46,6 +46,14 @@ 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; + // 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 + // 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; @@ -65,7 +73,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 @@ -236,7 +244,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/detail/debug.hpp b/src/detail/debug.hpp new file mode 100644 index 000000000..567cca741 --- /dev/null +++ b/src/detail/debug.hpp @@ -0,0 +1,44 @@ +#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.cbegin(), v.cend()); + std::cerr << "refined grid values:" << std::endl; + 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 new file mode 100644 index 000000000..b71a8eb6a --- /dev/null +++ b/src/detail/grid_refinement_helpers.hpp @@ -0,0 +1,230 @@ +// helper classes to store arrays associated with grid refinement +#include + +#pragma once + +namespace libcloudphxx +{ + namespace lgrngn + { + template + class ref_common + { + public: + const bool ref_flag; // true if refinement is actually done + + protected: + // ctor + ref_common(const bool ref_flag): + ref_flag(ref_flag) + {} + + ref_common() = delete; + }; + + + // for single values of type T + template + class ref_val : ref_common + { + private: + using parent_t = ref_common; + T val, val_ref; + + public: + ref_val(T val, T val_ref): + val(val), + val_ref(val_ref), + parent_t(true) + {} + ref_val() = delete; + + T& get() + { + 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 + class ref_arr_common : public ref_common + { + private: + using parent_t = ref_common; + using vec_t = typename std::conditional, thrust_device::vector>::type; + + protected: + using parent_t::parent_t; + vec_t arr, arr_ref; // actual data, arr_ref initialized only if refinement is actually done + + void resize(thrust_size_t size, thrust_size_t size_ref) + { + arr.resize(size); + 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(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 begin() + { + 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; + } + + 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 that need values + // in the normal grid and in the refined grid (or in only one of these) + template // value type and container type + class ref_grid : public ref_arr_common + { + using parent_t = ref_arr_common; + + public: + // ctor + 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); + } + + 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 that need + // two independent values: for the normal and for the refined grid (e.g. ijk) + template // value type and container type + class ref_part : public ref_arr_common + { + using parent_t = ref_arr_common; + + public: + // ctor + ref_part(const int &n_ref): + parent_t(n_ref > 1) + {assert(n_ref > 0);} + + ref_part() = delete; + + void resize(int n) + { + parent_t::resize(n, n); + } + void reserve(int n) + { + parent_t::reserve(n, n); + } + }; + }; +}; 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.ipp b/src/impl/particles_impl.ipp index b07d3332c..188ce887b 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 @@ -45,7 +44,8 @@ namespace libcloudphxx // member fields opts_init_t opts_init; // a copy const int n_dims; - const thrust_size_t n_cell; + const int nx_ref, ny_ref, nz_ref; // number of refined cells in each direction + 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,10 +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_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 + 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 @@ -101,38 +101,51 @@ 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, ijk, // Eulerian grid cell indices (always zero for 0D) - sorted_id, sorted_ijk; + 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 + // 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 lft, rgt, abv, blw, fre, hnd; // TODO: could be reused after advection! - // moment-counting stuff - thrust_device::vector + // moment-counting stuff on normal and refined grid + ref_grid count_ijk; // key-value pair for sorting particles by cell index - thrust_device::vector + ref_grid count_num; // number of particles in a given grid cell - thrust_device::vector + ref_grid count_mom; // statistical moment // TODO (perhaps tmp_device_real_cell could be referenced?) - thrust_size_t count_n; + ref_val count_n; // Eulerian-Lagrangian interface vars - thrust_device::vector - rhod, // dry air density + ref_grid th, // potential temperature (dry) rv, // 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 + 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 + dv; // cell volume + + 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; @@ -141,13 +154,16 @@ 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 - T, // temperature [K] - p, // pressure [Pa] - RH, // relative humisity - eta,// dynamic viscosity 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 thrust_device::vector SGS_mix_len; // SGS mixing length profile @@ -170,6 +186,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; @@ -205,10 +229,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, @@ -216,10 +237,15 @@ namespace libcloudphxx tmp_device_real_part3, tmp_device_real_part4, tmp_device_real_part5, + &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, - &u01; // uniform random numbers between 0 and 1 // TODO: use the tmp array as rand argument? + 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 @@ -298,11 +324,33 @@ 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 + // 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), + 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), + 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(n_cell), + rhod(n_cell), + dv(n_cell), zero(0), n_part(0), sorted(false), @@ -312,7 +360,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), @@ -338,7 +386,16 @@ 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), - 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 + 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) { // set 0 dev_count to mark that its not a multi_CUDA spawn @@ -379,7 +436,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); } @@ -467,7 +524,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(); @@ -480,6 +537,7 @@ namespace libcloudphxx void init_vterm(); void fill_outbuf(); + void fill_outbuf_ref(); void mpi_exchange(); // rename hskpng_ -> step_? @@ -487,9 +545,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(); @@ -525,12 +584,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( @@ -560,8 +621,9 @@ 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 &); void update_incloud_time(const real_t &dt); diff --git a/src/impl/particles_impl_ante_adding_SD.ipp b/src/impl/particles_impl_ante_adding_SD.ipp index b77502a11..9e5a603dd 100644 --- a/src/impl/particles_impl_ante_adding_SD.ipp +++ b/src/impl/particles_impl_ante_adding_SD.ipp @@ -14,16 +14,16 @@ 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.)); // drv = - tot_vol_bfr thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output + count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg + thrust::make_permutation_iterator(drv.begin_ref(), count_ijk.begin_ref()), // output thrust::negate() ); diff --git a/src/impl/particles_impl_chem_henry.ipp b/src/impl/particles_impl_chem_henry.ipp index fabaa8f44..b4f5e8d79 100644 --- a/src/impl/particles_impl_chem_henry.ipp +++ b/src/impl/particles_impl_chem_henry.ipp @@ -307,8 +307,8 @@ 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_new(tmp_device_real_cell1); + 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) { @@ -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); + 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 205ccacb2..de5b7b3c7 100644 --- a/src/impl/particles_impl_coal.ipp +++ b/src/impl/particles_impl_coal.ipp @@ -265,15 +265,15 @@ 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 - &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 @@ -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) 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) 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 @@ -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_cond.ipp b/src/impl/particles_impl_cond.ipp index c3fdfa618..d77017b2f 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -18,37 +18,39 @@ 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 + auto &lambda_D(tmp_device_real_cell1); + auto &lambda_K(tmp_device_real_cell2); // --- 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(); 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) 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, // input - 1st arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output + count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg + thrust::make_permutation_iterator(drv.begin_ref(), count_ijk.begin_ref()), // output thrust::negate() ); auto hlpr_zip_iter = thrust::make_zip_iterator(thrust::make_tuple( thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), - thrust::make_permutation_iterator(rv.begin(), ijk.begin()), - thrust::make_permutation_iterator(T.begin(), ijk.begin()), - thrust::make_permutation_iterator(eta.begin(), ijk.begin()), + 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 @@ -58,7 +60,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.begin_ref(), ijk.begin_ref()), RH_plus_ssp.begin(), arg::_1 + arg::_2 ); @@ -68,7 +70,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() ) ), @@ -82,8 +84,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.begin(), ijk.begin()) + thrust::make_permutation_iterator(p.begin_ref(), ijk.begin_ref()), + thrust::make_permutation_iterator(RH.begin_ref(), ijk.begin_ref()) ) ), rw2.begin(), // output @@ -93,13 +95,13 @@ 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 - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // input - 2nd arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output + count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg + 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_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_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_fill_outbuf.ipp b/src/impl/particles_impl_fill_outbuf.ipp index 2855fe550..ab78cfd4f 100644 --- a/src/impl/particles_impl_fill_outbuf.ipp +++ b/src/impl/particles_impl_fill_outbuf.ipp @@ -22,18 +22,46 @@ 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( - count_mom.begin(), // input - begin - count_mom.begin() + count_n, // 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() + { + 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_Tpr.ipp b/src/impl/particles_impl_hskpng_Tpr.ipp index 38433bcad..97c9692d7 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.begin(), th.end(), // input - first arg - rhod.begin(), // input - second arg - T.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,50 +187,98 @@ 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 - T.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() ); } { - 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); + auto it = thrust::make_zip_iterator(thrust::make_tuple( + rhod.begin_ref(), + rv.begin_ref(), + T.begin_ref() + )); 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 - p.begin(), // output + it, // input - begin + it + n_cell.get_ref(), // input - end + p.begin_ref(), // output detail::common__theta_dry__p() ); } + 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(), 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_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) ); } // dynamic viscosity { + // on normal grid: how if we dont know T, because we dont know th? + // 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( - T.begin(), T.end(), // 1st arg + 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() + ); } + // adjusting dv if using a parcel set-up (1kg of dry air) if (n_dims == 0) { diff --git a/src/impl/particles_impl_hskpng_count.ipp b/src/impl/particles_impl_hskpng_count.ipp index 2619094c1..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); + count_n.get() = it_pair.first - count_ijk.begin(); + assert(count_n.get() <= n_cell.get()); } }; }; diff --git a/src/impl/particles_impl_hskpng_ijk.ipp b/src/impl/particles_impl_hskpng_ijk.ipp index 2320462f2..a598c1125 100644 --- a/src/impl/particles_impl_hskpng_ijk.ipp +++ b/src/impl/particles_impl_hskpng_ijk.ipp @@ -15,36 +15,37 @@ 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); } }; }; // 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_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; @@ -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(), - 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 @@ -122,23 +123,61 @@ 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 + const real_t &vd, + const real_t pos_shift = 0 ) { thrust::transform( - vx.begin(), vx.end(), // input - vi.begin(), // 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? + vx.begin() + begin_shift, vx.end(), // input + vi.begin() + begin_shift, // output + 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? ); } - } helper; + + cell_idx_hlpr_t(const thrust_size_t begin_shift = 0): + begin_shift(begin_shift) + {} + }; + }; + + // 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) + { + if(opts_init.n_ref == 1) return; + + cell_idx_hlpr_t helper(begin_shift); + + 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); + + // 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_mfp.ipp b/src/impl/particles_impl_hskpng_mfp.ipp index 1cf02535b..e875436a1 100644 --- a/src/impl/particles_impl_hskpng_mfp.ipp +++ b/src/impl/particles_impl_hskpng_mfp.ipp @@ -41,11 +41,12 @@ 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 + auto &lambda_D(tmp_device_real_cell1); + auto &lambda_K(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.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_hskpng_remove.ipp b/src/impl/particles_impl_hskpng_remove.ipp index ae31eee6b..1d4e343e0 100644 --- a/src/impl/particles_impl_hskpng_remove.ipp +++ b/src/impl/particles_impl_hskpng_remove.ipp @@ -21,12 +21,13 @@ 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); + std::set*> n_t_vctrs_ref; + n_t_vctrs_ref.insert(&ijk); + namespace arg = thrust::placeholders; // remove chemical stuff @@ -75,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 d9f494bb7..7311688ed 100644 --- a/src/impl/particles_impl_hskpng_resize.ipp +++ b/src/impl/particles_impl_hskpng_resize.ipp @@ -15,8 +15,8 @@ namespace libcloudphxx } } { - thrust_device::vector *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 4f4abcd55..0a0e15f15 100644 --- a/src/impl/particles_impl_hskpng_sort.ipp +++ b/src/impl/particles_impl_hskpng_sort.ipp @@ -19,11 +19,16 @@ 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_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_hskpng_turb_ss.ipp b/src/impl/particles_impl_hskpng_turb_ss.ipp index f14d8ea74..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); @@ -55,7 +56,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 +68,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_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_hskpng_vterm.ipp b/src/impl/particles_impl_hskpng_vterm.ipp index 765d7768c..85438080a 100644 --- a/src/impl/particles_impl_hskpng_vterm.ipp +++ b/src/impl/particles_impl_hskpng_vterm.ipp @@ -157,9 +157,9 @@ 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(), 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 @@ -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.begin(), ijk.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(T.begin_ref(), ijk.begin_ref()), + 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 vt.begin(), // condition argument vt.begin(), // output @@ -207,9 +207,9 @@ 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(), 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) @@ -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.begin(), ijk.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(T.begin_ref(), ijk.begin_ref()), + 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 vt.begin(), // output detail::common__vterm__vt(opts_init.terminal_velocity) 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_count_num.ipp b/src/impl/particles_impl_init_count_num.ipp index ad3dcb414..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); @@ -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_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_grid.ipp b/src/impl/particles_impl_init_grid.ipp index bfdcf9319..23deb3f99 100644 --- a/src/impl/particles_impl_init_grid.ipp +++ b/src/impl/particles_impl_init_grid.ipp @@ -58,36 +58,33 @@ namespace libcloudphxx { namespace arg = thrust::placeholders; - // filling in sample volume data - dv.resize(n_cell); - int n_cell_halo(2 * halo_x); // number of halo cells if (n_dims > 0) { // 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 +93,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 +105,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 +117,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 +127,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 +139,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 2b92d6601..8c515316d 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -13,28 +13,31 @@ 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); + // 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 = 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); - 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); +// 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); - sstp_tmp_th.resize(n_cell); - sstp_tmp_rh.resize(n_cell); + 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()); } + + //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 ee2282a33..c80c97ab7 100644 --- a/src/impl/particles_impl_init_ijk.ipp +++ b/src/impl/particles_impl_init_ijk.ipp @@ -44,9 +44,9 @@ 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[n_part_old])) + detail::arbitrary_sequence(&(*(ijk.begin() + n_part_old))) ); } }; diff --git a/src/impl/particles_impl_init_n.ipp b/src/impl/particles_impl_init_n.ipp index 50c4880a7..2aa6f18e8 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 @@ -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_init_sanity_check.ipp b/src/impl/particles_impl_init_sanity_check.ipp index 5d309dafa..eba510ff2 100644 --- a/src/impl/particles_impl_init_sanity_check.ipp +++ b/src/impl/particles_impl_init_sanity_check.ipp @@ -141,6 +141,37 @@ 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"); + // 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/impl/particles_impl_init_sync.ipp b/src/impl/particles_impl_init_sync.ipp index 616792eab..f1765319d 100644 --- a/src/impl/particles_impl_init_sync.ipp +++ b/src/impl/particles_impl_init_sync.ipp @@ -13,15 +13,15 @@ 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.get()); +// 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); + 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_init_wet.ipp b/src/impl/particles_impl_init_wet.ipp index 77646cd47..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.begin(), ijk.begin() + n_part_old), - pi_t(T.begin(), ijk.begin() + 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_mass_dens.ipp b/src/impl/particles_impl_mass_dens.ipp index 7975ce2ac..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(); @@ -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); + 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..b145235cc 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,17 @@ 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()), + &_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,7 +223,7 @@ 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( @@ -221,37 +233,37 @@ namespace libcloudphxx detail::moment_counter(power) ), // output - keys - count_ijk.begin(), + _count_ijk.begin(), // output - values - count_mom.begin() + _count_mom.begin() ); - count_n = 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, 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 <= n_cell); + assert(_count_n <= _n_cell); 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, // 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, 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 +273,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, // 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, 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 +293,27 @@ 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, 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 << 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); 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); 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); 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 +323,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); } }; }; diff --git a/src/impl/particles_impl_post_adding_SD.ipp b/src/impl/particles_impl_post_adding_SD.ipp index 9bebe5352..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; @@ -24,9 +24,9 @@ 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 - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // 2nd arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output + count_mom.begin(), count_mom.begin() + count_n.get(), // input - 1st arg + 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_rcyc.ipp b/src/impl/particles_impl_rcyc.ipp index 053cf559a..b2d9e23a3 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( 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_src_dry_distros_matching.ipp b/src/impl/particles_impl_src_dry_distros_matching.ipp index 61faa20bb..c51e23cb2 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 { @@ -255,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_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/impl/particles_impl_sstp.ipp b/src/impl/particles_impl_sstp.ipp index 28ce6191e..6b7a6d5f3 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 @@ -45,8 +48,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.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( @@ -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) { @@ -108,16 +111,18 @@ namespace libcloudphxx void particles_t::impl::sstp_step_exact( const int &step ) - { + { if (sstp_cond == 1) return; namespace arg = thrust::placeholders; 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 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_update_incloud_time.ipp b/src/impl/particles_impl_update_incloud_time.ipp index 75dbec831..b85937a46 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.begin_ref(), + ijk.begin_ref() ) )), // 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..2caad0baa 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() ), @@ -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,36 +97,36 @@ 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 - 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, typename thrust_device::vector::iterator > it_pair = thrust::reduce_by_key( - sorted_ijk.begin(), sorted_ijk.end(), + sorted_ijk.begin_ref(), sorted_ijk.end_ref(), thrust::make_permutation_iterator(pdstate.begin(), sorted_id.begin()), - count_ijk.begin(), - count_mom.begin() + count_ijk.begin_ref(), + count_mom.begin_ref() ); - count_n = 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, // input - 1st arg - thrust::make_permutation_iterator(dstate.begin(), count_ijk.begin()), // 2nd arg - thrust::make_permutation_iterator(dstate.begin(), count_ijk.begin()), // output + 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() ); // 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() ); @@ -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/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( 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)"); diff --git a/src/particles.tpp b/src/particles.tpp index 477f7a295..2e048dc71 100644 --- a/src/particles.tpp +++ b/src/particles.tpp @@ -28,6 +28,8 @@ #include "detail/kernel_interpolation.hpp" #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" diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index a948143a3..35ad31216 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -113,14 +113,14 @@ 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 - pimpl->count_n = pimpl->n_cell; - thrust::sequence(pimpl->count_ijk.begin(), pimpl->count_ijk.end()); + pimpl->count_n.get() = pimpl->n_cell.get(); + thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } // records temperature @@ -130,14 +130,14 @@ namespace libcloudphxx pimpl->hskpng_Tpr(); thrust::copy( - pimpl->T.begin(), - pimpl->T.end(), - pimpl->count_mom.begin() + pimpl->T.begin_ref(), + pimpl->T.end_ref(), + pimpl->count_mom.begin_ref() ); // T defined in all cells - pimpl->count_n = pimpl->n_cell; - thrust::sequence(pimpl->count_ijk.begin(), pimpl->count_ijk.end()); + pimpl->count_n.get_ref() = pimpl->n_cell.get_ref(); + thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } // records relative humidity @@ -147,14 +147,14 @@ namespace libcloudphxx pimpl->hskpng_Tpr(); thrust::copy( - pimpl->RH.begin(), - pimpl->RH.end(), - pimpl->count_mom.begin() + pimpl->RH.begin_ref(), + pimpl->RH.end_ref(), + pimpl->count_mom.begin_ref() ); // RH defined in all cells - pimpl->count_n = pimpl->n_cell; - thrust::sequence(pimpl->count_ijk.begin(), pimpl->count_ijk.end()); + pimpl->count_n.get_ref() = pimpl->n_cell.get_ref(); + thrust::sequence(pimpl->count_ijk.begin_ref(), pimpl->count_ijk.end_ref()); } // records super-droplet concentration per grid cell @@ -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 @@ -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.begin_ref(), + pimpl->ijk.begin_ref() ), thrust::make_permutation_iterator( - pimpl->RH.begin(), - pimpl->ijk.begin() + pimpl->RH.begin_ref(), + pimpl->ijk.begin_ref() ) )), // 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.begin_ref(), + pimpl->ijk.begin_ref() ) )), // input - 2nd arg rc2.begin(), // output @@ -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.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); + 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 diff --git a/src/particles_init.ipp b/src/particles_init.ipp index 5fcf9554c..8c1843c67 100644 --- a/src/particles_init.ipp +++ b/src/particles_init.ipp @@ -25,47 +25,43 @@ 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); // 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(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; #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) 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); diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 0ffee2b0e..edf04fd0e 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) @@ -98,10 +93,10 @@ 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(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); @@ -109,19 +104,19 @@ 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"); 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(), pimpl->rv.end()) >= 0); - assert(*thrust::min_element(pimpl->th.begin(), pimpl->th.end()) >= 0); - assert(*thrust::min_element(pimpl->rhod.begin(), pimpl->rhod.end()) >= 0); + assert(*thrust::min_element(pimpl->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_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);} @@ -181,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) @@ -211,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(); @@ -264,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)