diff --git a/bindings/python/common.hpp b/bindings/python/common.hpp index c361d6f4..7e05e5fc 100644 --- a/bindings/python/common.hpp +++ b/bindings/python/common.hpp @@ -119,20 +119,22 @@ namespace libcloudphxx } template - real_t rw3_cr(const real_t &rd3, const real_t &kappa, const real_t &T) + real_t rw3_cr(const real_t &rd3, const real_t &rd3_insol, const real_t &kappa, const real_t &T) { return cmn::kappa_koehler::rw3_cr( rd3 * si::cubic_metres, + rd3_insol * si::cubic_metres, kappa * si::dimensionless(), T * si::kelvins ) / si::cubic_metres; } template - real_t S_cr(const real_t &rd3, const real_t &kappa, const real_t &T) + real_t S_cr(const real_t &rd3, const real_t &rd3_insol, const real_t &kappa, const real_t &T) { return cmn::kappa_koehler::S_cr( rd3 * si::cubic_metres, + rd3_insol * si::cubic_metres, kappa * si::dimensionless(), T * si::kelvins ); diff --git a/include/libcloudph++/common/ice_nucleation.hpp b/include/libcloudph++/common/ice_nucleation.hpp index 58340c84..1548029f 100644 --- a/include/libcloudph++/common/ice_nucleation.hpp +++ b/include/libcloudph++/common/ice_nucleation.hpp @@ -21,7 +21,7 @@ namespace libcloudphxx BOOST_GPU_ENABLED quantity T_freeze_CDF_inv( const INP_t& INP_type, // type of ice nucleating particle - const real_t rd2_insol, // radius squared of insoluble particle in m^2 + const real_t rd3_insol, // radius cubed of insoluble particle in m^3 const real_t rand // random number between [0, 1] ) { real_t A = real_t(4) @@ -30,7 +30,7 @@ namespace libcloudphxx #else * CUDART_PI #endif - * rd2_insol; // surface area of the insoluble particle + * pow(rd3_insol, real_t(2./3.)); // surface area of the insoluble particle if (INP_type == INP_t::mineral && A > real_t(1e-20)) { @@ -54,12 +54,12 @@ namespace libcloudphxx BOOST_GPU_ENABLED real_t operator()(const thrust::tuple &tpl) const { - const real_t &rd2_insol = thrust::get<0>(tpl); // from rd2 vector + const real_t &rd3_insol = thrust::get<0>(tpl); // from rd3 vector const real_t &rand = thrust::get<1>(tpl); // from rand vector return ice_nucleation::template T_freeze_CDF_inv( INP_type, - rd2_insol, + rd3_insol, rand ) / si::kelvin; } @@ -71,13 +71,13 @@ namespace libcloudphxx BOOST_GPU_ENABLED real_t p_freeze( const INP_t& INP_type, // type of ice nucleating particle - const real_t rd2_insol, // radius squared of insoluble particle in m^2 + const real_t rd3_insol, // radius cubed of insoluble particle in m^3 const real_t rw2, // wet radius squared in m^2 const real_t T, // temperature in kelvin const real_t dt // time step in seconds ) { - if (rd2_insol > real_t(0)) + if (rd3_insol > real_t(0)) { real_t A = real_t(4) #if !defined(__NVCC__) @@ -85,7 +85,7 @@ namespace libcloudphxx #else * CUDART_PI #endif - * rd2_insol; // surface area of the insoluble particle + * pow(rd3_insol, real_t(2./3.)); // surface area of the insoluble particle real_t d_aw = real_t(1) - const_cp::p_vsi(T * si::kelvin)/ const_cp::p_vs(T * si::kelvin); // water activity if (INP_type == INP_t::mineral) { @@ -127,13 +127,13 @@ namespace libcloudphxx BOOST_GPU_ENABLED real_t operator()(const thrust::tuple &tpl) const { - const real_t &rd2_insol = thrust::get<0>(tpl); // radius squared of insoluble particle + const real_t &rd3_insol = thrust::get<0>(tpl); // radius cubed of insoluble particle const real_t &rw2 = thrust::get<1>(tpl); // wet radius squared const real_t &T = thrust::get<2>(tpl); // temperature in kelvin return ice_nucleation::p_freeze( INP_type, - rd2_insol, + rd3_insol, rw2, T, dt diff --git a/include/libcloudph++/common/kappa_koehler.hpp b/include/libcloudph++/common/kappa_koehler.hpp index 4292937c..4896e096 100644 --- a/include/libcloudph++/common/kappa_koehler.hpp +++ b/include/libcloudph++/common/kappa_koehler.hpp @@ -29,15 +29,16 @@ namespace libcloudphxx template BOOST_GPU_ENABLED quantity rw3_eq_nokelvin( - quantity rd3, + quantity rd3, + quantity rd3_insol, quantity kappa, quantity RH ) - { + { assert(RH > .03); // kappa-koehler assumes well dissolved matter assert(RH < 1.); // no equilibrium over RH=100% - assert(kappa >= 0); - return rd3 * (1 - RH * (1 - kappa)) / (1 - RH); + assert(kappa >= 0); + return (rd3 * (1 - RH * (1 - kappa)) + rd3_insol) / (1 - RH); } /// @brief activity of water in solution (eqs. 1,6) in @copydetails Petters_and_Kreidenweis_2007 @@ -46,11 +47,13 @@ namespace libcloudphxx quantity a_w( quantity rw3, quantity rd3, + quantity rd3_insol, quantity kappa ) { assert(kappa > 0); - return (rw3 - rd3) / (rw3 - rd3 * (real_t(1) - kappa)); + assert(rw3 >= rd3 + rd3_insol - std::numeric_limits::epsilon() * si::cubic_metres); + return (rw3 - rd3 - rd3_insol) / (rw3 - rd3 * (real_t(1) - kappa) - rd3_insol); } // functor to be passed to the minimisation functions used below @@ -62,6 +65,7 @@ namespace libcloudphxx { const quantity RH; const quantity rd3; + const quantity rd3_insol; const quantity kappa; const quantity T; @@ -70,9 +74,10 @@ namespace libcloudphxx rw3_eq_minfun( const quantity &RH, const quantity &rd3, + const quantity &rd3_insol, const quantity &kappa, const quantity &T - ) : RH(RH), rd3(rd3), kappa(kappa), T(T) {} + ) : RH(RH), rd3(rd3), rd3_insol(rd3_insol), kappa(kappa), T(T) {} BOOST_GPU_ENABLED real_t operator()(real_t rw3) const @@ -80,10 +85,15 @@ namespace libcloudphxx #if !defined(__NVCC__) using std::pow; using std::cbrt; + using std::max; +#else + using thrust::max; #endif + + real_t rw3_safe = max(rw3, real_t((this->rd3 + this->rd3_insol) / si::cubic_metres)); return this->RH - - a_w(rw3 * si::cubic_metres, this->rd3, this->kappa) - * kelvin::klvntrm(real_t(cbrt(rw3)) * si::metres, this->T); + - a_w(rw3_safe * si::cubic_metres, this->rd3, this->rd3_insol, this->kappa) + * kelvin::klvntrm(real_t(cbrt(rw3_safe)) * si::metres, this->T); } }; @@ -92,6 +102,7 @@ namespace libcloudphxx struct rw3_cr_minfun { const quantity rd3; + const quantity rd3_insol; const quantity kappa; const quantity T; @@ -99,9 +110,10 @@ namespace libcloudphxx BOOST_GPU_ENABLED rw3_cr_minfun( const quantity &rd3, + const quantity &rd3_insol, const quantity &kappa, const quantity &T - ) : rd3(rd3), kappa(kappa), T(T) {} + ) : rd3(rd3), rd3_insol(rd3_insol), kappa(kappa), T(T) {} BOOST_GPU_ENABLED real_t operator()(real_t rw3) const @@ -111,8 +123,8 @@ namespace libcloudphxx using std::cbrt; #endif return (kelvin::A(T) - * (rd3 - rw3 * si::cubic_metres) - * ((kappa - 1) * rd3 + rw3 * si::cubic_metres) + * (rd3 + rd3_insol - rw3 * si::cubic_metres) + * ((kappa - 1) * rd3 - rd3_insol + rw3 * si::cubic_metres) + 3 * kappa * rd3 * rw3 * cbrt(rw3) * si::cubic_metres * si::metres ) / si::cubic_metres / si::cubic_metres / si::metres; } @@ -127,7 +139,8 @@ namespace libcloudphxx template BOOST_GPU_ENABLED quantity rw3_eq( - quantity rd3, + quantity rd3, + quantity rd3_insol, quantity kappa, quantity RH, quantity T @@ -136,12 +149,12 @@ namespace libcloudphxx assert(RH < 1); // no equilibrium over RH=100% assert(kappa >= 0); - if(kappa == 0) return rd3; + if(kappa == 0) return rd3 + rd3_insol; return common::detail::toms748_solve( - detail::rw3_eq_minfun(RH, rd3, kappa, T), // the above-defined functor - real_t(rd3 / si::cubic_metres), // min - real_t(rw3_eq_nokelvin(rd3, kappa, RH) / si::cubic_metres) // max + detail::rw3_eq_minfun(RH, rd3, rd3_insol, kappa, T), // the above-defined functor + real_t((rd3 + rd3_insol) / si::cubic_metres), // min + real_t(rw3_eq_nokelvin(rd3, rd3_insol, kappa, RH) / si::cubic_metres) // max ) * si::cubic_metres; } @@ -151,7 +164,8 @@ namespace libcloudphxx template BOOST_GPU_ENABLED quantity rw3_cr( - quantity rd3, + quantity rd3, + quantity rd3_insol, quantity kappa, quantity T ) @@ -159,12 +173,13 @@ namespace libcloudphxx assert(kappa > 0); quantity rd3_dbl = static_cast >(rd3); + quantity rd3_insol_dbl = static_cast >(rd3_insol); quantity T_dbl = static_cast >(T); return real_t(common::detail::toms748_solve( - detail::rw3_cr_minfun(rd3_dbl, double(kappa), T_dbl), // the above-defined functor - double(1e0 * (rd3 / si::cubic_metres)), // min - double(1e8 * (rd3 / si::cubic_metres)) // max + detail::rw3_cr_minfun(rd3_dbl, rd3_insol_dbl, double(kappa), T_dbl), // the above-defined functor + double(1e0 * ((rd3 + rd3_insol) / si::cubic_metres)), // min + double(1e8 * ((rd3 + rd3_insol) / si::cubic_metres)) // max )) * si::cubic_metres; } @@ -173,6 +188,7 @@ namespace libcloudphxx BOOST_GPU_ENABLED quantity S_cr( quantity rd3, + quantity rd3_insol, quantity kappa, quantity T ) @@ -182,9 +198,9 @@ namespace libcloudphxx using std::pow; using std::cbrt; #endif - quantity rw3 = rw3_cr(rd3, kappa, T); + quantity rw3 = rw3_cr(rd3, rd3_insol, kappa, T); - return a_w(rw3, rd3, kappa) * kelvin::klvntrm( + return a_w(rw3, rd3, rd3_insol, kappa) * kelvin::klvntrm( real_t(cbrt(real_t(rw3 / si::cubic_metres))) * si::metres, T ); diff --git a/src/impl/condensation/common/particles_impl_cond_common.ipp b/src/impl/condensation/common/particles_impl_cond_common.ipp index aa9da4df..f33d3d88 100644 --- a/src/impl/condensation/common/particles_impl_cond_common.ipp +++ b/src/impl/condensation/common/particles_impl_cond_common.ipp @@ -93,13 +93,14 @@ namespace libcloudphxx const quantity RH_max; const quantity lambda_D; const quantity lambda_K; + const quantity rd3_insol; // ctor BOOST_GPU_ENABLED advance_rw2_minfun( const real_t &dt, const real_t &rw2, - const thrust::tuple, real_t, real_t> &tpl, + const thrust::tuple, real_t, real_t> &tpl, const real_t &RH_max ) : dt(dt * si::seconds), @@ -115,6 +116,7 @@ namespace libcloudphxx RH( thrust::get<2>(tpl)), lambda_D(thrust::get<7>(thrust::get<0>(tpl)) * si::metres), lambda_K(thrust::get<8>(thrust::get<0>(tpl)) * si::metres), + rd3_insol(thrust::get<9>(thrust::get<0>(tpl)) * si::cubic_metres), RH_max(RH_max) {} @@ -157,7 +159,7 @@ namespace libcloudphxx T, p, RH > RH_max ? RH_max : RH, - a_w(rw3, rd3, kpa), + a_w(rw3, rd3, rd3_insol, kpa), klvntrm(rw, T) ); } @@ -186,7 +188,7 @@ namespace libcloudphxx BOOST_GPU_ENABLED real_t operator()( const real_t &rw2_old, - const thrust::tuple, real_t, real_t> &tpl + const thrust::tuple, real_t, real_t> &tpl ) const { #if !defined(__NVCC__) using std::min; @@ -221,7 +223,8 @@ namespace libcloudphxx "kpa: %g " "vt: %g " "lambda_D: %g " - "lambda_K: %g\n", + "lambda_K: %g\n" + "rd3_insol: %g ", drw2, rw2_old, dt, RH_max, thrust::get<0>(tpl_in), // rhod thrust::get<1>(tpl_in), // rv @@ -233,7 +236,8 @@ namespace libcloudphxx thrust::get<5>(tpl_in), // kpa thrust::get<6>(tpl_in), // vt thrust::get<7>(tpl_in), // lambda_D - thrust::get<8>(tpl_in) // lambda_K + thrust::get<8>(tpl_in), // lambda_K + thrust::get<9>(tpl_in) // rd3_insol ); assert(0); } @@ -282,10 +286,11 @@ namespace libcloudphxx "eta: %g " "rd3: %g " "kpa: %g " - "vt: %g\n", + "vt: %g\n" + "rd3_insol: %g ", a, b, drw2, rw2_old, rd2, dt, RH_max, thrust::get<0>(tpl_in),thrust::get<1>(tpl_in), thrust::get<2>(tpl_in),thrust::get<1>(tpl),thrust::get<2>(tpl),thrust::get<3>(tpl_in), - thrust::get<4>(tpl_in),thrust::get<5>(tpl_in),thrust::get<6>(tpl_in) + thrust::get<4>(tpl_in),thrust::get<5>(tpl_in),thrust::get<6>(tpl_in), thrust::get<9>(tpl_in) ); assert(0); } diff --git a/src/impl/condensation/percell/particles_impl_cond.ipp b/src/impl/condensation/percell/particles_impl_cond.ipp index 19d82878..fb29de15 100644 --- a/src/impl/condensation/percell/particles_impl_cond.ipp +++ b/src/impl/condensation/percell/particles_impl_cond.ipp @@ -54,7 +54,8 @@ namespace libcloudphxx 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_K.begin(), ijk.begin()), + rd3_insol.begin() )); // calculating drop growth in a timestep using backward Euler diff --git a/src/impl/condensation/perparticle/perparticle_advance_rw2.ipp b/src/impl/condensation/perparticle/perparticle_advance_rw2.ipp index 8386c1c5..be6d1a37 100644 --- a/src/impl/condensation/perparticle/perparticle_advance_rw2.ipp +++ b/src/impl/condensation/perparticle/perparticle_advance_rw2.ipp @@ -27,7 +27,8 @@ namespace libcloudphxx 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_K.begin(), ijk.begin()), + rd3_insol.begin() )); thrust::transform( diff --git a/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp b/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp index 6d0c2187..bc9f7a75 100644 --- a/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp +++ b/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp @@ -61,6 +61,7 @@ namespace libcloudphxx const real_t lambda_D = thrust::get<4>(thrust::get<1>(tpl)); const real_t lambda_K = thrust::get<5>(thrust::get<1>(tpl)); const real_t rd3 = thrust::get<6>(thrust::get<1>(tpl)); + const real_t rd3_insol = thrust::get<7>(thrust::get<1>(tpl)); const real_t kpa = thrust::get<0>(thrust::get<2>(tpl)); const real_t vt = thrust::get<1>(thrust::get<2>(tpl)); const real_t dot_ssp = turb_cond ? thrust::get<0>(thrust::get<1>(tpl)) : 0; @@ -155,7 +156,8 @@ namespace libcloudphxx kpa, vt, lambda_D, - lambda_K + lambda_K, + rd3_insol ), sstp_tmp_p, RH @@ -230,7 +232,8 @@ namespace libcloudphxx kpa, vt, lambda_D, - lambda_K + lambda_K, + rd3_insol ), sstp_tmp_p, RH @@ -314,7 +317,8 @@ namespace libcloudphxx thrust::make_permutation_iterator(dv.begin(), ijk.begin()), thrust::make_permutation_iterator(lambda_D.begin(), ijk.begin()), thrust::make_permutation_iterator(lambda_K.begin(), ijk.begin()), - rd3.begin() + rd3.begin(), + rd3_insol.begin() )), thrust::make_zip_iterator(thrust::make_tuple( kpa.begin(), diff --git a/src/impl/diagnose_SD_attributes/particles_impl_fill_outbuf.ipp b/src/impl/diagnose_SD_attributes/particles_impl_fill_outbuf.ipp index 82f8e3db..ff86b1cd 100644 --- a/src/impl/diagnose_SD_attributes/particles_impl_fill_outbuf.ipp +++ b/src/impl/diagnose_SD_attributes/particles_impl_fill_outbuf.ipp @@ -39,12 +39,12 @@ namespace libcloudphxx template std::vector particles_t::impl::fill_attr_outbuf(const std::string &name) { - const std::set attr_names = {"rw2", "rd3", "kappa", "rd2_insol", "T_freeze", "ice_a", "ice_c", "ice_rho", "x", "y", "z"}; // TODO implement "n" - it is n_t type and others are real_t + const std::set attr_names = {"rw2", "rd3", "kappa", "rd3_insol", "T_freeze", "ice_a", "ice_c", "ice_rho", "x", "y", "z"}; // TODO implement "n" - it is n_t type and others are real_t if (std::find(std::begin(attr_names), std::end(attr_names), name) == std::end(attr_names)) throw std::runtime_error("Unknown attribute name passed to get_attr."); if (!opts_init.ice_switch && - (name == "ice_a" || name == "ice_c" || name == "ice_rho" || name == "rd2_insol" || name == "T_freeze")) + (name == "ice_a" || name == "ice_c" || name == "ice_rho" || name == "T_freeze")) { throw std::runtime_error("Requested ice attribute '" + name + "' but ice_switch is off."); } @@ -57,7 +57,7 @@ namespace libcloudphxx name == "rw2" ? rw2 : name == "rd3" ? rd3 : name == "kappa" ? kpa : - name == "rd2_insol" ? rd2_insol : + name == "rd3_insol" ? rd3_insol : name == "T_freeze" ? T_freeze : name == "ice_a" ? ice_a : name == "ice_c" ? ice_c : diff --git a/src/impl/diagnose_SD_attributes/particles_impl_update_incloud_time.ipp b/src/impl/diagnose_SD_attributes/particles_impl_update_incloud_time.ipp index 66a6a06a..7d8922c1 100644 --- a/src/impl/diagnose_SD_attributes/particles_impl_update_incloud_time.ipp +++ b/src/impl/diagnose_SD_attributes/particles_impl_update_incloud_time.ipp @@ -43,7 +43,8 @@ namespace libcloudphxx thrust::transform( rd3.begin(), rd3.end(), // input - 1st arg thrust::make_zip_iterator(make_tuple( - kpa.begin(), + rd3_insol.begin(), + kpa.begin(), thrust::make_permutation_iterator( T.begin(), ijk.begin() diff --git a/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp b/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp index c28b7706..6ad3b1a8 100644 --- a/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp +++ b/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp @@ -21,6 +21,7 @@ namespace libcloudphxx thrust::transform_if( rd3.begin(), rd3.end(), // input - 1st arg thrust::make_zip_iterator(make_tuple( + rd3_insol.begin(), kpa.begin(), thrust::make_constant_iterator(opts_init.rc2_T + 273.15) )), // input - 2nd arg diff --git a/src/impl/ice/particles_impl_ice_nucl_melt.ipp b/src/impl/ice/particles_impl_ice_nucl_melt.ipp index 5f7f32cd..00661320 100644 --- a/src/impl/ice/particles_impl_ice_nucl_melt.ipp +++ b/src/impl/ice/particles_impl_ice_nucl_melt.ipp @@ -62,7 +62,7 @@ namespace libcloudphxx BOOST_GPU_ENABLED void operator()(thrust::tuple< real_t&, real_t&, real_t&, real_t&, // to be updated (rw2, a, c, rho_i) - const real_t&, const real_t&, const real_t& // rd2_insol, u01, T + const real_t&, const real_t&, const real_t& // rd3_insol, u01, T > tpl) const { auto& rw2 = thrust::get<0>(tpl); @@ -70,11 +70,11 @@ namespace libcloudphxx auto& c = thrust::get<2>(tpl); auto& rho_i = thrust::get<3>(tpl); - const real_t rd2_insol = thrust::get<4>(tpl); + const real_t rd3_insol = thrust::get<4>(tpl); const real_t u01 = thrust::get<5>(tpl); const real_t T = thrust::get<6>(tpl); - if (rw2 > real_t(0) && u01 < common::ice_nucleation::p_freeze(INP_type, rd2_insol, rw2, T, dt)) + if (rw2 > real_t(0) && u01 < common::ice_nucleation::p_freeze(INP_type, rd3_insol, rw2, T, dt)) { rho_i = common::moist_air::rho_i() * si::cubic_metres / si::kilograms; a = pow(rw2, real_t(0.5)) * pow(common::moist_air::rho_w() / common::moist_air::rho_i(), real_t(1./3.)); @@ -148,7 +148,7 @@ namespace libcloudphxx ice_a.begin(), ice_c.begin(), ice_rho.begin(), - rd2_insol.begin(), + rd3_insol.begin(), u01.begin(), thrust::make_permutation_iterator(T.begin(), ijk.begin()) )), @@ -157,7 +157,7 @@ namespace libcloudphxx ice_a.begin(), ice_c.begin(), ice_rho.begin(), - rd2_insol.begin(), + rd3_insol.begin(), u01.begin(), thrust::make_permutation_iterator(T.begin(), ijk.begin()) )) + n_part, 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 0585bcb9..64345155 100644 --- a/src/impl/initialization/particles_impl_init_SD_with_distros.ipp +++ b/src/impl/initialization/particles_impl_init_SD_with_distros.ipp @@ -54,10 +54,10 @@ namespace libcloudphxx { // init kappa init_kappa(kpa_rd_insol.kappa); + init_insol_dry_sizes(kpa_rd_insol.rd_insol); if (opts_init.ice_switch) { - init_insol_dry_sizes(kpa_rd_insol.rd_insol); init_a_c_rho_ice(); if (! opts_init.time_dep_ice_nucl) { 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 93c45049..19f3726f 100644 --- a/src/impl/initialization/particles_impl_init_SD_with_sizes.ipp +++ b/src/impl/initialization/particles_impl_init_SD_with_sizes.ipp @@ -47,10 +47,10 @@ namespace libcloudphxx // init kappa and rd_insol init_kappa(kappa); + init_insol_dry_sizes(rd_insol); if (opts_init.ice_switch) { - init_insol_dry_sizes(rd_insol); init_a_c_rho_ice(); if (! opts_init.time_dep_ice_nucl) { diff --git a/src/impl/initialization/particles_impl_init_T_freeze.ipp b/src/impl/initialization/particles_impl_init_T_freeze.ipp index ee0f5da5..9958342d 100644 --- a/src/impl/initialization/particles_impl_init_T_freeze.ipp +++ b/src/impl/initialization/particles_impl_init_T_freeze.ipp @@ -23,8 +23,8 @@ namespace libcloudphxx rand_u01(u01, n_part_to_init); thrust::transform( - thrust::make_zip_iterator(thrust::make_tuple(rd2_insol.begin() + n_part_old, u01.begin() + n_part_old)), - thrust::make_zip_iterator(thrust::make_tuple(rd2_insol.end(), u01.end())), + thrust::make_zip_iterator(thrust::make_tuple(rd3_insol.begin() + n_part_old, u01.begin())), + thrust::make_zip_iterator(thrust::make_tuple(rd3_insol.begin() + n_part_old, u01.begin())) + n_part_to_init, T_freeze.begin() + n_part_old, T_freeze_CDF_inv_functor(opts_init.inp_type) ); diff --git a/src/impl/initialization/particles_impl_init_insol_dry_sizes.ipp b/src/impl/initialization/particles_impl_init_insol_dry_sizes.ipp index 3facc1a9..074391b9 100644 --- a/src/impl/initialization/particles_impl_init_insol_dry_sizes.ipp +++ b/src/impl/initialization/particles_impl_init_insol_dry_sizes.ipp @@ -16,8 +16,8 @@ namespace libcloudphxx real_t radius ) { - real_t rad2 = radius * radius; - thrust::fill(rd2_insol.begin() + n_part_old, rd2_insol.end(), rad2); + real_t rad3 = radius * radius * radius; + thrust::fill(rd3_insol.begin() + n_part_old, rd3_insol.end(), rad3); } }; }; diff --git a/src/impl/initialization/particles_impl_init_wet.ipp b/src/impl/initialization/particles_impl_init_wet.ipp index 85dd9607..ed40d8eb 100644 --- a/src/impl/initialization/particles_impl_init_wet.ipp +++ b/src/impl/initialization/particles_impl_init_wet.ipp @@ -22,19 +22,20 @@ namespace libcloudphxx rw2_eq(const real_t &RH_max) : RH_max(RH_max) {} BOOST_GPU_ENABLED - real_t operator()(const thrust::tuple &tpl) + real_t operator()(const thrust::tuple &tpl) { #if !defined(__NVCC__) using std::min; using std::pow; #endif const quantity rd3 = thrust::get<0>(tpl) * si::cubic_metres; - const quantity kpa = thrust::get<1>(tpl); - const quantity RH = min(thrust::get<2>(tpl), RH_max); - const quantity T = thrust::get<3>(tpl) * si::kelvins; + const quantity rd3_insol = thrust::get<1>(tpl) * si::cubic_metres; + const quantity kpa = thrust::get<2>(tpl); + const quantity RH = min(thrust::get<3>(tpl), RH_max); + const quantity T = thrust::get<4>(tpl) * si::kelvins; return pow(common::kappa_koehler::rw3_eq( - rd3, kpa, RH, T + rd3, rd3_insol, kpa, RH, T ) / si::cubic_metres, real_t(2./3)); } }; @@ -54,6 +55,7 @@ namespace libcloudphxx typedef thrust::zip_iterator< thrust::tuple< + typename thrust_device::vector::iterator, typename thrust_device::vector::iterator, typename thrust_device::vector::iterator, pi_t, @@ -62,7 +64,8 @@ namespace libcloudphxx > zip_it_t; zip_it_t zip_it(thrust::make_tuple( - rd3.begin() + n_part_old, + rd3.begin() + n_part_old, + rd3_insol.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) diff --git a/src/impl/initialization/particles_impl_reserve_hskpng_npart.ipp b/src/impl/initialization/particles_impl_reserve_hskpng_npart.ipp index aa2a9d98..e43bb2b2 100644 --- a/src/impl/initialization/particles_impl_reserve_hskpng_npart.ipp +++ b/src/impl/initialization/particles_impl_reserve_hskpng_npart.ipp @@ -39,7 +39,6 @@ namespace libcloudphxx } if(opts_init.ice_switch) { - rd2_insol.reserve(opts_init.n_sd_max); ice_a.reserve(opts_init.n_sd_max); ice_c.reserve(opts_init.n_sd_max); ice_rho.reserve(opts_init.n_sd_max); @@ -64,6 +63,7 @@ namespace libcloudphxx rw2.reserve(opts_init.n_sd_max); n.reserve(opts_init.n_sd_max); kpa.reserve(opts_init.n_sd_max); + rd3_insol.reserve(opts_init.n_sd_max); if(opts_init.diag_incloud_time) incloud_time.reserve(opts_init.n_sd_max); diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 681bf689..12b81113 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -92,7 +92,7 @@ namespace libcloudphxx dv, // grid-cell volumes (per grid cell) incloud_time, // time this SD has been within a cloud rc2, // critical radius squared (estimated for temperature from opts_init.rc2_T) - rd2_insol, // dry radii squared of insoluble aerosol + rd3_insol, // dry radii cubed of insoluble aerosol T_freeze, // freezing temperature ice_a, // equatorial radius of ice ice_c, // polar radius of ice @@ -442,6 +442,7 @@ namespace libcloudphxx distmem_real_vctrs.insert({&rd3, detail::no_initial_value}); distmem_real_vctrs.insert({&rw2, detail::no_initial_value}); distmem_real_vctrs.insert({&kpa, detail::no_initial_value}); + distmem_real_vctrs.insert({&rd3_insol, detail::no_initial_value}); distmem_real_vctrs.insert({&vt, detail::invalid}); @@ -477,7 +478,6 @@ namespace libcloudphxx if(opts_init.ice_switch) { - distmem_real_vctrs.insert({&rd2_insol, detail::no_initial_value}); distmem_real_vctrs.insert({&ice_a, detail::no_initial_value}); distmem_real_vctrs.insert({&ice_c, detail::no_initial_value}); distmem_real_vctrs.insert({&ice_rho, detail::no_initial_value}); diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index 720a3935..5c95c676 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -41,17 +41,19 @@ namespace libcloudphxx struct rw3_cr { BOOST_GPU_ENABLED - real_t operator()(const real_t &rd3, const thrust::tuple &tpl) + real_t operator()(const real_t &rd3, const thrust::tuple &tpl) { - const quantity kpa = thrust::get<0>(tpl); - const quantity T = thrust::get<1>(tpl) * si::kelvins; + const quantity rd3_insol = thrust::get<0>(tpl) * si::cubic_meters; + const quantity kpa = thrust::get<1>(tpl); + const quantity T = thrust::get<2>(tpl) * si::kelvins; #if !defined(__NVCC__) using std::pow; #endif return pow( common::kappa_koehler::rw3_cr( - rd3 * si::cubic_metres, + rd3 * si::cubic_metres, + rd3_insol, kpa, T ) / si::cubic_metres @@ -96,14 +98,16 @@ namespace libcloudphxx struct RH_minus_Sc { BOOST_GPU_ENABLED - real_t operator()(const real_t &rd3, const thrust::tuple &tpl) + real_t operator()(const real_t &rd3, const thrust::tuple &tpl) { - const quantity kpa = thrust::get<0>(tpl); - const quantity T = thrust::get<1>(tpl) * si::kelvins; - const quantity RH = thrust::get<2>(tpl); + const quantity rd3_insol = thrust::get<0>(tpl) * si::cubic_meters; + const quantity kpa = thrust::get<1>(tpl); + const quantity T = thrust::get<2>(tpl) * si::kelvins; + const quantity RH = thrust::get<3>(tpl); return RH - common::kappa_koehler::S_cr( - rd3 * si::cubic_metres, + rd3 * si::cubic_metres, + rd3_insol, kpa, T ); @@ -361,7 +365,8 @@ namespace libcloudphxx thrust::transform( pimpl->rd3.begin(), pimpl->rd3.end(), // input - 1st arg thrust::make_zip_iterator(make_tuple( - pimpl->kpa.begin(), + pimpl->rd3_insol.begin(), + pimpl->kpa.begin(), thrust::make_permutation_iterator( pimpl->T.begin(), pimpl->ijk.begin() @@ -392,7 +397,8 @@ namespace libcloudphxx thrust::transform( pimpl->rd3.begin(), pimpl->rd3.end(), // input - 1st arg thrust::make_zip_iterator(make_tuple( - pimpl->kpa.begin(), + pimpl->rd3_insol.begin(), + pimpl->kpa.begin(), thrust::make_permutation_iterator( pimpl->T.begin(), pimpl->ijk.begin() diff --git a/tests/python/unit/api_common.py b/tests/python/unit/api_common.py index 360eed31..68756292 100644 --- a/tests/python/unit/api_common.py +++ b/tests/python/unit/api_common.py @@ -24,8 +24,8 @@ assert common.th_std2dry(common.th_dry2std(th, rv), rv) == th rd3 = (.2e-6)**3 -assert common.rw3_cr(rd3, .5, 300) > rd3 -assert common.S_cr(rd3, .5, 300) > 1 +assert common.rw3_cr(rd3, 0, .5, 300) > rd3 +assert common.S_cr(rd3, 0, .5, 300) > 1 # just testing if pressure at 200m is lower than at 100m assert common.p_hydro(100, 300, .01, 0, 100000) > common.p_hydro(200, 300, .01, 0, 100000) diff --git a/tests/python/unit/diag_incloud_time.py b/tests/python/unit/diag_incloud_time.py index a848f38d..de78274b 100644 --- a/tests/python/unit/diag_incloud_time.py +++ b/tests/python/unit/diag_incloud_time.py @@ -19,10 +19,15 @@ def lognormal(lnr): opts = lgrngn.opts_t() +kappa1 = 0.61 +kappa2 = 1.28 +rd_insol = 0. + opts_init = lgrngn.opts_init_t() -opts_init.dry_distros = {(.61, 0.):lognormal, (1.28, 0.):lognormal} +opts_init.dry_distros = {(kappa1, rd_insol):lognormal, (kappa2, rd_insol):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 = 0.1 opts_init.sd_conc = int(1e2) diff --git a/tests/python/unit/source.py b/tests/python/unit/source.py index 04f30660..26629246 100644 --- a/tests/python/unit/source.py +++ b/tests/python/unit/source.py @@ -30,7 +30,7 @@ def lognormal_src(lnr): ) / log(stdev) / sqrt(2*pi); kappa = .61 -rd_insol = 0.5e-6 +rd_insol = 0. def test(opts_init, opts): opts_init.rng_seed = int(time())