diff --git a/bindings/python/lgrngn.hpp b/bindings/python/lgrngn.hpp index 53394bfa..52b9efb1 100644 --- a/bindings/python/lgrngn.hpp +++ b/bindings/python/lgrngn.hpp @@ -234,11 +234,11 @@ namespace libcloudphxx return *arg->opts_init; } - // set dry distros from a dict with (kappa, soluble_fraction) as key + // set dry distros from a dict with (kappa, soluble_fraction, sd_conc, sd_const_multi) as key template void set_dd( lgr::opts_init_t *arg, - const bp::dict &kappa_func) // a dict keyed by (kappa, soluble_fraction) + const bp::dict &kappa_func) // a dict keyed by (kappa, soluble_fraction, sd_conc, sd_const_multi) { arg->dry_distros.clear(); for (int i = 0; i < len(kappa_func.keys()); ++i) @@ -246,34 +246,35 @@ namespace libcloudphxx bp::tuple key = bp::extract(kappa_func.keys()[i]); const real_t kappa = bp::extract(key[0]); const real_t soluble_fraction = bp::extract(key[1]); + const unsigned long long sd_conc = bp::extract(key[2]); + const unsigned long long sd_const_multi = bp::extract(key[3]); arg->dry_distros.emplace( - libcloudphxx::lgrngn::kappa_soluble_fraction_t{kappa, soluble_fraction}, - std::make_shared>(kappa_func.values()[i]) + std::make_tuple(kappa, soluble_fraction, sd_conc, sd_const_multi), + std::static_pointer_cast>( + std::make_shared>(kappa_func.values()[i])) ); } } - // src_dry_distros moved from opts_init to opts + // src_dry_distros moved from opts_init to opts; key is (kappa, soluble_fraction, sd_conc, sd_const_multi, supstp) template void set_sdd( // src_dry_distro lgr::opts_t *arg, - const bp::dict &kappa_func) // a dict keyed by (kappa, soluble_fraction) + const bp::dict &kappa_func) // a dict keyed by (kappa, soluble_fraction, sd_conc, sd_const_multi, supstp) { arg->src_dry_distros.clear(); for (int i = 0; i < len(kappa_func.keys()); ++i) { bp::tuple key = bp::extract(kappa_func.keys()[i]); - bp::tuple val = bp::extract(kappa_func.values()[i]); const real_t kappa = bp::extract(key[0]); const real_t soluble_fraction = bp::extract(key[1]); - - const int sd_conc = bp::extract(val[1]); - const int supstp = bp::extract(val[2]); + const int sd_conc = bp::extract(key[2]); + const int sd_const_multi = bp::extract(key[3]); + const int supstp = bp::extract(key[4]); arg->src_dry_distros.emplace( - libcloudphxx::lgrngn::kappa_soluble_fraction_t{kappa, soluble_fraction}, - std::make_tuple(std::static_pointer_cast>( - std::make_shared>(val[0])), - sd_conc, supstp) + std::make_tuple(kappa, soluble_fraction, sd_conc, sd_const_multi, supstp), + std::static_pointer_cast>( + std::make_shared>(kappa_func.values()[i])) ); } } @@ -311,7 +312,7 @@ namespace libcloudphxx const int count = bp::extract (conc_count_list[1]); size_conc_map[bp::extract(size_conc.keys()[i])] = std::make_pair(conc, count); } - arg->dry_sizes[libcloudphxx::lgrngn::kappa_soluble_fraction_t{kappa, soluble_fraction}] = size_conc_map; + arg->dry_sizes[std::make_tuple(kappa, soluble_fraction)] = size_conc_map; } } @@ -319,36 +320,36 @@ namespace libcloudphxx template void set_sds( // src_dry_sizes lgr::opts_t *arg, - const bp::dict &kappa_func // a dict keyed by (kappa, soluble_fraction) + const bp::dict &kappa_func // a dict keyed by (kappa, soluble_fraction, supstp) ) { arg->src_dry_sizes.clear(); if(len(kappa_func.keys()) == 0) return; - // loop over kappas and soluble_fraction + // loop over kappa, soluble_fraction, and source intervals for (int j = 0; j < len(kappa_func.keys()); ++j) { - // extract the key tuple (kappa, soluble_fraction) + // extract the key tuple (kappa, soluble_fraction, supstp) const bp::tuple key = bp::extract(kappa_func.keys()[j]); const real_t kappa = bp::extract(key[0]); const real_t soluble_fraction = bp::extract(key[1]); + const int supstp = bp::extract(key[2]); - // extract size : {conc, count, supstp_src} dict for this (kappa, soluble_fraction) + // extract size : {conc, count} dict for this (kappa, soluble_fraction, supstp) const bp::dict size_conc = bp::extract(kappa_func.values()[j]); - std::map> size_conc_map; + std::map> size_conc_map; - // turn the size : {conc, count, supstp_src} dict into a size : {conc, count, supstp_src} map + // turn the size : {conc, count} dict into a size : {conc, count} map for (int i = 0; i < len(size_conc.keys()); ++i) { const bp::list conc_count_list = bp::extract(size_conc.values()[i]); - assert(len(conc_count_list) == 3); + assert(len(conc_count_list) == 2); const real_t conc = bp::extract(conc_count_list[0]); const int count = bp::extract (conc_count_list[1]); - const int supstp = bp::extract (conc_count_list[2]); - size_conc_map[bp::extract(size_conc.keys()[i])] = std::make_tuple(conc, count, supstp); + size_conc_map[bp::extract(size_conc.keys()[i])] = std::make_pair(conc, count); } - arg->src_dry_sizes[libcloudphxx::lgrngn::kappa_soluble_fraction_t{kappa, soluble_fraction}] = size_conc_map; + arg->src_dry_sizes[std::make_tuple(kappa, soluble_fraction, supstp)] = size_conc_map; } } diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index 949f94ab..2de69c6b 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -331,10 +331,8 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("supstp_rlx", &lgr::opts_init_t::supstp_rlx) .def_readwrite("kernel", &lgr::opts_init_t::kernel) .def_readwrite("adve_scheme", &lgr::opts_init_t::adve_scheme) - .def_readwrite("sd_conc", &lgr::opts_init_t::sd_conc) .def_readwrite("sd_conc_large_tail", &lgr::opts_init_t::sd_conc_large_tail) .def_readwrite("aerosol_independent_of_rhod", &lgr::opts_init_t::aerosol_independent_of_rhod) - .def_readwrite("sd_const_multi", &lgr::opts_init_t::sd_const_multi) .def_readwrite("rlx_bins", &lgr::opts_init_t::rlx_bins) .def_readwrite("rlx_sd_per_bin", &lgr::opts_init_t::rlx_sd_per_bin) .def_readwrite("rlx_timescale", &lgr::opts_init_t::rlx_timescale) diff --git a/docs/API.md b/docs/API.md index bce276f5..52702303 100644 --- a/docs/API.md +++ b/docs/API.md @@ -76,7 +76,6 @@ opts_init.nz = 64; // grid cells in z opts_init.dx = 100.0; // cell size x [m] opts_init.dy = 100.0; // cell size y [m] opts_init.dz = 100.0; // cell size z [m] -opts_init.sd_conc = 64; // super-droplets per cell // Create particle system auto particles = libcloudphxx::lgrngn::factory( @@ -439,7 +438,7 @@ arrinfo_t info_3d_v(arr_3d, str); See [USER_OPTIONS.md](USER_OPTIONS.md#lagrangian-lgrngn-initialization-options) for complete documentation of initialization options including: - Domain configuration (nx, ny, nz, dx, dy, dz) -- Super-droplet configuration (sd_conc, sd_conc_mean) +- Super-droplet configuration (per-distribution sd_conc, sd_const_multi) - Aerosol size distributions (dry_distros, dry_sizes) - Chemistry options (chem_switch, chem_rho) - GPU settings (dev_count, dev_id) @@ -954,7 +953,6 @@ int main() { opts_init.dt = dt; opts_init.nx = nx; opts_init.ny = ny; opts_init.nz = nz; opts_init.dx = dx; opts_init.dy = dy; opts_init.dz = dz; - opts_init.sd_conc = 64; opts_init.dry_distros = {{0.04e-6, 60e6}}; // 40 nm mode opts_init.kappa = 0.61; diff --git a/docs/USER_OPTIONS.md b/docs/USER_OPTIONS.md index bbd685d9..4277f5d7 100644 --- a/docs/USER_OPTIONS.md +++ b/docs/USER_OPTIONS.md @@ -135,8 +135,7 @@ These options are set once at initialization and cannot be changed during the si | Option | Type | Default | Description | |--------|------|---------|-------------| -| `sd_conc` | `unsigned long long` | `0` | Number of super-droplets per cell | -| `sd_const_multi` | `unsigned long long` | `0` | Alternative to `sd_conc`: constant multiplicity for all SDs | +| `dry_distros[].sd_const_multi` | `unsigned long long` | `0` | Constant multiplicity for the distribution; cannot be used with that distribution's `sd_conc` | | `n_sd_max` | `unsigned long long` | `0` | Maximum number of super-droplets in the system (should account for sources) | | `sd_conc_large_tail` | `bool` | `false` | Add more SDs to better represent large tail of the distribution | | `rd_min`, `rd_max` | `real_t` | `-1` | Min/max dry radius of droplets [m]; negative = auto-detect | @@ -147,8 +146,8 @@ Two methods are available: **1. Distribution-based (recommended):** ```cpp -typedef std::unordered_map< - kappa_soluble_fraction_t, // (kappa - hygroscopicity parameter, soluble_fraction - volume fraction of soluble part) +typedef std::map< + std::tuple, // kappa, soluble_fraction, sd_conc, sd_const_multi std::shared_ptr> // n(ln(rd)) @ STP > dry_distros_t; dry_distros_t dry_distros; @@ -157,7 +156,7 @@ dry_distros_t dry_distros; **2. Size-number pairs:** ```cpp typedef std::map< - kappa_soluble_fraction_t // (kappa - hygroscopicity parameter, soluble_fraction - volume fraction of soluble part) + std::tuple, // kappa, soluble_fraction std::map // STP concentration [1/m^3], number of SDs > @@ -275,13 +274,11 @@ dry_sizes_t dry_sizes; | Option | Type | Default | Description | |--------|------|---------|-------------| | `src_type` | `src_t` | `off` | Type of CCN source (`off`, `constant`, etc.) | -| `src_dry_distros` | `dry_distros_t` | - | Source distribution per unit time | +| `src_dry_distros` | `src_dry_distros_t` | - | Source distribution per unit time; map key is `(kappa, soluble_fraction, sd_conc, sd_const_multi, supstp)` and simple sources require exactly one of `sd_conc` or `sd_const_multi` to be nonzero | | `src_dry_sizes` | `dry_sizes_t` | - | Alternative source specification using size-number pairs | -| `src_sd_conc` | `unsigned long long` | `0` | Number of SDs created per cell per source iteration | | `src_x0`, `src_x1` | `real_t` | `0` | Source box x-boundaries [m] (rounded to cell boundaries) | | `src_y0`, `src_y1` | `real_t` | `0` | Source box y-boundaries [m] | | `src_z0`, `src_z1` | `real_t` | `0` | Source box z-boundaries [m] | -| `supstp_src` | `int` | `1` | Timestep interval for applying source | #### Aerosol Relaxation Options @@ -385,7 +382,6 @@ opts_init.ny = 100; opts_init.nz = 100; opts_init.dx = opts_init.dy = opts_init.dz = 10; // 10 m grid spacing opts_init.dt = 1.0; // 1 s timestep -opts_init.sd_conc = 64; // 64 SDs per cell opts_init.sstp_cond = 4; // 4 condensation substeps opts_init.adaptive_sstp_cond = true; // Enable adaptive substepping opts_init.exact_sstp_cond = true; // Per-particle substepping @@ -400,7 +396,7 @@ auto lognormal = [](double lnr) { return n_tot * exp(-pow((lnr - log(mean_r)), 2) / 2 / pow(log(stdev), 2)) / log(stdev) / sqrt(2 * M_PI); }; -opts_init.dry_distros[0.61] = std::make_shared(lognormal); +opts_init.dry_distros[{0.61, 1., 64, 0}] = std::make_shared(lognormal); // kappa, soluble_fraction, sd_conc, sd_const_multi // Runtime options (can change each step) libcloudphxx::lgrngn::opts_t opts; diff --git a/include/libcloudph++/lgrngn/distro_t.hpp b/include/libcloudph++/lgrngn/distro_t.hpp index b09ab21c..b4d09211 100644 --- a/include/libcloudph++/lgrngn/distro_t.hpp +++ b/include/libcloudph++/lgrngn/distro_t.hpp @@ -6,35 +6,19 @@ namespace libcloudphxx { using common::unary_function; - template - struct kappa_soluble_fraction_t { - real_t kappa; - real_t soluble_fraction; // volume fraction of the soluble part (0-1) - - kappa_soluble_fraction_t(real_t kappa_, real_t soluble_fraction_) - : kappa(kappa_), soluble_fraction(soluble_fraction_) - {} - - bool operator<(const kappa_soluble_fraction_t &other) const - { - if (kappa != other.kappa) return kappa < other.kappa; - return soluble_fraction < other.soluble_fraction; - } - }; - // initial dry sizes of aerosol // defined with a distribution // uses shared_ptr to make opts_init copyable template using dry_distros_t = std::map< - kappa_soluble_fraction_t, // (kappa, soluble_fraction); dry_distros defines total dry radius + std::tuple, // kappa, soluble_fraction, sd_conc, sd_const_multi std::shared_ptr> // n(ln(rd)) @ STP; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true >; // defined with a size-number pair template using dry_sizes_t = std::map< - kappa_soluble_fraction_t, // (kappa, soluble_fraction); dry_sizes defines total dry radius + std::tuple, // (kappa, soluble_fraction); dry_sizes defines total dry radius std::map // STP_concentration [1/m^3], number of SD that represent this radius kappa and concentration > @@ -43,16 +27,16 @@ namespace libcloudphxx // similar, but for sources of aerosols after initialization template using src_dry_distros_t = std::map< - kappa_soluble_fraction_t, // (kappa, soluble_fraction); src_dry_distros defines total dry radius - std::tuple>, int, int> // 1st: n(ln(rd)) @ STP created per second; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true; 2nd: sd_conc for this distribution ; 3rd: supstp for this aerosol (interval in timesteps beween addition of these aerosols) + std::tuple, // kappa, soluble_fraction, sd_conc, sd_const_multi, supstp + std::shared_ptr> // n(ln(rd)) @ STP created per second; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true >; // defined with a size-number pair template using src_dry_sizes_t = std::map< - kappa_soluble_fraction_t, // (kappa, soluble_fraction); src_dry_sizes defines total dry radius + std::tuple, // kappa, soluble_fraction, supstp std::map // STP_concentration [1/m^3] created per second, number of SD that represent this radius kappa and concentration, supstp + std::pair // STP_concentration [1/m^3] created per second, number of SD that represent this radius kappa and concentration > >; diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index f4ae0ee6..2d43a3e7 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -52,9 +52,6 @@ namespace libcloudphxx // Lagrangian domain extents real_t x0, y0, z0, x1, y1, z1; - // no. of super-droplets per cell - unsigned long long sd_conc; - // should more SDs be added to better represent large tail of the distribution bool sd_conc_large_tail; @@ -65,9 +62,6 @@ 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 - unsigned long long sd_const_multi; - // max no. of super-droplets in the system // should be enough to store particles from sources unsigned long long n_sd_max; @@ -187,10 +181,8 @@ namespace libcloudphxx dx(1), dy(1), dz(1), x0(0), y0(0), z0(0), x1(1), y1(1), z1(1), - sd_conc(0), sd_conc_large_tail(false), aerosol_independent_of_rhod(false), - sd_const_multi(0), dt(0), sstp_cond(1), sstp_coal(1), sstp_chem(1), sstp_cond_act(1), chem_switch(false), // chemical reactions turned off by default diff --git a/src/impl/initialization/particles_impl_init_SD_with_distros.ipp b/src/impl/initialization/particles_impl_init_SD_with_distros.ipp index 001c2cdc..67ddec0a 100644 --- a/src/impl/initialization/particles_impl_init_SD_with_distros.ipp +++ b/src/impl/initialization/particles_impl_init_SD_with_distros.ipp @@ -14,35 +14,26 @@ namespace libcloudphxx template void particles_t::impl::init_SD_with_distros() { - // calc sum of ln(rd) ranges of all distributions - real_t tot_lnrd_rng = 0.; - if(opts_init.sd_conc > 0) - for (auto ddi = opts_init.dry_distros.cbegin(); ddi != opts_init.dry_distros.cend(); ++ddi) - { - init_dist_analysis_sd_conc( - *(ddi->second), - opts_init.sd_conc - ); - tot_lnrd_rng += log_rd_max - log_rd_min; - } - // initialize SDs of each kappa-type for (auto ddi = opts_init.dry_distros.cbegin(); ddi != opts_init.dry_distros.cend(); ++ddi) { - if(opts_init.sd_conc > 0) + const auto &distro = ddi->second; + const auto sd_conc = std::get<2>(ddi->first); + const auto sd_const_multi = std::get<3>(ddi->first); + if(sd_conc > 0) { - init_SD_with_distros_sd_conc(*(ddi->second), tot_lnrd_rng); + init_SD_with_distros_sd_conc(*distro, sd_conc); init_SD_with_distros_finalize(ddi->first); if(opts_init.sd_conc_large_tail) { - init_SD_with_distros_tail(*(ddi->second), log_rd_max); + init_SD_with_distros_tail(*distro, log_rd_max); init_SD_with_distros_finalize(ddi->first); } } - if(opts_init.sd_const_multi > 0) + if(sd_const_multi > 0) { - init_SD_with_distros_const_multi(*(ddi->second)); + init_SD_with_distros_const_multi(*distro, sd_const_multi); init_SD_with_distros_finalize(ddi->first); } } @@ -50,16 +41,16 @@ namespace libcloudphxx // final inits common for tail/sd_conc/const_multi template - void particles_t::impl::init_SD_with_distros_finalize(const kappa_soluble_fraction_t &kpa_sol_frac, const bool unravel_ijk_switch) + void particles_t::impl::init_SD_with_distros_finalize(const std::tuple &kpa_sol_frac, const bool unravel_ijk_switch) { // dry_distros defines total dry radius; insoluble part determined by soluble_fraction // init kappa - init_kappa(kpa_sol_frac.kappa, kpa_sol_frac.soluble_fraction); + init_kappa(std::get<0>(kpa_sol_frac), std::get<1>(kpa_sol_frac)); if (opts_init.ice_switch) { - init_insol(kpa_sol_frac.soluble_fraction); + init_insol(std::get<1>(kpa_sol_frac)); init_a_c_rho_ice(); if (! opts_init.time_dep_ice_nucl) diff --git a/src/impl/initialization/particles_impl_init_SD_with_distros_const_multi.ipp b/src/impl/initialization/particles_impl_init_SD_with_distros_const_multi.ipp index 8e365756..26834574 100644 --- a/src/impl/initialization/particles_impl_init_SD_with_distros_const_multi.ipp +++ b/src/impl/initialization/particles_impl_init_SD_with_distros_const_multi.ipp @@ -11,7 +11,7 @@ namespace libcloudphxx namespace lgrngn { template - void particles_t::impl::init_SD_with_distros_const_multi(const common::unary_function &fun) + void particles_t::impl::init_SD_with_distros_const_multi(const common::unary_function &fun, const n_t sd_const_multi) { // analyze the distribution, TODO: just did it init_dist_analysis_const_multi(fun); @@ -19,7 +19,7 @@ namespace libcloudphxx throw std::runtime_error(detail::formatter() << "Distribution analysis error: rd_min(" << exp(log_rd_min) << ") >= rd_max(" << exp(log_rd_max) << ")"); // init number of SDs of this kappa in cells, TODO: due to rounding, we might end up with not exactly sd_conc SDs per cell... - init_count_num_const_multi(fun); + init_count_num_const_multi(fun, sd_const_multi); // update no of particles // TODO: move to a separate function @@ -35,7 +35,7 @@ namespace libcloudphxx init_dry_const_multi(fun); // init multiplicities - init_n_const_multi(opts_init.sd_const_multi); + init_n_const_multi(sd_const_multi); } }; }; diff --git a/src/impl/initialization/particles_impl_init_SD_with_distros_sd_conc.ipp b/src/impl/initialization/particles_impl_init_SD_with_distros_sd_conc.ipp index 29921f04..f26f37a9 100644 --- a/src/impl/initialization/particles_impl_init_SD_with_distros_sd_conc.ipp +++ b/src/impl/initialization/particles_impl_init_SD_with_distros_sd_conc.ipp @@ -11,21 +11,18 @@ namespace libcloudphxx namespace lgrngn { template - void particles_t::impl::init_SD_with_distros_sd_conc(const common::unary_function &fun, const real_t &tot_lnrd_rng) + void particles_t::impl::init_SD_with_distros_sd_conc(const common::unary_function &fun, const n_t sd_conc) { // analyze the distribution, TODO: just did it in init_SD_with_distros init_dist_analysis_sd_conc( fun, - opts_init.sd_conc + sd_conc ); if(log_rd_min >= log_rd_max) throw std::runtime_error(detail::formatter() << "Distribution analysis error: rd_min(" << exp(log_rd_min) << ") >= rd_max(" << exp(log_rd_max) << ")"); - // init number of SDs of this kappa in cells, TODO: due to rounding, we might end up with not exactly sd_conc SDs per cell... - // adjust the multiplicity init coefficient to smaller number of SDs representing this kappa-type - real_t fraction = (log_rd_max - log_rd_min) / tot_lnrd_rng; - multiplier *= opts_init.sd_conc / int(fraction * opts_init.sd_conc + 0.5); - init_count_num_sd_conc(fraction); + // init number of SDs of this kappa in each cell + init_count_num_sd_conc(sd_conc); // update no of particles // TODO: move to a separate function diff --git a/src/impl/initialization/particles_impl_init_SD_with_sizes.ipp b/src/impl/initialization/particles_impl_init_SD_with_sizes.ipp index f14283ab..9e297c0f 100644 --- a/src/impl/initialization/particles_impl_init_SD_with_sizes.ipp +++ b/src/impl/initialization/particles_impl_init_SD_with_sizes.ipp @@ -22,8 +22,8 @@ namespace libcloudphxx // loop over (kappa, soluble_fraction) pairs for (auto dsi = opts_init.dry_sizes.cbegin(); dsi != opts_init.dry_sizes.cend(); ++dsi) { - const real_t &kappa(dsi->first.kappa); - const real_t &soluble_fraction(dsi->first.soluble_fraction); + const real_t &kappa(std::get<0>(dsi->first)); + const real_t &soluble_fraction(std::get<1>(dsi->first)); const auto &size_number_map(dsi->second); // loop over the "size : {concentration, count}" pairs for this (kappa, soluble_fraction) pair diff --git a/src/impl/initialization/particles_impl_init_count_num.ipp b/src/impl/initialization/particles_impl_init_count_num.ipp index f2aa1710..41ef0e83 100644 --- a/src/impl/initialization/particles_impl_init_count_num.ipp +++ b/src/impl/initialization/particles_impl_init_count_num.ipp @@ -29,9 +29,9 @@ namespace libcloudphxx // init number of SDs to be initialized per cell template - void particles_t::impl::init_count_num_sd_conc(const real_t &ratio) + void particles_t::impl::init_count_num_sd_conc(const n_t sd_conc) { - thrust::fill(count_num.begin(), count_num.end(), ratio * opts_init.sd_conc); + thrust::fill(count_num.begin(), count_num.end(), sd_conc); } // calculate number of droplets in a cell from concentration [1/m^3], taking into account cell volume and air density @@ -97,15 +97,6 @@ namespace libcloudphxx //init_count_num_hlpr(conc_multi.first, conc_multi.second); } - template - void particles_t::impl::init_count_num_const_multi( - const common::unary_function &n_of_lnrd_stp - ) - { - const real_t integral = detail::integrate(n_of_lnrd_stp, log_rd_min, log_rd_max, config.bin_precision); - init_count_num_hlpr(integral, opts_init.sd_const_multi); - } - template void particles_t::impl::init_count_num_const_multi( const common::unary_function &n_of_lnrd_stp, @@ -167,5 +158,52 @@ namespace libcloudphxx break; } } + + template + void particles_t::impl::init_count_num_src_const_multi( + const common::unary_function &n_of_lnrd_stp, + const thrust_size_t &const_multi, + const real_t &dt + ) + { + const real_t integral = detail::integrate(n_of_lnrd_stp, log_rd_min, log_rd_max, config.bin_precision) * dt; + init_count_num_hlpr(integral, const_multi); + + namespace arg = thrust::placeholders; + const thrust_size_t i0 = opts_init.src_x0 / opts_init.dx + 0.5; + const thrust_size_t i1 = opts_init.src_x1 / opts_init.dx + 0.5; + const thrust_size_t j0 = opts_init.src_y0 / opts_init.dy + 0.5; + const thrust_size_t j1 = opts_init.src_y1 / opts_init.dy + 0.5; + const thrust_size_t k0 = opts_init.src_z0 / opts_init.dz + 0.5; + const thrust_size_t k1 = opts_init.src_z1 / opts_init.dz + 0.5; + + switch(n_dims) + { + case 0 : throw std::runtime_error("libcloudph++: init_count_num_src_const_multi called in 0D"); + case 1 : throw std::runtime_error("libcloudph++: init_count_num_src_const_multi called in 1D"); + case 2: + thrust::transform( + zero, zero + n_cell, count_num.begin(), count_num.begin(), + arg::_2 + * ((arg::_1 % opts_init.nz) < k1) + * ((arg::_1 % opts_init.nz) >= k0) + * ((arg::_1 / opts_init.nz) < i1) + * ((arg::_1 / opts_init.nz) >= i0) + ); + break; + case 3: + thrust::transform( + zero, zero + n_cell, count_num.begin(), count_num.begin(), + arg::_2 + * ((arg::_1 % opts_init.nz) < k1) + * ((arg::_1 % opts_init.nz) >= k0) + * ((arg::_1 / (opts_init.nz * opts_init.ny)) < i1) + * ((arg::_1 / (opts_init.nz * opts_init.ny)) >= i0) + * (((arg::_1 / opts_init.nz) % opts_init.ny) < j1) + * (((arg::_1 / opts_init.nz) % opts_init.ny) >= j0) + ); + break; + } + } }; }; diff --git a/src/impl/initialization/particles_impl_init_sanity_check.ipp b/src/impl/initialization/particles_impl_init_sanity_check.ipp index 8b84abf1..b9258147 100644 --- a/src/impl/initialization/particles_impl_init_sanity_check.ipp +++ b/src/impl/initialization/particles_impl_init_sanity_check.ipp @@ -70,17 +70,32 @@ namespace libcloudphxx throw std::runtime_error("libcloudph++: Both dry_distros and dry_sizes are undefined"); for(auto &dist : opts_init.dry_distros) - if(dist.first.soluble_fraction < 0 || dist.first.soluble_fraction > 1) + if(std::get<1>(dist.first) < 0 || std::get<1>(dist.first) > 1) throw std::runtime_error("libcloudph++: soluble_fraction in dry_distros must be in [0, 1]"); for(auto &size : opts_init.dry_sizes) - if(size.first.soluble_fraction < 0 || size.first.soluble_fraction > 1) + if(std::get<1>(size.first) < 0 || std::get<1>(size.first) > 1) throw std::runtime_error("libcloudph++: soluble_fraction in dry_sizes must be in [0, 1]"); - if(opts_init.sd_conc_large_tail && opts_init.sd_conc == 0) - throw std::runtime_error("libcloudph++: Sd_conc_large_tail make sense only with sd_conc init (i.e. sd_conc>0)"); + bool has_sd_conc = false; + bool has_sd_const_multi = false; + for(const auto &dist : opts_init.dry_distros) + { + const auto sd_conc = std::get<2>(dist.first); + const auto sd_const_multi = std::get<3>(dist.first); + if(sd_conc > 0 && sd_const_multi > 0) + throw std::runtime_error("libcloudph++: specify either sd_conc or sd_const_multi for each dry distribution, not both"); + has_sd_conc = has_sd_conc || sd_conc > 0; + has_sd_const_multi = has_sd_const_multi || sd_const_multi > 0; + } + + if(opts_init.sd_conc_large_tail) + { + if(!has_sd_conc) + throw std::runtime_error("libcloudph++: sd_conc_large_tail requires sd_conc > 0 in at least one dry distribution"); + } - if(opts_init.sd_const_multi > 0 && opts_init.src_type!=src_t::off) + if(has_sd_const_multi && opts_init.src_type!=src_t::off) throw std::runtime_error("libcloudph++: aerosol source and constant multiplicity option are not compatible"); // NOTE: why not? if (n_dims > 0) @@ -101,8 +116,7 @@ namespace libcloudphxx } if (opts_init.dt == 0) throw std::runtime_error("libcloudph++: please specify opts_init.dt"); - if (opts_init.sd_conc * opts_init.sd_const_multi != 0) throw std::runtime_error("libcloudph++: specify either opts_init.sd_conc or opts_init.sd_const_multi, not both"); - if (opts_init.sd_conc == 0 && opts_init.sd_const_multi == 0 && opts_init.dry_sizes.size() == 0) throw std::runtime_error("libcloudph++: please specify opts_init.sd_conc, opts_init.sd_const_multi or opts_init.dry_sizes"); + if (!has_sd_conc && !has_sd_const_multi && opts_init.dry_sizes.size() == 0) throw std::runtime_error("libcloudph++: please specify sd_conc or sd_const_multi in dry_distros, or dry_sizes"); if (opts_init.coal_switch) { if(opts_init.terminal_velocity == vt_t::undefined) throw std::runtime_error("libcloudph++: please specify opts_init.terminal_velocity or turn off opts_init.coal_switch"); @@ -153,11 +167,11 @@ namespace libcloudphxx if(opts_init.chem_switch) { for(auto &dist : opts_init.dry_distros) - if(dist.first.soluble_fraction < 1) + if(std::get<1>(dist.first) < 1) throw std::runtime_error("libcloudph++: insoluble aerosol (defined in opts_init.dry_distros) does not work with chemistry"); for(auto &size : opts_init.dry_sizes) - if(size.first.soluble_fraction < 1) + if(std::get<1>(size.first) < 1) throw std::runtime_error("libcloudph++: insoluble aerosol (defined in opts_init.dry_sizes) does not work with chemistry"); } diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index f4c9e8f0..dcdf693e 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -383,7 +383,7 @@ namespace libcloudphxx sstp_chem(_opts_init.sstp_chem), // sstp_cond_act(std::max(_opts_init.sstp_cond_act, _opts_init.sstp_cond)), sstp_cond_act(_opts_init.sstp_cond_act), - pure_const_multi (((_opts_init.sd_conc) == 0) && (_opts_init.sd_const_multi > 0 || _opts_init.dry_sizes.size() > 0)), // coal prob can be greater than one only in sd_conc simulations + pure_const_multi (std::none_of(_opts_init.dry_distros.cbegin(), _opts_init.dry_distros.cend(), [](const auto &distro) { return std::get<2>(distro.first) > 0; }) && (std::any_of(_opts_init.dry_distros.cbegin(), _opts_init.dry_distros.cend(), [](const auto &distro) { return std::get<3>(distro.first) > 0; }) || _opts_init.dry_sizes.size() > 0)), // coal prob can be greater than one only in sd_conc simulations //tmp_device_real_part(6), tmp_host_real_part("tmp_host_real_part"), tmp_host_real_grid("tmp_host_real_grid"), @@ -547,10 +547,10 @@ namespace libcloudphxx void sanity_checks(); void init_SD_with_distros(); - void init_SD_with_distros_sd_conc(const common::unary_function &, const real_t &); + void init_SD_with_distros_sd_conc(const common::unary_function &, const n_t); void init_SD_with_distros_tail(const common::unary_function &, const real_t); - void init_SD_with_distros_const_multi(const common::unary_function &); - void init_SD_with_distros_finalize(const kappa_soluble_fraction_t &, const bool unravel_ijk = true); + void init_SD_with_distros_const_multi(const common::unary_function &, const n_t); + void init_SD_with_distros_finalize(const std::tuple &, const bool unravel_ijk = true); void init_SD_with_sizes(); void init_sanity_check( const arrinfo_t, const arrinfo_t, const arrinfo_t, @@ -587,12 +587,12 @@ namespace libcloudphxx void init_T_freeze(); void init_a_c_rho_ice(); void init_incloud_time(); - void init_count_num_sd_conc(const real_t & = 1); - void init_count_num_const_multi(const common::unary_function &); + void init_count_num_sd_conc(const n_t); void init_count_num_const_multi(const common::unary_function &, const thrust_size_t &); void init_count_num_dry_sizes(const std::pair &); void init_count_num_hlpr(const real_t &, const thrust_size_t &); void init_count_num_src(const thrust_size_t &); + void init_count_num_src_const_multi(const common::unary_function &, const thrust_size_t &, const real_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); diff --git a/src/impl/sources_and_relaxation_of_SDs/particles_impl_rlx_dry_distros.ipp b/src/impl/sources_and_relaxation_of_SDs/particles_impl_rlx_dry_distros.ipp index d4b8c1a7..290ae15d 100644 --- a/src/impl/sources_and_relaxation_of_SDs/particles_impl_rlx_dry_distros.ipp +++ b/src/impl/sources_and_relaxation_of_SDs/particles_impl_rlx_dry_distros.ipp @@ -303,7 +303,7 @@ namespace libcloudphxx n_part_to_init = n_part - n_part_old; // hskpng_resize_npart(); - init_SD_with_distros_finalize(lgrngn::kappa_soluble_fraction_t{kappa, real_t(1)}, false); // no need to unravel ijk there, because i j k are already initialized + init_SD_with_distros_finalize(std::make_tuple(kappa, real_t(1), 0ull, 0ull), false); // no need to unravel ijk there, because i j k are already initialized // TODO: we assume that relaxation produces water (not ice), and that soluble_fraction==1 for produced aerosols // TODO: asserts of newly added SD parameters? e.g. how many SD, how big is multiplicity etc. diff --git a/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros.ipp b/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros.ipp index d3c57a40..f4e4f7ce 100644 --- a/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros.ipp +++ b/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros.ipp @@ -19,7 +19,7 @@ namespace libcloudphxx if (sdd.size() > 1) throw std::runtime_error("libcloudph++: src_dry_distros can only have a single kappa value."); - if (opts_init.src_type == src_t::matching && !sdd.empty() && sdd.begin()->first.kappa != opts_init.dry_distros.begin()->first.kappa) + if (opts_init.src_type == src_t::matching && !sdd.empty() && std::get<0>(sdd.begin()->first) != std::get<0>(opts_init.dry_distros.begin()->first)) throw std::runtime_error("libcloudph++: For 'matching' CCN source, kappa of the source has to be the same as that of the initial profile (no kappa matching done)"); if(opts_init.src_type == src_t::matching) diff --git a/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros_matching.ipp b/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros_matching.ipp index 1b174bda..d13555b5 100644 --- a/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros_matching.ipp +++ b/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros_matching.ipp @@ -51,16 +51,18 @@ namespace libcloudphxx void particles_t::impl::src_dry_distros_matching(const src_dry_distros_t &sdd) { auto p_sdd = sdd.cbegin(); - assert(p_sdd->first.soluble_fraction == 1); // partial solubility in matching source not implemented yet + assert(std::get<1>(p_sdd->first) == 1); // partial solubility in matching source not implemented yet + if(get<3>(p_sdd->first) > 0) + throw std::runtime_error("libcloudph++: constant multiplicity source distributions are only supported with src_t::simple"); // add the source only once every number of steps - assert(get<2>(p_sdd->second) > 0); - if(src_stp_ctr % get<2>(p_sdd->second) != 0) return; + assert(get<4>(p_sdd->first) > 0); + if(src_stp_ctr % get<4>(p_sdd->first) != 0) return; - const real_t sup_dt = get<2>(p_sdd->second) * opts_init.dt; + const real_t sup_dt = get<4>(p_sdd->first) * opts_init.dt; // set number of SDs to init; use count_num as storage - init_count_num_src(get<1>(p_sdd->second)); + init_count_num_src(get<2>(p_sdd->first)); // --- sort already existing SDs; primary key ijk, secondary rd --- // TODO: do all of this only on SDs in cells below src_z1? @@ -102,8 +104,8 @@ namespace libcloudphxx // TODO: this could be done once at the beginning of the simulation // TODO: take rd_insol into account here? init_dist_analysis_sd_conc( - *get<0>(p_sdd->second), - get<1>(p_sdd->second), + *p_sdd->second, + get<2>(p_sdd->first), sup_dt ); @@ -221,12 +223,12 @@ namespace libcloudphxx // init other properties of SDs that didnt have a match init_kappa( - p_sdd->first.kappa, - p_sdd->first.soluble_fraction + std::get<0>(p_sdd->first), + std::get<1>(p_sdd->first) ); if (opts_init.ice_switch) { - init_insol(p_sdd->first.soluble_fraction); + init_insol(std::get<1>(p_sdd->first)); init_a_c_rho_ice(); if (! opts_init.time_dep_ice_nucl) { @@ -238,7 +240,7 @@ namespace libcloudphxx init_incloud_time(); init_n_sd_conc( - *get<0>(p_sdd->second) + *p_sdd->second ); // TODO: document that n_of_lnrd_stp is expected! // init rw @@ -390,7 +392,7 @@ namespace libcloudphxx // init n of the copied SDs, but using the src distribution init_n_sd_conc( - *get<0>(p_sdd->second) + *p_sdd->second ); // TODO: document that n_of_lnrd_stp is expected! // add the just-initialized multiplicities to the old ones diff --git a/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros_simple.ipp b/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros_simple.ipp index 5262cbd6..9f8b69f4 100644 --- a/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros_simple.ipp +++ b/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_distros_simple.ipp @@ -22,25 +22,33 @@ namespace libcloudphxx // TODO: add a loop to allow sdd.size>1 auto p_sdd = sdd.cbegin(); - if(p_sdd->first.soluble_fraction < 0 || p_sdd->first.soluble_fraction > 1) + if(std::get<1>(p_sdd->first) < 0 || std::get<1>(p_sdd->first) > 1) throw std::runtime_error("libcloudph++: soluble_fraction in opts.src_dry_distros must be in [0, 1]"); - // add the source only once every number of steps - assert(get<2>(p_sdd->second) > 0); - if(src_stp_ctr % get<2>(p_sdd->second) != 0) return; + const int sd_conc = get<2>(p_sdd->first); + const int sd_const_multi = get<3>(p_sdd->first); + const int supstp = get<4>(p_sdd->first); + if(sd_conc > 0 && sd_const_multi > 0) + throw std::runtime_error("libcloudph++: specify either sd_conc or sd_const_multi for each source dry distribution, not both"); - const real_t sup_dt = get<2>(p_sdd->second) * opts_init.dt; + // add the source only once every number of steps + assert(supstp > 0); + if(src_stp_ctr % supstp != 0) return; - // set number of SDs to init; use count_num as storage - init_count_num_src(get<1>(p_sdd->second)); + const real_t sup_dt = supstp * opts_init.dt; - // analyze distribution to get rd_min and max needed for bin sizes - // TODO: this could be done once at the beginning of the simulation - init_dist_analysis_sd_conc( - *(get<0>(p_sdd->second)), - get<1>(p_sdd->second), - sup_dt - ); + if(sd_conc > 0) + { + init_count_num_src(sd_conc); + init_dist_analysis_sd_conc(*p_sdd->second, sd_conc, sup_dt); + } + else if(sd_const_multi > 0) + { + init_dist_analysis_const_multi(*p_sdd->second); + init_count_num_src_const_multi(*p_sdd->second, sd_const_multi, sup_dt); + } + else + throw std::runtime_error("libcloudph++: specify either sd_conc or sd_const_multi for each source dry distribution"); namespace arg = thrust::placeholders; @@ -55,20 +63,25 @@ namespace libcloudphxx // init ijk and rd3 of new particles init_ijk(); - init_dry_sd_conc(); - - init_n_sd_conc( - *get<0>(p_sdd->second) - ); + if(sd_conc > 0) + { + init_dry_sd_conc(); + init_n_sd_conc(*p_sdd->second); + } + else + { + init_dry_const_multi(*p_sdd->second); + init_n_const_multi(sd_const_multi); + } // init other properties of SDs init_kappa( - p_sdd->first.kappa, - p_sdd->first.soluble_fraction + std::get<0>(p_sdd->first), + std::get<1>(p_sdd->first) ); if (opts_init.ice_switch) { - init_insol(p_sdd->first.soluble_fraction); + init_insol(std::get<1>(p_sdd->first)); init_a_c_rho_ice(); if (! opts_init.time_dep_ice_nucl) { diff --git a/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_sizes.ipp b/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_sizes.ipp index 130c190a..b7ec880a 100644 --- a/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_sizes.ipp +++ b/src/impl/sources_and_relaxation_of_SDs/particles_impl_src_dry_sizes.ipp @@ -15,25 +15,26 @@ namespace libcloudphxx void particles_t::impl::src_dry_sizes(const src_dry_sizes_t &sds) { - // loop over (kappa, soluble_fraction) pairs + // loop over (kappa, soluble_fraction, supstp) tuples // for (typename dry_sizes_t::const_iterator dsi = opts.src_dry_sizes.begin(); dsi != opts.src_dry_sizes.end(); ++dsi) for (auto dsi = sds.cbegin(); dsi != sds.cend(); ++dsi) { - const real_t &kappa(dsi->first.kappa); - const real_t &soluble_fraction(dsi->first.soluble_fraction); + const real_t &kappa(std::get<0>(dsi->first)); + const real_t &soluble_fraction(std::get<1>(dsi->first)); + const int supstp(std::get<2>(dsi->first)); const auto &size_number_map(dsi->second); if(soluble_fraction < 0 || soluble_fraction > 1) throw std::runtime_error("libcloudph++: soluble_fraction in opts.src_dry_sizes must be in [0, 1]"); - // loop over the "size : {concentration per second, multiplicity, supstp}" for this (kappa, soluble_fraction) pair + // loop over the "size : {concentration per second, multiplicity}" for this (kappa, soluble_fraction, supstp) tuple for (auto sni = size_number_map.cbegin(); sni != size_number_map.cend(); ++sni) { // add the source only once every number of steps - assert(get<2>(sni->second) > 0); - if(src_stp_ctr % get<2>(sni->second) != 0) continue; + assert(supstp > 0); + if(src_stp_ctr % supstp != 0) continue; - const real_t sup_dt = get<2>(sni->second) * opts_init.dt; + const real_t sup_dt = supstp * opts_init.dt; // init number of SDs of this kappa in cells init_count_num_src(get<1>(sni->second)); @@ -65,7 +66,7 @@ namespace libcloudphxx } // init multiplicities - init_n_dry_sizes(get<0>(sni->second)*sup_dt, get<1>(sni->second)); + init_n_dry_sizes(sni->second.first * sup_dt, sni->second.second); // initialising wet radii init_wet(); diff --git a/tests/mpi/mpi_adve_test.cpp b/tests/mpi/mpi_adve_test.cpp index b6c50923..0cd3f8d7 100644 --- a/tests/mpi/mpi_adve_test.cpp +++ b/tests/mpi/mpi_adve_test.cpp @@ -90,8 +90,7 @@ void test(backend_t backend, std::string back_name, int ndims, bool dir, int n_d //opts_init.nx = nx_min; //int nx_total = nx_min * size; opts_init.x1 = opts_init.nx * opts_init.dx;// nx_factor/2*(rank/2+1); - opts_init.sd_conc = 64; - opts_init.n_sd_max = 1000*opts_init.sd_conc; + opts_init.n_sd_max = 1000 * 64; opts_init.rng_seed = 4444 + rank; if(ndims>1) { @@ -108,7 +107,6 @@ void test(backend_t backend, std::string back_name, int ndims, bool dir, int n_d opts_init.dev_id = rank%n_devices; //opts_init.dev_id = rank; std::cout << "device id: " << opts_init.dev_id << std::endl; -// opts_init.sd_const_multi = 1; /* boost::assign::ptr_map_insert< @@ -121,7 +119,7 @@ void test(backend_t backend, std::string back_name, int ndims, bool dir, int n_d */ opts_init.dry_distros.emplace( - libcloudphxx::lgrngn::kappa_soluble_fraction_t{double(0.001), double(1.)}, // kappa, soluble_fraction + std::make_tuple(double(0.001), double(1.), 64ull, 0ull), // kappa, soluble_fraction, sd_conc, sd_const_multi std::make_shared>() // distribution ); diff --git a/tests/python/physics/coalescence_golovin.py b/tests/python/physics/coalescence_golovin.py index aaa8d0ae..40b649ad 100644 --- a/tests/python/physics/coalescence_golovin.py +++ b/tests/python/physics/coalescence_golovin.py @@ -70,7 +70,7 @@ def golovin(v,t,n0,v0,b): kappa = 1e-10 soluble_fraction = 1. -opts_init.dry_distros = {(kappa, soluble_fraction):expvolumelnr} +opts_init.dry_distros = {(kappa, soluble_fraction, pow(2,14), 0):expvolumelnr} opts_init.kernel = lgrngn.kernel_t.golovin opts_init.terminal_velocity = lgrngn.vt_t.beard77 @@ -112,12 +112,11 @@ def calc_golovin(res,t,n0,v0,b): for i in range(0,2): #loop to test sd_conc and const_multi options for opts_dt in [-1, simulation_time/2.]: # -1 means opts_init.dt (==simulation_time) is used if(i==0): - opts_init.sd_conc = pow(2,14) opts_init.n_sd_max = pow(2,14) else: - opts_init.sd_conc = 0 - opts_init.sd_const_multi = 1000 - opts_init.n_sd_max = int(float(n_zero) / opts_init.sd_const_multi + 10) + sd_const_multi = 1000 + opts_init.dry_distros = {(kappa, soluble_fraction, 0, sd_const_multi):expvolumelnr} + opts_init.n_sd_max = int(float(n_zero) / sd_const_multi + 10) opts.dt = opts_dt if opts_dt < 0: diff --git a/tests/python/physics/coalescence_hall_davis_no_waals.py b/tests/python/physics/coalescence_hall_davis_no_waals.py index 2c896ef9..f4d16478 100644 --- a/tests/python/physics/coalescence_hall_davis_no_waals.py +++ b/tests/python/physics/coalescence_hall_davis_no_waals.py @@ -54,10 +54,9 @@ def expvolumelnr(lnr): kappa = 0 #1e-10 soluble_fraction = 1. -opts_init.dry_distros = {(kappa, soluble_fraction):expvolumelnr} +opts_init.dry_distros = {(kappa, soluble_fraction, pow(2,14), 0):expvolumelnr} -opts_init.sd_conc = pow(2,14) -opts_init.n_sd_max = opts_init.sd_conc +opts_init.n_sd_max = pow(2,14) opts_init.kernel = lgrngn.kernel_t.hall_davis_no_waals opts = lgrngn.opts_t() diff --git a/tests/python/physics/coalescence_onishi_hall.py b/tests/python/physics/coalescence_onishi_hall.py index 4f447e47..0d013068 100755 --- a/tests/python/physics/coalescence_onishi_hall.py +++ b/tests/python/physics/coalescence_onishi_hall.py @@ -34,7 +34,6 @@ def expvolumelnr(lnr): opts_init.dt = 1. opts_init.terminal_velocity = lgrngn.vt_t.beard77fast -opts_init.sd_conc = 1024 opts_init.n_sd_max = 1024 t10_arr = np.zeros((2, n_runs)) @@ -45,7 +44,7 @@ def expvolumelnr(lnr): rhod = 1.22419 * np.ones((1,)) diss_rate = epsilon * np.ones((1,)) -opts_init.dry_distros = {0.:expvolumelnr} +opts_init.dry_distros = {(0., 1024, 0):expvolumelnr} Opts = lgrngn.opts_t() Opts.adve = False diff --git a/tests/python/physics/lgrngn_cond.py b/tests/python/physics/lgrngn_cond.py index 040c6f98..72ccd101 100644 --- a/tests/python/physics/lgrngn_cond.py +++ b/tests/python/physics/lgrngn_cond.py @@ -31,14 +31,13 @@ def lognormal(lnr): opts_init = lgrngn.opts_init_t() kappa = .61 soluble_fraction = 1. -opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} +opts_init.dry_distros = {(kappa, soluble_fraction, int(1e2), 0):lognormal} opts_init.coal_switch = False opts_init.sedi_switch = False opts_init.ice_switch = False opts_init.RH_max = 0.999 # to comply with the assert(RH<1) at init opts_init.dt = 1 -opts_init.sd_conc = int(1e2) -opts_init.n_sd_max = opts_init.sd_conc +opts_init.n_sd_max = int(1e2) backend = lgrngn.backend_t.serial diff --git a/tests/python/physics/lgrngn_cond_substepping.py b/tests/python/physics/lgrngn_cond_substepping.py index 5c0f745d..d1d91646 100644 --- a/tests/python/physics/lgrngn_cond_substepping.py +++ b/tests/python/physics/lgrngn_cond_substepping.py @@ -49,15 +49,14 @@ def lognormal2(lnr): opts = lgrngn.opts_t() opts_init = lgrngn.opts_init_t() -opts_init.dry_distros = {(.61,1.):lognormal, (1.28,1.):lognormal2} # normal mode + GCCNs +opts_init.dry_distros = {(.61,1., 500, 0):lognormal, (1.28,1., 500, 0):lognormal2} # normal mode + GCCNs opts_init.coal_switch = False opts_init.sedi_switch = False # opts_init.RH_max = 1.0001 opts_init.RH_max = 0.95 opts_init.ice_switch = False opts_init.dt = 1 -opts_init.sd_conc = int(1e3) -opts_init.n_sd_max = opts_init.sd_conc +opts_init.n_sd_max = int(1e3) opts_init.rc2_T = 10 # results are the same for 0C to 100C opts_init.sstp_cond_adapt_drw2_eps = 1e-3 #1e-3 @@ -302,7 +301,7 @@ def _write_csv(path, rows): # save results to a CSV file for refdata comparison and for plotting for r in records: - r['sd_conc'] = opts_init.sd_conc + r['sd_conc'] = 1000 r['RH_max'] = opts_init.RH_max r['dt'] = 1 diff --git a/tests/python/physics/puddle.py b/tests/python/physics/puddle.py index 7958b1d5..ed225a6c 100644 --- a/tests/python/physics/puddle.py +++ b/tests/python/physics/puddle.py @@ -18,7 +18,7 @@ def lognormal(lnr): Opts_init = lgrngn.opts_init_t() kappa = .61 soluble_fraction = 1. -Opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} +Opts_init.dry_distros = {(kappa, soluble_fraction, 10000, 0):lognormal} Opts_init.coal_switch = False Opts_init.sedi_switch = True Opts_init.ice_switch = False @@ -34,8 +34,7 @@ def lognormal(lnr): Opts_init.x1 = Opts_init.nx * Opts_init.dx Opts_init.rng_seed = int(time()) -Opts_init.sd_conc = 10000 -Opts_init.n_sd_max = Opts_init.sd_conc * (Opts_init.nx * Opts_init.nz) +Opts_init.n_sd_max = 10000 * (Opts_init.nx * Opts_init.nz) Opts_init.dev_count = 2 diff --git a/tests/python/physics/test_coal.py b/tests/python/physics/test_coal.py index a875dcea..4ef6ffc5 100644 --- a/tests/python/physics/test_coal.py +++ b/tests/python/physics/test_coal.py @@ -31,9 +31,8 @@ def expvolumelnr(lnr): kappa2 = 0.9 soluble_fraction = 1. -opts_init.dry_distros = {(kappa1, soluble_fraction):expvolumelnr, (kappa2, soluble_fraction):expvolumelnr} +opts_init.dry_distros = {(kappa1, soluble_fraction, pow(2,13), 0):expvolumelnr, (kappa2, soluble_fraction, pow(2,13), 0):expvolumelnr} -opts_init.sd_conc = pow(2,14) opts_init.n_sd_max = pow(2,14) opts_init.kernel = lgrngn.kernel_t.geometric diff --git a/tests/python/unit/SD_removal.py b/tests/python/unit/SD_removal.py index 67f75ff8..f66f4fe0 100755 --- a/tests/python/unit/SD_removal.py +++ b/tests/python/unit/SD_removal.py @@ -39,9 +39,8 @@ def expvolumelnr(lnr): kappa = .01 soluble_fraction = 1. -opts_init.dry_distros = {(kappa, soluble_fraction):expvolumelnr} +opts_init.dry_distros = {(kappa, soluble_fraction, 64, 0):expvolumelnr} -opts_init.sd_conc = 64 opts_init.n_sd_max = 64 opts_init.chem_switch = True opts_init.sedi_switch = False diff --git a/tests/python/unit/adve_scheme.py b/tests/python/unit/adve_scheme.py index df310a37..67c3b72c 100644 --- a/tests/python/unit/adve_scheme.py +++ b/tests/python/unit/adve_scheme.py @@ -23,9 +23,8 @@ def lognormal(lnr): kappa = .61 soluble_fraction = 1. -opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} +opts_init.dry_distros = {(kappa, soluble_fraction, 50, 0):lognormal} -opts_init.sd_conc = 50 opts_init.n_sd_max = 50 opts_init.coal_switch = False opts_init.sedi_switch = False diff --git a/tests/python/unit/api_lgrngn.py b/tests/python/unit/api_lgrngn.py index 2f3d865f..11dc49ba 100644 --- a/tests/python/unit/api_lgrngn.py +++ b/tests/python/unit/api_lgrngn.py @@ -26,12 +26,12 @@ def lognormal(lnr): kappa3 = 0.8 soluble_fraction = 0.5 rho_stp = 1.2248 -opts_init.dry_distros = {(kappa1, soluble_fraction): lognormal} +sd_conc = 64 +opts_init.dry_distros = {(kappa1, soluble_fraction, sd_conc, 0):lognormal} opts_init.kernel = lgrngn.kernel_t.geometric opts_init.terminal_velocity = lgrngn.vt_t.beard76 opts_init.adve_scheme = lgrngn.as_t.euler opts_init.dt = 1 -opts_init.sd_conc = 64 opts_init.n_sd_max = int(1e6) # some space for tail SDs opts_init.rng_seed = 396 opts_init.rng_seed_init = 456 @@ -53,12 +53,12 @@ def lognormal(lnr): print("y1 = ", opts_init.y1) print("z1 = ", opts_init.z1) -print("sd_conc = ", opts_init.sd_conc) +print("sd_conc = ", sd_conc) print("RH_max = ", opts_init.RH_max) print("rng_seed = ", opts_init.rng_seed) print("rng_seed_init = ", opts_init.rng_seed_init) print("kernel =", opts_init.kernel) -print("sd_conc =", opts_init.sd_conc) +print("sd_conc =", sd_conc) print("terminal_velocity =", opts_init.terminal_velocity) print("adve_scheme =", opts_init.adve_scheme) print("chem_rho =", opts_init.chem_rho) @@ -106,7 +106,7 @@ def lognormal(lnr): print("chem_gas[SO2] = ", opts.chem_gas[lgrngn.chem_species_t.SO2]) print("chem_gas = ", opts.chem_gas) -opts.src_dry_distros = {(kappa1, soluble_fraction):(lognormal, 64, 1)} +opts.src_dry_distros = {(kappa1, soluble_fraction, 64, 0, 1):lognormal} # --------- test runs ----------- @@ -147,7 +147,7 @@ def lognormal(lnr): prtcls.diag_all() prtcls.diag_sd_conc() print(frombuffer(prtcls.outbuf())) -assert frombuffer(prtcls.outbuf()) == opts_init.sd_conc # parcel set-up +assert frombuffer(prtcls.outbuf()) == sd_conc # parcel set-up # ---------- # 0D (parcel) with explicit calls to sync_in and step_cond @@ -184,7 +184,7 @@ def lognormal(lnr): prtcls.diag_all() prtcls.diag_sd_conc() print(frombuffer(prtcls.outbuf())) -assert frombuffer(prtcls.outbuf()) == opts_init.sd_conc # parcel set-up +assert frombuffer(prtcls.outbuf()) == sd_conc # parcel set-up # ---------- @@ -203,17 +203,16 @@ def lognormal(lnr): assert len(frombuffer(prtcls.outbuf())) == 1 print(frombuffer(prtcls.outbuf())) assert (frombuffer(prtcls.outbuf()) > 0).all() -assert sum(frombuffer(prtcls.outbuf())) >= opts_init.sd_conc +assert sum(frombuffer(prtcls.outbuf())) >= sd_conc # ---------- # 0D const multi - number of SDs and number of particles print("0D const multi") -sd_conc_old = opts_init.sd_conc -opts_init.sd_conc = 0 prtcls_per_cell = 2 * n_tot / rho_stp #rhod=1; 2* because of two distributions -opts_init.sd_const_multi = int(prtcls_per_cell / 64) +sd_const_multi = int(prtcls_per_cell / 64) +opts_init.dry_distros = {(kappa1, soluble_fraction, 0, sd_const_multi):lognormal} prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) prtcls.diag_all() @@ -233,9 +232,8 @@ def lognormal(lnr): prtcls.diag_all() prtcls.diag_wet_mom(0) prtcls_tot = frombuffer(prtcls.outbuf()).sum() -assert ((prtcls_tot / sd_tot) == opts_init.sd_const_multi) -opts_init.sd_const_multi = 0 -opts_init.sd_conc = sd_conc_old +assert ((prtcls_tot / sd_tot) == sd_const_multi) +opts_init.dry_distros = {(kappa1, soluble_fraction, sd_conc, 0):lognormal} @@ -246,8 +244,6 @@ def lognormal(lnr): opts_init.dry_sizes = {(kappa1, soluble_fraction) : {1.e-6 : [30. * rho_stp, 15], 15.e-6 : [10. * rho_stp, 10]}, (kappa2, soluble_fraction) : {1.2e-6 : [20. * rho_stp, 10], 12.e-6 : [15. * rho_stp, 15]}} -sd_conc_old = opts_init.sd_conc -opts_init.sd_conc = 0 prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) @@ -278,9 +274,8 @@ def lognormal(lnr): # go back to distros init -opts_init.sd_conc = sd_conc_old opts_init.dry_sizes = dict() -opts_init.dry_distros = {(kappa1, soluble_fraction):lognormal, (kappa2, soluble_fraction):lognormal} +opts_init.dry_distros = {(kappa1, soluble_fraction, 64, 0):lognormal, (kappa2, soluble_fraction, 64, 0):lognormal} @@ -295,7 +290,7 @@ def lognormal(lnr): prtcls.diag_all() prtcls.diag_sd_conc() print(frombuffer(prtcls.outbuf())) -assert frombuffer(prtcls.outbuf())[0] == 84 # 64 from dry_distro and 20 from sizes +assert frombuffer(prtcls.outbuf())[0] == 148 # 128 from dry_distro and 20 from sizes # go back to distros init opts_init.dry_sizes = dict() @@ -314,7 +309,7 @@ def lognormal(lnr): prtcls.diag_all() prtcls.diag_sd_conc() print(frombuffer(prtcls.outbuf())) -assert frombuffer(prtcls.outbuf())[0] > 84 # 64 from dry_distro and 20 from sizes + tail +assert frombuffer(prtcls.outbuf())[0] > 148 # 128 from dry_distro and 20 from sizes + tail # go back to distros init opts_init.sd_conc_large_tail = 0 @@ -326,9 +321,10 @@ def lognormal(lnr): # 0D dry_sizes + const_multi init print("0D dry_sizes + const_multi") opts_init.dry_sizes = {(kappa3, soluble_fraction) : {1.e-6 : [30. * rho_stp, 15], 15.e-6 : [10. * rho_stp, 5]}} -opts_init.sd_conc = 0 +opts_init.dry_distros = {(kappa1, soluble_fraction, 0, 0):lognormal, (kappa2, soluble_fraction, 0, 0):lognormal} prtcls_per_cell = 2 * n_tot / rho_stp #rhod=1; 2* because of two distributions -opts_init.sd_const_multi = int(prtcls_per_cell / 64) +sd_const_multi = int(prtcls_per_cell / 64) +opts_init.dry_distros = {(kappa1, soluble_fraction, 0, sd_const_multi):lognormal, (kappa2, soluble_fraction, 0, sd_const_multi):lognormal} prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) @@ -339,8 +335,7 @@ def lognormal(lnr): assert frombuffer(prtcls.outbuf())[0] == 84 # 64 from dry_distro and 20 from sizes # go back to distros init -opts_init.sd_conc = sd_conc_old -opts_init.sd_const_multi = 0 +opts_init.dry_distros = {(kappa1, soluble_fraction, 32, 0):lognormal, (kappa2, soluble_fraction, 32, 0):lognormal} opts_init.dry_sizes = dict() @@ -426,7 +421,7 @@ def lognormal(lnr): assert len(frombuffer(prtcls.outbuf())) == opts_init.nx print(frombuffer(prtcls.outbuf())) assert (frombuffer(prtcls.outbuf()) > 0).all() - assert sum(frombuffer(prtcls.outbuf())) == opts_init.nx * opts_init.sd_conc + assert sum(frombuffer(prtcls.outbuf())) == opts_init.nx * sd_conc # ---------- # 2D (periodic horizontal domain) @@ -476,7 +471,7 @@ def lognormal(lnr): assert len(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx print(frombuffer(prtcls.outbuf())) assert (frombuffer(prtcls.outbuf()) > 0).all() - assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.sd_conc + assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * sd_conc assert opts_init.nx == prtcls.opts_init.nx #TODO: test profile vs. 2D array @@ -497,7 +492,7 @@ def lognormal(lnr): assert len(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx print(frombuffer(prtcls.outbuf())) assert (frombuffer(prtcls.outbuf()) > 0).all() -assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.sd_conc +assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * sd_conc assert opts_init.nx == prtcls.opts_init.nx opts_init.turb_adve_switch=False @@ -539,7 +534,7 @@ def lognormal(lnr): assert len(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny print(frombuffer(prtcls.outbuf())) assert (frombuffer(prtcls.outbuf()) > 0).all() - assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny * opts_init.sd_conc + assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny * sd_conc print("3D turb") eps = arr_t([eps, eps]) @@ -557,7 +552,7 @@ def lognormal(lnr): assert len(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny print(frombuffer(prtcls.outbuf())) assert (frombuffer(prtcls.outbuf()) > 0).all() -assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny * opts_init.sd_conc +assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny * sd_conc opts_init.turb_adve_switch=False opts.turb_adve=False @@ -579,17 +574,18 @@ def lognormal(lnr): assert len(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny print(frombuffer(prtcls.outbuf())) assert (frombuffer(prtcls.outbuf()) > 0).all() -assert sum(frombuffer(prtcls.outbuf())) >= opts_init.nz * opts_init.nx * opts_init.ny * opts_init.sd_conc +assert sum(frombuffer(prtcls.outbuf())) >= opts_init.nz * opts_init.nx * opts_init.ny * sd_conc # ---------- # 3D const multi - number of SDs and number of particles print("3D const multi") -opts_init.sd_conc = 0 +opts_init.dry_distros = {(kappa1, soluble_fraction, 0, 0):lognormal, (kappa2, soluble_fraction, 0, 0):lognormal} cell_vol = opts_init.dx * opts_init.dy * opts_init.dz prtcls_per_cell = 2 * n_tot * cell_vol / rho_stp #rhod=1; 2* because of two distributions -opts_init.sd_const_multi = int(prtcls_per_cell / 64) +sd_const_multi = int(prtcls_per_cell / 64) +opts_init.dry_distros = {(kappa1, soluble_fraction, 0, sd_const_multi):lognormal, (kappa2, soluble_fraction, 0, sd_const_multi):lognormal} n_cell = opts_init.nz * opts_init.nx * opts_init.ny prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) @@ -611,7 +607,7 @@ def lognormal(lnr): prtcls.diag_all() prtcls.diag_wet_mom(0) prtcls_tot = frombuffer(prtcls.outbuf()).sum() -assert ((prtcls_tot / sd_tot) * cell_vol == opts_init.sd_const_multi) +assert ((prtcls_tot / sd_tot) * cell_vol == sd_const_multi) @@ -652,10 +648,8 @@ def lognormal(lnr): # 3D dry_sizes + sd_conc init print("3D dry_sizes + sd_conc") soluble_fraction = 1 # no insoluble aerosol from now on -opts_init.dry_distros = {(kappa1, soluble_fraction):lognormal, (kappa2, soluble_fraction):lognormal} +opts_init.dry_distros = {(kappa1, soluble_fraction, 64, 0):lognormal, (kappa2, soluble_fraction, 64, 0):lognormal} opts_init.dry_sizes = {(kappa1, soluble_fraction) : {1.e-6 : [30./ cell_vol * rho_stp, 15], 15.e-6 : [10. / cell_vol * rho_stp, 5]}} -opts_init.sd_conc = sd_conc_old -opts_init.sd_const_multi = 0 prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) @@ -663,13 +657,13 @@ def lognormal(lnr): prtcls.diag_all() prtcls.diag_sd_conc() print(frombuffer(prtcls.outbuf())) -assert (frombuffer(prtcls.outbuf()) == 84).all() # 64 from dry_distro and 20 from sizes +assert (frombuffer(prtcls.outbuf()) == 148).all() # 128 from dry_distro and 20 from sizes # test if get_attr work and if kappas are set correctly kappa = asarray(prtcls.get_attr("kappa")) # assert (kappa[:(32*opts_init.nx*opts_init.ny*opts_init.nz)] == kappa2).all() # assert (kappa[(32*opts_init.nx*opts_init.ny*opts_init.nz):] == kappa1).all() -n = 32 * opts_init.nx * opts_init.ny * opts_init.nz +n = 64 * opts_init.nx * opts_init.ny * opts_init.nz assert (kappa[:n] == kappa1).all() assert (kappa[n:2*n] == kappa2).all() assert (kappa[2*n:] == kappa1).all() @@ -687,7 +681,7 @@ def lognormal(lnr): prtcls.diag_all() prtcls.diag_sd_conc() print(frombuffer(prtcls.outbuf())) -assert (frombuffer(prtcls.outbuf())[0] > 64 + 20).all() # 64 from dry_distro and 20 from sizes + tail +assert (frombuffer(prtcls.outbuf())[0] > 128 + 20).all() # 128 from dry_distro and 20 from sizes + tail # go back to distros init @@ -700,9 +694,10 @@ def lognormal(lnr): # 3D dry_sizes + const_multi init print("3D dry_sizes + const_multi") opts_init.dry_sizes = {(kappa1, soluble_fraction) : {1.e-6 : [30./ cell_vol * rho_stp, 15], 15.e-6 : [10. / cell_vol * rho_stp, 5]}} -opts_init.sd_conc = 0 +opts_init.dry_distros = {(kappa1, soluble_fraction, 0, 0):lognormal, (kappa2, soluble_fraction, 0, 0):lognormal} prtcls_per_cell = 2 * n_tot * cell_vol / rho_stp #rhod=1; 2* because of two distributions -opts_init.sd_const_multi = int(prtcls_per_cell / 64) +sd_const_multi = int(prtcls_per_cell / 64) +opts_init.dry_distros = {(kappa1, soluble_fraction, 0, sd_const_multi):lognormal, (kappa2, soluble_fraction, 0, sd_const_multi):lognormal} prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) @@ -713,8 +708,7 @@ def lognormal(lnr): assert (frombuffer(prtcls.outbuf())[0] == 84).all() # 64 from dry_distro and 20 from sizes # go back to distros init -opts_init.sd_conc = sd_conc_old -opts_init.sd_const_multi = 0 +opts_init.dry_distros = {(kappa1, soluble_fraction, 32, 0):lognormal, (kappa2, soluble_fraction, 32, 0):lognormal} opts_init.dry_sizes = dict() @@ -736,7 +730,7 @@ def lognormal(lnr): assert len(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny print(frombuffer(prtcls.outbuf())) assert (frombuffer(prtcls.outbuf()) > 0).all() -assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny * opts_init.sd_conc +assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny * sd_conc opts_init.ice_switch = False opts_init.coal_switch = True diff --git a/tests/python/unit/chem_coal.py b/tests/python/unit/chem_coal.py index 8ee9a814..da565a56 100755 --- a/tests/python/unit/chem_coal.py +++ b/tests/python/unit/chem_coal.py @@ -36,9 +36,8 @@ def expvolumelnr(lnr): kappa = 0.1 soluble_fraction = 1. -opts_init.dry_distros = {(kappa, soluble_fraction):expvolumelnr} +opts_init.dry_distros = {(kappa, soluble_fraction, 64, 0):expvolumelnr} -opts_init.sd_conc = 64 opts_init.n_sd_max = 64 opts_init.chem_switch = True diff --git a/tests/python/unit/col_kernels.py b/tests/python/unit/col_kernels.py index 4d418d10..9a7ce556 100644 --- a/tests/python/unit/col_kernels.py +++ b/tests/python/unit/col_kernels.py @@ -31,8 +31,7 @@ def lognormal(lnr): print(kernel) opts_init = lgrngn.opts_init_t() opts_init.dt = 1 - opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} - opts_init.sd_conc = 50 + opts_init.dry_distros = {(kappa, soluble_fraction, 50, 0):lognormal} opts_init.n_sd_max = 50 opts_init.terminal_velocity=lgrngn.vt_t.beard76 opts_init.kernel = kernel diff --git a/tests/python/unit/diag_incloud_time.py b/tests/python/unit/diag_incloud_time.py index 50a27467..b122a9a6 100644 --- a/tests/python/unit/diag_incloud_time.py +++ b/tests/python/unit/diag_incloud_time.py @@ -20,13 +20,12 @@ def lognormal(lnr): opts = lgrngn.opts_t() opts_init = lgrngn.opts_init_t() -opts_init.dry_distros = {(.61, 1.):lognormal, (1.28, 1.):lognormal} +opts_init.dry_distros = {(.61, 1., 50, 0):lognormal, (1.28, 1., 50, 0):lognormal} opts_init.coal_switch = False opts_init.sedi_switch = False opts_init.RH_max = 0.999 # to comply with the assert(RH<1) at init opts_init.dt = 0.1 -opts_init.sd_conc = int(1e2) -opts_init.n_sd_max = opts_init.sd_conc +opts_init.n_sd_max = int(1e2) opts_init.diag_incloud_time = True backend = lgrngn.backend_t.serial diff --git a/tests/python/unit/ice_SD.py b/tests/python/unit/ice_SD.py index 5ffd7045..03c7f258 100644 --- a/tests/python/unit/ice_SD.py +++ b/tests/python/unit/ice_SD.py @@ -21,9 +21,8 @@ def lognormal(lnr): kappa = .61 soluble_fraction = 0.7 -opts_init.dry_distros = {(kappa, soluble_fraction): lognormal} +opts_init.dry_distros = {(kappa, soluble_fraction, 100, 0):lognormal} opts_init.dt = 0.1 -opts_init.sd_conc = 100 opts_init.n_sd_max = 1000 opts_init.RH_max = 0.95 diff --git a/tests/python/unit/lgrngn_adve.py b/tests/python/unit/lgrngn_adve.py index 5b975f93..c6de6afe 100644 --- a/tests/python/unit/lgrngn_adve.py +++ b/tests/python/unit/lgrngn_adve.py @@ -25,7 +25,7 @@ def lognormal(lnr): Opts_init = lgrngn.opts_init_t() kappa = .61 soluble_fraction = 1. -Opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} +Opts_init.dry_distros = {(kappa, soluble_fraction, 10, 0):lognormal} Opts_init.coal_switch = False Opts_init.sedi_switch = False @@ -38,7 +38,6 @@ def lognormal(lnr): Opts_init.z1 = Opts_init.nz * Opts_init.dz Opts_init.x1 = Opts_init.nx * Opts_init.dx -Opts_init.sd_conc = 10 Opts_init.n_sd_max = 10 * (Opts_init.nx * Opts_init.nz) Backend = lgrngn.backend_t.serial diff --git a/tests/python/unit/lgrngn_subsidence.py b/tests/python/unit/lgrngn_subsidence.py index 4b5b516e..c35cbe28 100644 --- a/tests/python/unit/lgrngn_subsidence.py +++ b/tests/python/unit/lgrngn_subsidence.py @@ -18,7 +18,7 @@ def lognormal(lnr): Opts_init = lgrngn.opts_init_t() kappa = .61 soluble_fraction = 1. -Opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} +Opts_init.dry_distros = {(kappa, soluble_fraction, 1000, 0):lognormal} Opts_init.coal_switch = False Opts_init.sedi_switch = False Opts_init.subs_switch = True @@ -35,8 +35,7 @@ def lognormal(lnr): Opts_init.x1 = Opts_init.nx * Opts_init.dx Opts_init.rng_seed = int(time()) -Opts_init.sd_conc = 1000 -Opts_init.n_sd_max = Opts_init.sd_conc * (Opts_init.nx * Opts_init.nz) +Opts_init.n_sd_max = 1000 * (Opts_init.nx * Opts_init.nz) Opts_init.w_LS = np.array([0.,1.,2.,0.,1.,2.]) # 1/s large-scale subsidence Backend = lgrngn.backend_t.serial diff --git a/tests/python/unit/lgrngn_turb_adve.py b/tests/python/unit/lgrngn_turb_adve.py index 49d821e7..ea428cf7 100644 --- a/tests/python/unit/lgrngn_turb_adve.py +++ b/tests/python/unit/lgrngn_turb_adve.py @@ -17,7 +17,7 @@ def lognormal(lnr): Opts_init = lgrngn.opts_init_t() kappa = .61 -Opts_init.dry_distros = {kappa:lognormal} +Opts_init.dry_distros = {(kappa, 1., 100, 0):lognormal} Opts_init.coal_switch = False Opts_init.sedi_switch = False Opts_init.turb_adve_switch = True @@ -32,8 +32,7 @@ def lognormal(lnr): Opts_init.x1 = Opts_init.nx * Opts_init.dx Opts_init.rng_seed = int(time()) -Opts_init.sd_conc = 100 -Opts_init.n_sd_max = Opts_init.sd_conc * (Opts_init.nx * Opts_init.nz) +Opts_init.n_sd_max = 100 * (Opts_init.nx * Opts_init.nz) Backend = lgrngn.backend_t.serial diff --git a/tests/python/unit/multiple_kappas.py b/tests/python/unit/multiple_kappas.py index 7316c9a5..deb6f0d5 100644 --- a/tests/python/unit/multiple_kappas.py +++ b/tests/python/unit/multiple_kappas.py @@ -35,12 +35,11 @@ def check_kappa_conc(prtcls, eps): kappa2 = 1.28 soluble_fraction = 1. rho_stp = 1.2248 -opts_init.dry_distros = {(kappa1, soluble_fraction):lognormal, (kappa2, soluble_fraction):lognormal} +opts_init.dry_distros = {(kappa1, soluble_fraction, 32, 0):lognormal, (kappa2, soluble_fraction, 64, 0):lognormal} opts_init.kernel = lgrngn.kernel_t.geometric opts_init.terminal_velocity = lgrngn.vt_t.beard76 opts_init.dt = 1 -opts_init.sd_conc = 64 -opts_init.n_sd_max = 512 +opts_init.n_sd_max = 1024 opts_init.rng_seed = 396 opts_init.src_z1 = opts_init.dz opts_init.sedi_switch = False @@ -79,15 +78,16 @@ def check_kappa_conc(prtcls, eps): prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) -check_kappa_conc(prtcls, 5e-3) +check_kappa_conc(prtcls, 7e-3) # 3D const multi - number of SDs and number of particles -opts_init.sd_conc = 0 cell_vol = opts_init.dx * opts_init.dy * opts_init.dz prtcls_per_cell = 2 * n_tot * cell_vol / rho_stp #rhod=1 -opts_init.sd_const_multi = int(prtcls_per_cell / 64) +sd_const_multi1 = int(prtcls_per_cell / 64) +sd_const_multi2 = int(prtcls_per_cell / 32) +opts_init.dry_distros = {(kappa1, soluble_fraction, 0, sd_const_multi1):lognormal, (kappa2, soluble_fraction, 0, sd_const_multi2):lognormal} n_cell = opts_init.nz * opts_init.nx * opts_init.ny -opts_init.n_sd_max = int(n_cell * prtcls_per_cell / opts_init.sd_const_multi) # 2* because of two distributions +opts_init.n_sd_max = int(n_cell * (n_tot * cell_vol / rho_stp / sd_const_multi1 + n_tot * cell_vol / rho_stp / sd_const_multi2)) prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) diff --git a/tests/python/unit/relax.py b/tests/python/unit/relax.py index 59e92954..ec231817 100644 --- a/tests/python/unit/relax.py +++ b/tests/python/unit/relax.py @@ -105,9 +105,8 @@ def test(opts_init): kappa = .61 soluble_fraction = 1. opts_init = lgrngn.opts_init_t() -opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} +opts_init.dry_distros = {(kappa, soluble_fraction, 1024, 0):lognormal} opts_init.rlx_dry_distros = {kappa: [lognormal_rlx, [0,2],[0,opts_init.dz]]} -opts_init.sd_conc = 1024 opts_init.rlx_bins = 1024 opts_init.rlx_timescale = 4 # whole simulation time is 2, so this means we should get half of the droplets added @@ -117,7 +116,7 @@ def test(opts_init): print(' --- dry_distros rlx sd_per_bin = 1 ---') opts_init.rlx_sd_per_bin = 1 -opts_init.n_sd_max = int((opts_init.sd_conc * 2 + opts_init.rlx_bins * opts_init.rlx_sd_per_bin * 2) * 2) # assuming nx=nz=2 +opts_init.n_sd_max = int((1024 * 2 + opts_init.rlx_bins * opts_init.rlx_sd_per_bin * 2) * 2) # assuming nx=nz=2 sd_conc, wet_mom0, wet_mom1 = test(opts_init) @@ -147,7 +146,7 @@ def test(opts_init): print(' --- dry_distros rlx sd_per_bin = 10 ---') opts_init.rlx_sd_per_bin = 10 -opts_init.n_sd_max = int((opts_init.sd_conc * 2 + opts_init.rlx_bins * opts_init.rlx_sd_per_bin * 2) * 2) # assuming nx=nz=2 +opts_init.n_sd_max = int((1024 * 2 + opts_init.rlx_bins * opts_init.rlx_sd_per_bin * 2) * 2) # assuming nx=nz=2 sd_conc, wet_mom0, wet_mom1 = test(opts_init) diff --git a/tests/python/unit/segfault_20150216.py b/tests/python/unit/segfault_20150216.py index b65618be..5222751c 100644 --- a/tests/python/unit/segfault_20150216.py +++ b/tests/python/unit/segfault_20150216.py @@ -26,8 +26,7 @@ def lognormal(lnr): kappa = .61 soluble_fraction = 1. -opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} -opts_init.sd_conc = 50 +opts_init.dry_distros = {(kappa, soluble_fraction, 50, 0):lognormal} opts_init.n_sd_max = 50 try: diff --git a/tests/python/unit/source.py b/tests/python/unit/source.py index af7fc923..76fefaeb 100644 --- a/tests/python/unit/source.py +++ b/tests/python/unit/source.py @@ -97,12 +97,11 @@ def test(opts_init, opts): print(' --- dry_distros simple src ---') opts_init = lgrngn.opts_init_t() opts = lgrngn.opts_t() -opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} -opts_init.sd_conc = 1024 +opts_init.dry_distros = {(kappa, soluble_fraction, 1024, 0):lognormal} src_sd_conc = 512 supstp_src = 50 -opts.src_dry_distros = {(kappa, soluble_fraction):(lognormal_src, src_sd_conc, supstp_src)} -opts_init.n_sd_max = int((opts_init.sd_conc * 2 + src_sd_conc * 2) * 2) # assuming nx=nz=2 +opts.src_dry_distros = {(kappa, soluble_fraction, src_sd_conc, 0, supstp_src):lognormal_src} +opts_init.n_sd_max = int((1024 * 2 + src_sd_conc * 2) * 2) # assuming nx=nz=2 opts_init.src_type = lgrngn.src_t.simple sd_conc, wet_mom0, wet_mom1 = test(opts_init, opts) @@ -121,16 +120,34 @@ def test(opts_init, opts): if (abs( (7.84 / 2.12) - (wet_mom1[0] + wet_mom1[2]) / (wet_mom1[1] + wet_mom1[3]) ) > 0.015): raise Exception("incorrect radius after source") +# ----------- test source with dry_distros simple, constant multiplicity ------------------ +print(' --- dry_distros simple constant multiplicity src ---') +opts_init = lgrngn.opts_init_t() +opts = lgrngn.opts_t() +opts_init.dry_distros = {(kappa, soluble_fraction, 1024, 0):lognormal} +src_const_multi = 100000 +supstp_src = 50 +opts.src_dry_distros = {(kappa, soluble_fraction, 0, src_const_multi, supstp_src):lognormal_src} +opts_init.n_sd_max = int((1024 * 2 + 300 * 2) * 2) # assuming nx=nz=2 +opts_init.src_type = lgrngn.src_t.simple + +sd_conc, wet_mom0, wet_mom1 = test(opts_init, opts) + +print('diag_sd_conc', sd_conc) +if not(sd_conc[0] == 1514 and sd_conc[2] == 1514): + raise Exception("wrong amount of constant-multiplicity SDs were added") +if not(sd_conc[1] == 1024 and sd_conc[3] == 1024): + raise Exception("constant-multiplicity SDs were added in wrong cells") + # --------------- test source with dry_distros matching ------------------ print(' --- dry_distros matching src ---') opts_init = lgrngn.opts_init_t() opts = lgrngn.opts_t() -opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} -opts_init.sd_conc = 1024 +opts_init.dry_distros = {(kappa, soluble_fraction, 1024, 0):lognormal} src_sd_conc = 512 supstp_src = 50 -opts.src_dry_distros = {(kappa, soluble_fraction):(lognormal_src, src_sd_conc, supstp_src)} -opts_init.n_sd_max = int((opts_init.sd_conc * 2 + src_sd_conc * 2) * 2) # assuming nx=nz=2 +opts.src_dry_distros = {(kappa, soluble_fraction, src_sd_conc, 0, supstp_src):lognormal_src} +opts_init.n_sd_max = int((1024 * 2 + src_sd_conc * 2) * 2) # assuming nx=nz=2 opts_init.src_type = lgrngn.src_t.matching sd_conc, wet_mom0, wet_mom1 = test(opts_init, opts) @@ -155,7 +172,7 @@ def test(opts_init, opts): opts = lgrngn.opts_t() supstp_src = 50 opts_init.dry_sizes = {(kappa, soluble_fraction) : {1.e-6 : [30., 20], 15.e-6 : [10., 10]}} -opts.src_dry_sizes = {(kappa, soluble_fraction) : {1.e-6 : [0.3, 10, supstp_src], 15.e-6 : [0.1, 5, supstp_src]}} +opts.src_dry_sizes = {(kappa, soluble_fraction, supstp_src) : {1.e-6 : [0.3, 10], 15.e-6 : [0.1, 5]}} opts_init.n_sd_max=240 opts_init.src_type = lgrngn.src_t.simple # dry sizes works the same for simple and matching (no matching done) diff --git a/tests/python/unit/sstp_cond.py b/tests/python/unit/sstp_cond.py index 60371b8a..77a1cbce 100644 --- a/tests/python/unit/sstp_cond.py +++ b/tests/python/unit/sstp_cond.py @@ -21,12 +21,11 @@ def test(turb_cond): opts_init = lgrngn.opts_init_t() kappa = .61 soluble_fraction = 1. - opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} + opts_init.dry_distros = {(kappa, soluble_fraction, 64, 0):lognormal} opts_init.coal_switch=0 opts_init.sedi_switch=0 opts_init.ice_switch=0 opts_init.dt = 1 - opts_init.sd_conc = 64 opts_init.n_sd_max = 512 opts_init.rng_seed = 396 opts_init.exact_sstp_cond = True # test would fail with per-cell sstp logic diff --git a/tests/python/unit/terminal_velocities.py b/tests/python/unit/terminal_velocities.py index e49b7823..6b2c8f4f 100644 --- a/tests/python/unit/terminal_velocities.py +++ b/tests/python/unit/terminal_velocities.py @@ -28,9 +28,8 @@ def lognormal(lnr): kappa = .61 soluble_fraction = 1. -opts_init.dry_distros = {(kappa, soluble_fraction):lognormal} +opts_init.dry_distros = {(kappa, soluble_fraction, 50, 0):lognormal} -opts_init.sd_conc = 50 opts_init.n_sd_max = 50 opts_init.kernel = lgrngn.kernel_t.geometric diff --git a/tests/python/unit/uniform_init.py b/tests/python/unit/uniform_init.py index f0111dfd..34ef84eb 100644 --- a/tests/python/unit/uniform_init.py +++ b/tests/python/unit/uniform_init.py @@ -44,10 +44,9 @@ def expvolumelnr(lnr): kappa = 1e-6 soluble_fraction = 1. -opts_init.dry_distros = {(kappa, soluble_fraction):expvolumelnr} +opts_init.dry_distros = {(kappa, soluble_fraction, 64, 0):expvolumelnr} -opts_init.sd_conc = 64 -opts_init.n_sd_max = opts_init.sd_conc * opts_init.nx * opts_init.ny * opts_init.nz +opts_init.n_sd_max = 64 * opts_init.nx * opts_init.ny * opts_init.nz try: prtcls = lgrngn.factory(lgrngn.backend_t.OpenMP, opts_init) @@ -67,9 +66,9 @@ def expvolumelnr(lnr): relative difference between water content in one of the cells and mean value greater than 15%: " \ + str(abs(water_content - mean_water_content_sd_conc)/water_content) + " > 0.15") -opts_init.sd_conc = 0 -opts_init.sd_const_multi = 100000 -opts_init.n_sd_max = int(opts_init.nx * opts_init.ny * opts_init.nz * (n_zero / opts_init.sd_const_multi + 100)); #TODO: why do we need to add this 100? integral not correct? +sd_const_multi = 100000 +opts_init.dry_distros = {(kappa, soluble_fraction, 0, sd_const_multi):expvolumelnr} +opts_init.n_sd_max = int(opts_init.nx * opts_init.ny * opts_init.nz * (n_zero / sd_const_multi + 100)); #TODO: why do we need to add this 100? integral not correct? try: prtcls = lgrngn.factory(lgrngn.backend_t.OpenMP, opts_init)