diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index 504e5f8b..5d1e6ae4 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -53,7 +53,7 @@ namespace libcloudphxx // Lagrangian domain extents real_t x0, y0, z0, x1, y1, z1; - // no. of super-droplets per cell + // no. of super-droplets per cell (or in the entire domain if domain_sd_init==true) unsigned long long sd_conc; // should more SDs be added to better represent large tail of the distribution @@ -128,9 +128,9 @@ namespace libcloudphxx real_t rd_min, rd_max; // min/max dry radius of droplets [m] - bool no_ccn_at_init; // if true, no ccn / SD are put at the start of the simulation - - bool open_side_walls, // if true, side walls are "open", i.e. SD are removed at contact. Periodic otherwise. + bool no_ccn_at_init, // if true, no ccn / SD are put at the start of the simulation + domain_sd_init, // if true, SDs are initialized in the entire domain as if it was a single cell (experimental, could mess with SD sources and other things...) + open_side_walls, // if true, side walls are "open", i.e. SD are removed at contact. Periodic otherwise. periodic_topbot_walls; // if true, top and bot walls are periodic. Open otherwise @@ -234,6 +234,7 @@ namespace libcloudphxx rd_max(-1), diag_incloud_time(false), no_ccn_at_init(false), + domain_sd_init(false), open_side_walls(false), periodic_topbot_walls(false), variable_dt_switch(false), diff --git a/src/impl/particles_impl_dist_analysis.ipp b/src/impl/particles_impl_dist_analysis.ipp index ac63c974..6e03e4e1 100644 --- a/src/impl/particles_impl_dist_analysis.ipp +++ b/src/impl/particles_impl_dist_analysis.ipp @@ -28,8 +28,9 @@ namespace libcloudphxx / sd_conc * dt * (n_dims == 0 - ? dv[0] - : (opts_init.dx * opts_init.dy * opts_init.dz) + ? dv[0] // 0D (box/parcel) + : (opts_init.dx * opts_init.dy * opts_init.dz) // 1D 2D 3D + * (opts_init.domain_sd_init ? opts_init.nx * opts_init.ny * opts_init.nz : 1) // if initializing in the entire domain, not per-cell ); log_rd_min = log(rd_min); @@ -52,8 +53,9 @@ namespace libcloudphxx / sd_conc * dt * (n_dims == 0 - ? dv[0] - : (opts_init.dx * opts_init.dy * opts_init.dz) + ? dv[0] // 0D (box/parcel) + : (opts_init.dx * opts_init.dy * opts_init.dz) // 1D 2D 3D + * (opts_init.domain_sd_init ? opts_init.nx * opts_init.ny * opts_init.nz : 1) // if initializing in the entire domain, not per-cell ); log_rd_min = log(rd_min); diff --git a/src/impl/particles_impl_init_SD_with_distros.ipp b/src/impl/particles_impl_init_SD_with_distros.ipp index 0a20ba6c..0e32901e 100644 --- a/src/impl/particles_impl_init_SD_with_distros.ipp +++ b/src/impl/particles_impl_init_SD_with_distros.ipp @@ -52,6 +52,17 @@ 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(); + + // get correct ijk if domain-initializing + if(opts_init.domain_sd_init) + hskpng_ijk(); + // init kappa init_kappa(kappa); @@ -77,13 +88,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_count_num.ipp b/src/impl/particles_impl_init_count_num.ipp index 73fe5c35..8a6a7ef4 100644 --- a/src/impl/particles_impl_init_count_num.ipp +++ b/src/impl/particles_impl_init_count_num.ipp @@ -31,7 +31,13 @@ namespace libcloudphxx template void particles_t::impl::init_count_num_sd_conc(const real_t &ratio) { - thrust::fill(count_num.begin(), count_num.end(), ratio * opts_init.sd_conc); + if(!opts_init.domain_sd_init) // usual per-cell init + thrust::fill(count_num.begin(), count_num.end(), ratio * opts_init.sd_conc); + else // if initializing in the entire domain, pretend that the first cell is the entire domain + { + thrust::fill(count_num.begin(), count_num.end(), 0); + *(count_num.begin()) = ratio * opts_init.sd_conc; + } } // calculate number of droplets in a cell from concentration [1/m^3], taking into account cell volume and air density @@ -44,22 +50,49 @@ namespace libcloudphxx namespace arg = thrust::placeholders; using common::earth::rho_stp; - // cell volume - thrust::transform( - dv.begin(), dv.end(), - arr.begin(), - arr.begin(), - arg::_2 * arg::_1 - ); - - // correct for density with respect to STP - if(!opts_init.aerosol_independent_of_rhod) + if(!opts_init.domain_sd_init || n_dims==0) // usual per-cell init + { + // cell volume thrust::transform( - rhod.begin(), rhod.end(), + dv.begin(), dv.end(), arr.begin(), arr.begin(), - arg::_1 / real_t(rho_stp() / si::kilograms * si::cubic_metres) * arg::_2 + arg::_2 * arg::_1 ); + + // correct for density with respect to STP + if(!opts_init.aerosol_independent_of_rhod) + thrust::transform( + rhod.begin(), rhod.end(), + arr.begin(), + arr.begin(), + arg::_1 / real_t(rho_stp() / si::kilograms * si::cubic_metres) * arg::_2 + ); + } + else // full-domain init 1D 2D 3D + { + // all SDs are initialized as if they were in 0-th cell, so set zero number in other cells + thrust::fill(arr.begin()+1, arr.end(), 0); + + // domain volume + thrust::transform( + arr.begin(), arr.begin()+1, + arr.begin(), + arg::_1 + * (opts_init.x1 - opts_init.x0) + * (n_dims >= 2 ? opts_init.z1 - opts_init.z0 : 1) + * (n_dims == 3 ? opts_init.y1 - opts_init.y0 : 1) + ); + + // correct for density with respect to STP; TODO: 0-th cell density is used, although later these SDs are spread around entire domain... + if(!opts_init.aerosol_independent_of_rhod) + thrust::transform( + rhod.begin(), rhod.begin()+1, + arr.begin(), + arr.begin(), + arg::_1 / real_t(rho_stp() / si::kilograms * si::cubic_metres) * arg::_2 + ); + } } template diff --git a/src/impl/particles_impl_init_n.ipp b/src/impl/particles_impl_init_n.ipp index 392c1d1c..48044fbf 100644 --- a/src/impl/particles_impl_init_n.ipp +++ b/src/impl/particles_impl_init_n.ipp @@ -93,7 +93,7 @@ namespace libcloudphxx // adjust to cell volume - if(n_dims > 0) + if(n_dims > 0 && !opts_init.domain_sd_init) { // rhod not needed anymore, reuse it to store dv on host thrust::copy( diff --git a/src/impl/particles_impl_init_sanity_check.ipp b/src/impl/particles_impl_init_sanity_check.ipp index 6f587771..5950fba2 100644 --- a/src/impl/particles_impl_init_sanity_check.ipp +++ b/src/impl/particles_impl_init_sanity_check.ipp @@ -141,6 +141,11 @@ namespace libcloudphxx throw std::runtime_error("libcloudph++: rlx_timescale <= 0"); if(opts_init.rlx_switch && opts_init.chem_switch) throw std::runtime_error("libcloudph++: CCN relaxation does not work with chemistry"); + + if(opts_init.domain_sd_init && opts_init.dry_sizes.size() > 0) + throw std::runtime_error("libcloudph++: entire domain initialization (domain_sd_init==true) does not work with initializing from a dry sizes dictionary (dry_sizes.size > 0)"); + if(opts_init.domain_sd_init && opts_init.sd_const_multi > 0) + throw std::runtime_error("libcloudph++: entire domain initialization (domain_sd_init==true) does not work with the constant-mutliplicity init (sd_const_multi > 0)"); } }; }; diff --git a/src/impl/particles_impl_init_wet.ipp b/src/impl/particles_impl_init_wet.ipp index 77646cd4..9480b5d6 100644 --- a/src/impl/particles_impl_init_wet.ipp +++ b/src/impl/particles_impl_init_wet.ipp @@ -44,7 +44,7 @@ namespace libcloudphxx { // initialising values of rw2 { - // calculating rw_eq + // calculating rw_eq { typedef thrust::permutation_iterator< typename thrust_device::vector::iterator, @@ -67,11 +67,11 @@ namespace libcloudphxx pi_t(T.begin(), ijk.begin() + n_part_old) )); - thrust::transform( - zip_it, zip_it + n_part_to_init, // input - rw2.begin() + n_part_old, // output + thrust::transform( + zip_it, zip_it + n_part_to_init, // input + rw2.begin() + n_part_old, // output detail::rw2_eq(opts_init.RH_max) - ); + ); } } } diff --git a/src/impl/particles_impl_init_xyz.ipp b/src/impl/particles_impl_init_xyz.ipp index 23978d31..cbda9b0c 100644 --- a/src/impl/particles_impl_init_xyz.ipp +++ b/src/impl/particles_impl_init_xyz.ipp @@ -29,7 +29,7 @@ namespace libcloudphxx using std::min; using std::max; #endif - + return u01 * min(p1, (ii+1) * dp) + (1. - u01) * max(p0, ii * dp); } }; @@ -55,17 +55,27 @@ namespace libcloudphxx // tossing random numbers rand_u01(n_part_to_init); - // shifting from [0,1] to random position within respective cell - // TODO: now the rand range is [0,1), include this here + // shifting from [0,1] to random position within respective cell + // TODO: now the rand range is [0,1), include this here + if(!opts_init.domain_sd_init) { - namespace arg = thrust::placeholders; - thrust::transform( - u01.begin(), - u01.begin() + n_part_to_init, + thrust::transform( + u01.begin(), + u01.begin() + n_part_to_init, ii[ix]->begin() + n_part_old, - v[ix]->begin() + n_part_old, + v[ix]->begin() + n_part_old, detail::pos_lgrngn_domain(a[ix], b[ix], d[ix]) - ); + ); + } + else // random position in the entire domain + { + namespace arg = thrust::placeholders; + thrust::transform( + u01.begin(), + u01.begin() + n_part_to_init, + v[ix]->begin() + n_part_old, + a[ix] + arg::_1 * (b[ix] - a[ix]) + ); } } }