From 213f53e4be419d96a4f1070a696ff55fb5b6a286 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Tue, 3 Feb 2026 16:45:50 +0100 Subject: [PATCH 01/21] ice growth formulas --- .../libcloudph++/common/ice_deposition.hpp | 115 ++++++++++ .../ice/particles_impl_ice_dep_common.ipp | 200 ++++++++++++++++++ 2 files changed, 315 insertions(+) create mode 100644 include/libcloudph++/common/ice_deposition.hpp create mode 100644 src/impl/ice/particles_impl_ice_dep_common.ipp diff --git a/include/libcloudph++/common/ice_deposition.hpp b/include/libcloudph++/common/ice_deposition.hpp new file mode 100644 index 00000000..afba553c --- /dev/null +++ b/include/libcloudph++/common/ice_deposition.hpp @@ -0,0 +1,115 @@ +#pragma once + +#include "units.hpp" +#include "const_cp.hpp" +#include + +#if defined(__NVCC__) +# include +#endif + +namespace libcloudphxx +{ + namespace common + { + namespace ice_deposition + { + + // capacitance + template + BOOST_GPU_ENABLED + quantity capacitance( + const quantity ice_a, // ice equatorial radius + const quantity ice_c // ice polar radius + ) + { + if (ice_a < std::numeric_limits::epsilon) throw std::runtime_error("ice_a too small for calculation of capacitance"); + real_t phi = ice_c / ice_a; + real_t e = std::sqrt(real_t(1) - phi * phi); // eccentricity + return phi < real_t(1) ? + ice_a * e / std::asin(e) : ice_a * e / std::log(1 + e) / phi; + } + + + // rate of change of ice mass + template + BOOST_GPU_ENABLED + quantity::type, real_t> dmi_dt( + const quantity D, // D + const quantity K, // K + const quantity rho_v, // ambient water vapour density + const quantity T, // ambient temperature + const quantity p, // ambient pressure + const quantity RH_i, // p_v/p_vsi = relative humidity w.r.t. ice + const quantity ice_a, // ice equatorial radius + const quantity ice_c // ice polar radius + ) + { + using moist_air::R_v; + + quantity::type, real_t> + l_s = const_cp::l_s(T); + + real_t C = capacitance(ice_a, ice_c); + + return real_t(4) * pi() * C * (real_t(1) - real_t(1) / RH_i) + / ( + real_t(1) + / D + / rho_v + + + l_s + / K + / RH_i + / T + * (l_s / R_v() / T - real_t(1)) + ) + ; + } + + + // deposition density Chen and Lamb + template + BOOST_GPU_ENABLED + quantity rho_dep( + const quantity D, // D + const quantity K, // K + const quantity rho_v, // ambient water vapour density + const quantity T, // ambient temperature + const quantity p, // ambient pressure + const quantity RH_i, // p_v/p_vsi = relative humidity w.r.t. ice + const quantity ice_a, // ice equatorial radius + const quantity ice_c // ice polar radius + ) + { + + return (lambda(T)< real_t(1) && ice_a < real_t(1e-4)*si::meters) ? + moist_air::rho_i<>() : rho_dep_CL94(S, T, p, RH_i); + } + + + // deposition density + template + BOOST_GPU_ENABLED + quantity rho_dep( + const quantity D, // D + const quantity K, // K + const quantity rho_v, // ambient water vapour density + const quantity T, // ambient temperature + const quantity p, // ambient pressure + const quantity RH_i, // p_v/p_vsi = relative humidity w.r.t. ice + const quantity ice_a, // ice equatorial radius + const quantity ice_c // ice polar radius + ) + { + + return moist_air::rho_i<>() * std::exp(- real_t(3) * std::max(d_rhoi/si::kilograms*si::cubic_meters - real_t(5e-5), real_t(0))/gamma(T)); + } + + + // gamma(T) + + + }; + }; +}; diff --git a/src/impl/ice/particles_impl_ice_dep_common.ipp b/src/impl/ice/particles_impl_ice_dep_common.ipp new file mode 100644 index 00000000..d51bb322 --- /dev/null +++ b/src/impl/ice/particles_impl_ice_dep_common.ipp @@ -0,0 +1,200 @@ +// // vim:filetype=cpp +// /** @file +// * @copyright University of Warsaw +// * @section LICENSE +// * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) +// */ +// +// // #include +// #include +// #include +// #include +// #include +// #include +// #include +// #include +// +// namespace libcloudphxx +// { +// namespace lgrngn +// { +// namespace detail +// { +// +// template +// struct advance_ice_a_minfun +// { +// const quantity ice_a_old; +// const quantity dt; +// const quantity rhod; +// const quantity rv; +// const quantity T; +// const quantity p; +// const quantity RH_i; +// const quantity eta; +// const quantity rd3; +// const quantity kpa; +// const quantity vt; +// const quantity RH_max; +// const quantity lambda_D; +// const quantity lambda_K; +// +// // ctor +// BOOST_GPU_ENABLED +// advance_ice_a_minfun( +// const real_t &dt, +// const real_t &ice_a, +// const thrust::tuple, real_t, real_t> &tpl, +// const real_t &RH_max +// ) : +// dt(dt * si::seconds), +// ice_a_old(ice_a * si::square_metres), +// rhod( thrust::get<0>(thrust::get<0>(tpl)) * si::kilograms / si::cubic_metres), +// rv( thrust::get<1>(thrust::get<0>(tpl))), +// T( thrust::get<2>(thrust::get<0>(tpl)) * si::kelvins), +// eta( thrust::get<3>(thrust::get<0>(tpl)) * si::pascals * si::seconds), +// rd3( thrust::get<4>(thrust::get<0>(tpl)) * si::cubic_metres), +// kpa( thrust::get<5>(thrust::get<0>(tpl))), +// vt( thrust::get<6>(thrust::get<0>(tpl)) * si::metres_per_second), +// p( thrust::get<1>(tpl) * si::pascals), +// RH_i( 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), +// RH_max(RH_max) +// {} +// +// BOOST_GPU_ENABLED +// quantity::type, real_t> d_ice_a_dt(const quantity &ice_a) const +// { +// using namespace common::maxwell_mason; +// using namespace common::kappa_koehler; +// using namespace common::kelvin; +// using common::moist_air::D_0; +// using common::moist_air::K_0; +// using common::moist_air::c_pd; +// using common::transition_regime::beta; +// using common::ventil::Sh; +// using common::ventil::Nu; +// #if !defined(__NVCC__) +// using std::sqrt; +// #endif +// +// const quantity +// Re = common::ventil::Re(vt, ice_a, rhod, eta), +// Sc = common::ventil::Sc(eta, rhod, D_0()), // TODO? cache +// Pr = common::ventil::Pr(eta, c_pd(), K_0()); // TODO? cache +// +// const quantity +// D = D_0() * beta(lambda_D / ice_a) * (Sh(Sc, Re) / 2); +// +// const quantity +// K = K_0() * beta(lambda_K / ice_a) * (Nu(Pr, Re) / 2); +// +// return da_dt( +// D, +// K, +// rhod * rv, +// T, +// p, +// RH_i > RH_max ? RH_max : RH_i +// ); +// } +// +// BOOST_GPU_ENABLED +// real_t operator()(const real_t &ice_a_unitless) const +// { +// const quantity ice_a = ice_a_unitless * si::metres; +// return (ice_a_old + dt * d_ice_a_dt(ice_a) - ice_a) / si::metres; +// } +// }; +// +// +// +// template +// struct advance_ice_axis +// { +// const real_t dt, RH_max; +// +// advance_ice_axis(const real_t &dt, const real_t &RH_max) +// : dt(dt), RH_max(RH_max) {} +// +// BOOST_GPU_ENABLED +// thrust::tuple operator()( +// const thrust::tuple &ac_old, +// const thrust::tuple< +// thrust::tuple, +// real_t, real_t> &tpl +// ) const +// { +// #if !defined(__NVCC__) +// using std::max; +// using std::isnan; +// using std::isinf; +// #endif +// const real_t a_old = thrust::get<0>(ac_old); +// const real_t c_old = thrust::get<1>(ac_old); +// +// // Skip liquid droplets +// if (a_old <= real_t(0) || c_old <= real_t(0)) +// return ac_old; +// +// advance_ice_a_minfun f_a(dt, a_old * a_old, tpl, RH_max); +// advance_ice_c_minfun f_c(dt, c_old * c_old, tpl, RH_max); +// +// const real_t da_dt = (f_a.drw2_dt(a_old * a_old * si::square_metres) / (2 * a_old * si::metres)) +// * si::seconds / si::metres; +// const real_t dc_dt = (f_c.drw2_dt(c_old * c_old * si::square_metres) / (2 * c_old * si::metres)) +// * si::seconds / si::metres; +// +// // to store the result +// real_t a_new, c_new; +// +// if (da_dt == 0 && dc_dt == 0) return ac_old; +// +// const real_t rd = cbrt(thrust::get<4>(tpl_in)); +// +// const real_t +// a = max(rd2, rw2_old + min(real_t(0), cond_mlt * drw2)), +// b = rw2_old + max(real_t(0), cond_mlt * drw2); +// +// // numerics (drw2 != 0 but a==b) +// if (a == b) +// { +// if constexpr (apply) +// return rw2_old; +// else +// return real_t(0); +// } +// +// real_t fa, fb; +// +// if (drw2 > 0) +// { +// fa = drw2; // for implicit Euler its equal to min_fun(x_old) +// fb = f(b); +// } +// else +// { +// fa = f(a); +// fb = drw2; // for implicit Euler its equal to min_fun(x_old) +// } +// +// // root-finding ill posed => explicit Euler +// if (fa * fb > 0) rw2_new = rw2_old + drw2; +// // otherwise implicit Euler +// else +// { +// auto _n_iter = n_iter; // we need a copy because toms748_solve expects non-const ref (n_iter is modified by it) +// rw2_new = common::detail::toms748_solve(f, a, b, fa, fb, eps_tolerance, _n_iter); +// } +// +// // check if it doesn't evaporate too much +// if(rw2_new < rd2) rw2_new = rd2; +// +// return thrust::make_tuple(a_new, c_new, rho_new); +// } +// }; +// +// }; +// }; +// }; From ae1d70d5767915d149903ccf913b19049dec1376 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Wed, 4 Feb 2026 18:19:37 +0100 Subject: [PATCH 02/21] dep density formulas --- .../libcloudph++/common/ice_deposition.hpp | 63 ++++++++++++------- 1 file changed, 42 insertions(+), 21 deletions(-) diff --git a/include/libcloudph++/common/ice_deposition.hpp b/include/libcloudph++/common/ice_deposition.hpp index afba553c..549d681d 100644 --- a/include/libcloudph++/common/ice_deposition.hpp +++ b/include/libcloudph++/common/ice_deposition.hpp @@ -15,12 +15,12 @@ namespace libcloudphxx namespace ice_deposition { - // capacitance + // capacitance of ice crystals template BOOST_GPU_ENABLED - quantity capacitance( - const quantity ice_a, // ice equatorial radius - const quantity ice_c // ice polar radius + quantity ice_capacitance( + const quantity ice_a, // equatorial radius + const quantity ice_c // polar radius ) { if (ice_a < std::numeric_limits::epsilon) throw std::runtime_error("ice_a too small for calculation of capacitance"); @@ -50,9 +50,8 @@ namespace libcloudphxx quantity::type, real_t> l_s = const_cp::l_s(T); - real_t C = capacitance(ice_a, ice_c); - - return real_t(4) * pi() * C * (real_t(1) - real_t(1) / RH_i) + return real_t(4) * pi() * ice_capacitance(ice_a, ice_c) + * (real_t(1) - real_t(1) / RH_i) / ( real_t(1) / D @@ -68,7 +67,7 @@ namespace libcloudphxx } - // deposition density Chen and Lamb + // deposition density from Shima et al. (2020) template BOOST_GPU_ENABLED quantity rho_dep( @@ -76,38 +75,60 @@ namespace libcloudphxx const quantity K, // K const quantity rho_v, // ambient water vapour density const quantity T, // ambient temperature - const quantity p, // ambient pressure const quantity RH_i, // p_v/p_vsi = relative humidity w.r.t. ice - const quantity ice_a, // ice equatorial radius - const quantity ice_c // ice polar radius + const quantity ice_a // ice equatorial radius ) { - return (lambda(T)< real_t(1) && ice_a < real_t(1e-4)*si::meters) ? - moist_air::rho_i<>() : rho_dep_CL94(S, T, p, RH_i); + return (ice_growth_ratio(T) < real_t(1) && ice_a < real_t(1e-4)*si::meters) ? + moist_air::rho_i<>() : rho_dep_CL94(D, K, rho_v, T, RH_i); } - // deposition density + // deposition density from Chen and Lamb (1994) template BOOST_GPU_ENABLED - quantity rho_dep( + quantity rho_dep_CL94( const quantity D, // D const quantity K, // K const quantity rho_v, // ambient water vapour density const quantity T, // ambient temperature - const quantity p, // ambient pressure - const quantity RH_i, // p_v/p_vsi = relative humidity w.r.t. ice - const quantity ice_a, // ice equatorial radius - const quantity ice_c // ice polar radius + const quantity RH_i // p_v/p_vsi = relative humidity w.r.t. ice ) { + using moist_air::R_v; + quantity::type, real_t> + l_s = const_cp::l_s(T); + + using const_cp::p_vs, const_cp::p_vsi; + + quantity d_rhoi = (std::min(RH_i, p_vs/p_vsi) - real_t(1)) + / RH_i + / D + / ( + real_t(1) + / D + / rho_v + + + l_s + / K + / RH_i + / T + * (l_s / R_v() / T - real_t(1)) + ); - return moist_air::rho_i<>() * std::exp(- real_t(3) * std::max(d_rhoi/si::kilograms*si::cubic_meters - real_t(5e-5), real_t(0))/gamma(T)); + return moist_air::rho_i<>() + * std::exp(- real_t(3) * std::max(d_rhoi/si::kilograms*si::cubic_meters - real_t(5e-5), real_t(0)) + / ice_growth_ratio(T)); } - // gamma(T) + // inherent ice growth ratio + template + BOOST_GPU_ENABLED + quantity ice_growth_ratio(const quantity T) + { + } }; From fa2aff0114ff0e4c77c9ff7219900474a946e9c8 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Fri, 17 Apr 2026 11:22:07 +0200 Subject: [PATCH 03/21] kinematic_2D ice options --- models/kinematic_2D/cases/icmw8_case1.hpp | 1 + .../kinematic_2D/src/kin_cloud_2d_lgrngn.hpp | 35 +++++++++++++++++-- models/kinematic_2D/src/opts_lgrngn.hpp | 25 +++++++++---- 3 files changed, 52 insertions(+), 9 deletions(-) diff --git a/models/kinematic_2D/cases/icmw8_case1.hpp b/models/kinematic_2D/cases/icmw8_case1.hpp index 3f44a4f7..76133a3f 100644 --- a/models/kinematic_2D/cases/icmw8_case1.hpp +++ b/models/kinematic_2D/cases/icmw8_case1.hpp @@ -46,6 +46,7 @@ namespace config //aerosol chemical composition parameters (needed for activation) // for lgrngn: quantity kappa; // CCN-derived value from Table 1 in Petters and Kreidenweis 2007 + quantity rd_insol; // insoluble dry radius // for blk_2m: quantity chem_b; //ammonium sulphate //chem_b = 1.33; // sodium chloride // for lagrangian simulations with aq. chemistry diff --git a/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp b/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp index 29c96498..b228fd48 100644 --- a/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp +++ b/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp @@ -27,7 +27,7 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common // member fields std::unique_ptr> prtcls; - bool coal, sedi; + bool coal, sedi, ice_switch, ice_nucl, time_dep_ice_nucl, depo; // helper methods void diag() @@ -79,6 +79,32 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common } rng_num++; } + // ice_a + int rng_num = 0; + for (auto &rng_moms : params.out_ice) + { + auto &rng(rng_moms.first); + prtcls->diag_ice_a_rng(rng.first / si::metres, rng.second / si::metres); + for (auto &mom : rng_moms.second) + { + prtcls->diag_ice_a_mom(mom); + this->record_aux(aux_name("ice_a", rng_num, mom), prtcls->outbuf()); + } + rng_num++; + } + // ice_c + int rng_num = 0; + for (auto &rng_moms : params.out_ice) + { + auto &rng(rng_moms.first); + prtcls->diag_ice_c_rng(rng.first / si::metres, rng.second / si::metres); + for (auto &mom : rng_moms.second) + { + prtcls->diag_ice_c_mom(mom); + this->record_aux(aux_name("ice_c", rng_num, mom), prtcls->outbuf()); + } + rng_num++; + } } { // rw3(rd) @@ -129,6 +155,10 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common { coal = params.cloudph_opts.coal; sedi = params.cloudph_opts.sedi; + ice_switch = params.cloudph_opts_init.ice_switch; + ice_nucl = params.cloudph_opts.ice_nucl; + time_dep_ice_nucl = params.cloudph_opts_init.time_dep_ice_nucl; + depo = params.cloudph_opts.depo; parent_t::hook_ante_loop(nt); @@ -151,6 +181,7 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common this->record_aux_const("n1_stp", this->setup.n1_stp * si::cubic_metres); this->record_aux_const("n2_stp", this->setup.n2_stp * si::cubic_metres); this->record_aux_const("kappa", this->setup.kappa); + this->record_aux_const("rd_insol", this->setup.rd_insol / si::metres); assert(params.backend != -1); assert(params.dt != 0); @@ -302,7 +333,7 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common bool async = true; libcloudphxx::lgrngn::opts_t cloudph_opts; libcloudphxx::lgrngn::opts_init_t cloudph_opts_init; - outmom_t out_dry, out_wet; + outmom_t out_dry, out_wet, out_ice; }; protected: diff --git a/models/kinematic_2D/src/opts_lgrngn.hpp b/models/kinematic_2D/src/opts_lgrngn.hpp index 924921d2..d6635168 100644 --- a/models/kinematic_2D/src/opts_lgrngn.hpp +++ b/models/kinematic_2D/src/opts_lgrngn.hpp @@ -76,7 +76,8 @@ void parse_moms( std::vector> min_maxnum; outmom_t &moms = opt == "out_dry" ? rt_params.out_dry : - rt_params.out_wet; + opt == "out_wet" ? rt_params.out_wet : + rt_params.out_ice; const bool result = qi::phrase_parse(first, last, *( @@ -140,9 +141,10 @@ void parse_moms( outmom_t &moms = opt == "out_dry" ? rt_params.out_dry : - opt == "out_wet" ? rt_params.out_wet : - opt == "out_chem" ? rt_params.out_chem : - rt_params.out_wet_pH; + opt == "out_wet" ? rt_params.out_wet : + opt == "out_ice" ? rt_params.out_ice : + opt == "out_chem" ? rt_params.out_chem : + rt_params.out_wet_pH; const bool result = qi::phrase_parse(first, last, *( @@ -198,7 +200,7 @@ void setopts_micro_chem( {solver_t::ix::th, {"th", "[K]"}}, {solver_t::ix::rv, {"rv", "[kg kg-1]"}} }; - out_set = {"out_dry", "out_wet"}; + out_set = {"out_dry", "out_wet", "out_ice"}; } template @@ -217,7 +219,7 @@ void setopts_micro_chem( {solver_t::ix::NH3g, {"NH3g", "[dimesnionless]"}}, {solver_t::ix::HNO3g, {"HNO3g","[dimensionless]"}} }; - out_set = {"out_dry", "out_wet", "out_chem", "out_wet_pH"}; + out_set = {"out_dry", "out_wet", "out_ice", "out_chem", "out_wet_pH"}; } // simulation and output parameters for micro=lgrngn @@ -246,6 +248,10 @@ void setopts_micro( ("chem_dsc", po::value()->default_value(rt_params.cloudph_opts.chem_dsc) , "dissociation (1=on, 0=off)") ("chem_rct", po::value()->default_value(rt_params.cloudph_opts.chem_rct) , "aqueous chemistry (1=on, 0=off)") ("chem_switch", po::value()->default_value(rt_params.cloudph_opts_init.chem_switch) , "aqueous chemistry (1=on, 0=off)") + ("ice_switch", po::value()->default_value(rt_params.cloudph_opts_init.ice_switch) , "enable ice (1=on, 0=off)") + ("ice_nucl", po::value()->default_value(rt_params.cloudph_opts.ice_nucl) , "ice nucleation (1=on, 0=off)") + ("time_dep_ice_nucl", po::value()->default_value(rt_params.cloudph_opts_init.time_dep_ice_nucl) , "time dependent ice nucleation (1=on, 0=off)") + ("depo", po::value()->default_value(rt_params.cloudph_opts.depo) , "ice depositional growth (1=on, 0=off)") // free parameters ("sstp_cond", po::value()->default_value(rt_params.cloudph_opts_init.sstp_cond), "no. of substeps for condensation") ("sstp_coal", po::value()->default_value(rt_params.cloudph_opts_init.sstp_coal), "no. of substeps for coalescence") @@ -255,6 +261,7 @@ void setopts_micro( // output ("out_dry", po::value()->default_value("0:1|0"), "dry radius ranges and moment numbers (r1:r2|n1,n2...;...)") ("out_wet", po::value()->default_value(".5e-6:25e-6|0,1,2,3;25e-6:1|0,3,6"), "wet radius ranges and moment numbers (r1:r2|n1,n2...;...)") + ("out_ice", po::value()->default_value(".1e-6:1|0,1,2,3"), "ice semi-axis ranges and moment numbers (r1:r2|n1,n2...;...)") ("out_wet_pH", po::value()->default_value("0:1|0"), "wet radius ranges for output of H+ and S_VI)") ("out_chem", po::value()->default_value("0:1|0"), "dry radius ranges for which chem mass is outputted") // collision and sedimentation @@ -297,7 +304,7 @@ void setopts_micro( */ rt_params.cloudph_opts_init.dry_distros.emplace( - libcloudphxx::lgrngn::kappa_rd_insol_t{setup.kappa, config::real_t(0.)}, // kappa, rd_insol + libcloudphxx::lgrngn::kappa_rd_insol_t{setup.kappa, setup.rd_insol / si::meters}, // kappa, rd_insol std::make_shared> (setup) ); @@ -306,6 +313,10 @@ void setopts_micro( rt_params.cloudph_opts.sedi = vm["sedi"].as(); rt_params.cloudph_opts.cond = vm["cond"].as(); rt_params.cloudph_opts.coal = vm["coal"].as(); + rt_params.cloudph_opts_init.ice_switch = vm["ice_switch"].as(); + rt_params.cloudph_opts.ice_nucl = vm["ice_nucl"].as(); + rt_params.cloudph_opts_init.time_dep_ice_nucl = vm["time_dep_ice_nucl"].as(); + rt_params.cloudph_opts.depo = vm["depo"].as(); rt_params.cloudph_opts.rcyc = vm["rcyc"].as(); rt_params.cloudph_opts.chem_dsl = vm["chem_dsl"].as(); From ed2331c1a859e302c0d4d631335fb3fc0d015646 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Fri, 17 Apr 2026 15:04:53 +0200 Subject: [PATCH 04/21] fix u01 for T_freeze --- src/impl/initialization/particles_impl_init_T_freeze.ipp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/impl/initialization/particles_impl_init_T_freeze.ipp b/src/impl/initialization/particles_impl_init_T_freeze.ipp index ee0f5da5..ce2097db 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(rd2_insol.begin() + n_part_old, u01.begin())), + thrust::make_zip_iterator(thrust::make_tuple(rd2_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) ); From 0e3a860bf8d02c3a1142443e2c9fcb2830194c40 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Wed, 22 Apr 2026 13:25:21 +0200 Subject: [PATCH 05/21] rd_insol in Kohler functions --- include/libcloudph++/common/kappa_koehler.hpp | 43 ++++++++++++------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/include/libcloudph++/common/kappa_koehler.hpp b/include/libcloudph++/common/kappa_koehler.hpp index 4292937c..76182287 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(rd3_insol < rw3); + 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 @@ -82,7 +87,7 @@ namespace libcloudphxx using std::cbrt; #endif return this->RH - - a_w(rw3 * si::cubic_metres, this->rd3, this->kappa) + - a_w(rw3 * si::cubic_metres, this->rd3, this->rd3_insol, this->kappa) * kelvin::klvntrm(real_t(cbrt(rw3)) * si::metres, this->T); } }; @@ -92,6 +97,7 @@ namespace libcloudphxx struct rw3_cr_minfun { const quantity rd3; + const quantity rd3_insol; const quantity kappa; const quantity T; @@ -99,9 +105,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 +118,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 +134,8 @@ namespace libcloudphxx template BOOST_GPU_ENABLED quantity rw3_eq( - quantity rd3, + quantity rd3, + quantity rd3_insol, quantity kappa, quantity RH, quantity T @@ -139,7 +147,7 @@ namespace libcloudphxx if(kappa == 0) return rd3; return common::detail::toms748_solve( - detail::rw3_eq_minfun(RH, rd3, kappa, T), // the above-defined functor + detail::rw3_eq_minfun(RH, rd3, rd3_insol, 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 ) * si::cubic_metres; @@ -151,7 +159,8 @@ namespace libcloudphxx template BOOST_GPU_ENABLED quantity rw3_cr( - quantity rd3, + quantity rd3, + quantity rd3_insol, quantity kappa, quantity T ) @@ -159,10 +168,11 @@ 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 + detail::rw3_cr_minfun(rd3_dbl, rd3_insol_dbl, double(kappa), T_dbl), // the above-defined functor double(1e0 * (rd3 / si::cubic_metres)), // min double(1e8 * (rd3 / si::cubic_metres)) // max )) * si::cubic_metres; @@ -173,6 +183,7 @@ namespace libcloudphxx BOOST_GPU_ENABLED quantity S_cr( quantity rd3, + quantity rd3_insol, quantity kappa, quantity T ) @@ -182,9 +193,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 ); From c1c421dc1f0a8b104b1e012b3cd6ab04690517b2 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Thu, 23 Apr 2026 12:34:01 +0200 Subject: [PATCH 06/21] init wet with rd_insol --- src/impl/initialization/particles_impl_init_wet.ipp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/impl/initialization/particles_impl_init_wet.ipp b/src/impl/initialization/particles_impl_init_wet.ipp index 85dd9607..6e266409 100644 --- a/src/impl/initialization/particles_impl_init_wet.ipp +++ b/src/impl/initialization/particles_impl_init_wet.ipp @@ -22,7 +22,7 @@ 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; @@ -32,9 +32,10 @@ namespace libcloudphxx 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<4>(tpl) * si::cubic_metres; return pow(common::kappa_koehler::rw3_eq( - rd3, kpa, RH, T + rd3, rd3_insol, kpa, RH, T ) / si::cubic_metres, real_t(2./3)); } }; From 6ef77f800899ceb5dfe462a818eb96ae774cae30 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Thu, 23 Apr 2026 18:32:56 +0200 Subject: [PATCH 07/21] update calls to functions with rd_insol --- bindings/python/common.hpp | 6 ++++-- .../common/particles_impl_cond_common.ipp | 6 ++++-- .../particles_impl_update_incloud_time.ipp | 3 ++- .../housekeeping/particles_impl_hskpng_rc2.ipp | 3 ++- src/particles_diag.ipp | 14 +++++++++----- 5 files changed, 21 insertions(+), 11 deletions(-) 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/src/impl/condensation/common/particles_impl_cond_common.ipp b/src/impl/condensation/common/particles_impl_cond_common.ipp index 49b61bc9..9dd3c734 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) ); } 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..5a3eec8d 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 @@ -47,7 +47,8 @@ namespace libcloudphxx thrust::make_permutation_iterator( T.begin(), ijk.begin() - ) + ), + rd3_insol.begin() )), // input - 2nd arg rc2.begin(), // output detail::rw3_cr() // op diff --git a/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp b/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp index c28b7706..02a36727 100644 --- a/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp +++ b/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp @@ -22,7 +22,8 @@ namespace libcloudphxx rd3.begin(), rd3.end(), // input - 1st arg thrust::make_zip_iterator(make_tuple( kpa.begin(), - thrust::make_constant_iterator(opts_init.rc2_T + 273.15) + thrust::make_constant_iterator(opts_init.rc2_T + 273.15), + rd3_insol.begin() )), // input - 2nd arg rc2.begin(), // condition argument rc2.begin(), // output diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index 720a3935..2b1fd2c7 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<2>(tpl) * si::cubic_meters; #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,15 @@ 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 real_t &rd3_insol, 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); return RH - common::kappa_koehler::S_cr( - rd3 * si::cubic_metres, + rd3 * si::cubic_metres, + rd3_insol * si::cubic_metres, kpa, T ); @@ -396,7 +399,8 @@ namespace libcloudphxx thrust::make_permutation_iterator( pimpl->T.begin(), pimpl->ijk.begin() - ) + ), + pimpl->rd3_insol.begin() )), // input - 2nd arg rc2.begin(), // output detail::rw3_cr() // op From 0b74c794df291eb342af55e11810e2ba22c99826 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Fri, 24 Apr 2026 14:21:46 +0200 Subject: [PATCH 08/21] changing rd2_insol to rd3_insol --- include/libcloudph++/common/ice_nucleation.hpp | 18 +++++++++--------- .../particles_impl_fill_outbuf.ipp | 6 +++--- src/impl/ice/particles_impl_ice_nucl_melt.ipp | 10 +++++----- .../particles_impl_init_T_freeze.ipp | 4 ++-- .../particles_impl_init_insol_dry_sizes.ipp | 4 ++-- .../particles_impl_reserve_hskpng_npart.ipp | 2 +- src/impl/particles_impl.ipp | 4 ++-- 7 files changed, 24 insertions(+), 24 deletions(-) diff --git a/include/libcloudph++/common/ice_nucleation.hpp b/include/libcloudph++/common/ice_nucleation.hpp index fec85b0e..e816e8f3 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/src/impl/diagnose_SD_attributes/particles_impl_fill_outbuf.ipp b/src/impl/diagnose_SD_attributes/particles_impl_fill_outbuf.ipp index 82f8e3db..beade301 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 == "rd3_insol" || 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/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_T_freeze.ipp b/src/impl/initialization/particles_impl_init_T_freeze.ipp index ce2097db..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())), - thrust::make_zip_iterator(thrust::make_tuple(rd2_insol.begin() + n_part_old, u01.begin())) + n_part_to_init, + 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_reserve_hskpng_npart.ipp b/src/impl/initialization/particles_impl_reserve_hskpng_npart.ipp index aa2a9d98..97a34d80 100644 --- a/src/impl/initialization/particles_impl_reserve_hskpng_npart.ipp +++ b/src/impl/initialization/particles_impl_reserve_hskpng_npart.ipp @@ -39,7 +39,7 @@ namespace libcloudphxx } if(opts_init.ice_switch) { - rd2_insol.reserve(opts_init.n_sd_max); + rd3_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); diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 681bf689..121ab92a 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 @@ -477,7 +477,7 @@ namespace libcloudphxx if(opts_init.ice_switch) { - distmem_real_vctrs.insert({&rd2_insol, detail::no_initial_value}); + distmem_real_vctrs.insert({&rd3_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}); From c645b355b4f4fe132f395b4de476afabdb1c6e2a Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Fri, 24 Apr 2026 17:01:21 +0200 Subject: [PATCH 09/21] debugging --- .../common/particles_impl_cond_common.ipp | 13 ++++++----- .../percell/particles_impl_cond.ipp | 3 ++- .../perparticle/perparticle_advance_rw2.ipp | 3 ++- ...erparticle_nomixing_adaptive_sstp_cond.ipp | 6 +++-- .../particles_impl_init_wet.ipp | 11 +++++----- src/particles_diag.ipp | 22 ++++++++++--------- 6 files changed, 34 insertions(+), 24 deletions(-) diff --git a/src/impl/condensation/common/particles_impl_cond_common.ipp b/src/impl/condensation/common/particles_impl_cond_common.ipp index 9dd3c734..b03b3dfe 100644 --- a/src/impl/condensation/common/particles_impl_cond_common.ipp +++ b/src/impl/condensation/common/particles_impl_cond_common.ipp @@ -188,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; @@ -223,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 @@ -235,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); } @@ -284,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 4651b429..1983ad2e 100644 --- a/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp +++ b/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp @@ -151,7 +151,8 @@ namespace libcloudphxx kpa, vt, lambda_D, - lambda_K + lambda_K, + rd3_insol ), sstp_tmp_p, RH @@ -226,7 +227,8 @@ namespace libcloudphxx kpa, vt, lambda_D, - lambda_K + lambda_K, + rd3_insol ), sstp_tmp_p, RH diff --git a/src/impl/initialization/particles_impl_init_wet.ipp b/src/impl/initialization/particles_impl_init_wet.ipp index 6e266409..0e8e2129 100644 --- a/src/impl/initialization/particles_impl_init_wet.ipp +++ b/src/impl/initialization/particles_impl_init_wet.ipp @@ -29,10 +29,10 @@ namespace libcloudphxx 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<4>(tpl) * si::cubic_metres; + 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, rd3_insol, kpa, RH, T @@ -63,7 +63,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/particles_diag.ipp b/src/particles_diag.ipp index 2b1fd2c7..7496aa93 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -44,8 +44,8 @@ namespace libcloudphxx 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<2>(tpl) * si::cubic_meters; + const quantity rd3_insol = thrust::get<1>(tpl) * si::cubic_meters; + const quantity T = thrust::get<2>(tpl) * si::kelvins; #if !defined(__NVCC__) using std::pow; @@ -98,15 +98,16 @@ namespace libcloudphxx struct RH_minus_Sc { BOOST_GPU_ENABLED - real_t operator()(const real_t &rd3, const real_t &rd3_insol, 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<1>(tpl) * si::cubic_meters; + 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_insol * si::cubic_metres, + rd3_insol, kpa, T ); @@ -364,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->kpa.begin(), + pimpl->rd3_insol.begin(), thrust::make_permutation_iterator( pimpl->T.begin(), pimpl->ijk.begin() @@ -395,12 +397,12 @@ namespace libcloudphxx thrust::transform( pimpl->rd3.begin(), pimpl->rd3.end(), // input - 1st arg thrust::make_zip_iterator(make_tuple( - pimpl->kpa.begin(), + pimpl->kpa.begin(), + pimpl->rd3_insol.begin(), thrust::make_permutation_iterator( pimpl->T.begin(), pimpl->ijk.begin() - ), - pimpl->rd3_insol.begin() + ) )), // input - 2nd arg rc2.begin(), // output detail::rw3_cr() // op From a41d59e8b3c97229eb3265e21f771fd5d6607788 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Mon, 27 Apr 2026 13:44:27 +0200 Subject: [PATCH 10/21] more debugging --- include/libcloudph++/common/kappa_koehler.hpp | 2 +- .../perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp | 4 +++- src/impl/initialization/particles_impl_init_wet.ipp | 1 + 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/include/libcloudph++/common/kappa_koehler.hpp b/include/libcloudph++/common/kappa_koehler.hpp index 76182287..981b3a5b 100644 --- a/include/libcloudph++/common/kappa_koehler.hpp +++ b/include/libcloudph++/common/kappa_koehler.hpp @@ -149,7 +149,7 @@ namespace libcloudphxx return common::detail::toms748_solve( detail::rw3_eq_minfun(RH, rd3, rd3_insol, 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 + real_t(rw3_eq_nokelvin(rd3, rd3_insol, kappa, RH) / si::cubic_metres) // max ) * si::cubic_metres; } 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 1983ad2e..59a8639a 100644 --- a/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp +++ b/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp @@ -57,6 +57,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; @@ -312,7 +313,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/initialization/particles_impl_init_wet.ipp b/src/impl/initialization/particles_impl_init_wet.ipp index 0e8e2129..ed40d8eb 100644 --- a/src/impl/initialization/particles_impl_init_wet.ipp +++ b/src/impl/initialization/particles_impl_init_wet.ipp @@ -55,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, From 64c270f02f271f863e0a27bb68d0eb09f044ff2b Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Mon, 27 Apr 2026 15:02:34 +0200 Subject: [PATCH 11/21] always initialize rd insol --- src/impl/initialization/particles_impl_init_SD_with_distros.ipp | 2 +- src/impl/initialization/particles_impl_init_SD_with_sizes.ipp | 2 +- src/impl/particles_impl.ipp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) 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/particles_impl.ipp b/src/impl/particles_impl.ipp index 121ab92a..12b81113 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -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({&rd3_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}); From d3ae28a29a668418a45af9bff1d1e61fa09b02a2 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Wed, 29 Apr 2026 15:24:37 +0200 Subject: [PATCH 12/21] debugging wet init --- include/libcloudph++/common/kappa_koehler.hpp | 16 ++++++++++++++-- .../particles_impl_reserve_hskpng_npart.ipp | 2 +- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/include/libcloudph++/common/kappa_koehler.hpp b/include/libcloudph++/common/kappa_koehler.hpp index 981b3a5b..bb3ac40a 100644 --- a/include/libcloudph++/common/kappa_koehler.hpp +++ b/include/libcloudph++/common/kappa_koehler.hpp @@ -52,7 +52,7 @@ namespace libcloudphxx ) { assert(kappa > 0); - assert(rd3_insol < rw3); + assert(rw3 >= (rd3 + rd3_insol)); return (rw3 - rd3 - rd3_insol) / (rw3 - rd3 * (real_t(1) - kappa) - rd3_insol); } @@ -86,6 +86,18 @@ namespace libcloudphxx using std::pow; using std::cbrt; #endif + // std::cerr << "==============================" << std::endl; + // std::cerr << "toms for min : " << real_t((rd3 + rd3_insol) / si::cubic_metres) * real_t(1.000001) << std::endl; + // std::cerr << "max: " << real_t(rw3_eq_nokelvin(rd3, rd3_insol, kappa, RH) / si::cubic_metres) << std::endl; + // std::cerr << "T: " << real_t(T / si::kelvin) << std::endl; + // std::cerr << "RH: " << real_t(RH) << std::endl; + // std::cerr << "Seq: " << real_t(a_w(rw3 * si::cubic_metres, this->rd3, this->rd3_insol, this->kappa) * kelvin::klvntrm(real_t(cbrt(rw3)) * si::metres, this->T)) << std::endl; + // std::cerr << "minfun: " << real_t(RH) - real_t(a_w(rw3 * si::cubic_metres, this->rd3, this->rd3_insol, this->kappa) * kelvin::klvntrm(real_t(cbrt(rw3)) * si::metres, this->T)) << std::endl; + // std::cerr << "rw3: " << real_t(rw3) << std::endl; + // std::cerr << "rd3_insol: " << real_t(rd3_insol / si::cubic_meters) << std::endl; + // std::cerr << "rd3: " << real_t(rd3 / si::cubic_meters) << std::endl; + // std::cerr << "kappa: " << real_t(kappa) << std::endl; + // std::cerr << "Seq no insol: " << real_t(a_w(rw3 * si::cubic_metres, this->rd3, real_t(0) * si::cubic_meters, this->kappa) * kelvin::klvntrm(real_t(cbrt(rw3)) * si::metres, this->T)) << std::endl; return this->RH - a_w(rw3 * si::cubic_metres, this->rd3, this->rd3_insol, this->kappa) * kelvin::klvntrm(real_t(cbrt(rw3)) * si::metres, this->T); @@ -148,7 +160,7 @@ namespace libcloudphxx return common::detail::toms748_solve( detail::rw3_eq_minfun(RH, rd3, rd3_insol, kappa, T), // the above-defined functor - real_t(rd3 / si::cubic_metres), // min + real_t((rd3 + rd3_insol) / si::cubic_metres) * real_t(1.000001), // min real_t(rw3_eq_nokelvin(rd3, rd3_insol, kappa, RH) / si::cubic_metres) // max ) * si::cubic_metres; } diff --git a/src/impl/initialization/particles_impl_reserve_hskpng_npart.ipp b/src/impl/initialization/particles_impl_reserve_hskpng_npart.ipp index 97a34d80..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) { - rd3_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); From 0b4e2046708c2cb10fbc4b96e8cdb7cb785f0bb7 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Wed, 29 Apr 2026 15:41:01 +0200 Subject: [PATCH 13/21] correct limits for toms --- include/libcloudph++/common/kappa_koehler.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/libcloudph++/common/kappa_koehler.hpp b/include/libcloudph++/common/kappa_koehler.hpp index bb3ac40a..6f596d05 100644 --- a/include/libcloudph++/common/kappa_koehler.hpp +++ b/include/libcloudph++/common/kappa_koehler.hpp @@ -160,7 +160,7 @@ namespace libcloudphxx return common::detail::toms748_solve( detail::rw3_eq_minfun(RH, rd3, rd3_insol, kappa, T), // the above-defined functor - real_t((rd3 + rd3_insol) / si::cubic_metres) * real_t(1.000001), // min + 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; } @@ -185,8 +185,8 @@ namespace libcloudphxx return real_t(common::detail::toms748_solve( detail::rw3_cr_minfun(rd3_dbl, rd3_insol_dbl, double(kappa), T_dbl), // the above-defined functor - double(1e0 * (rd3 / si::cubic_metres)), // min - double(1e8 * (rd3 / si::cubic_metres)) // max + double(1e0 * ((rd3 + rd3_insol) / si::cubic_metres)), // min + double(1e8 * ((rd3 + rd3_insol) / si::cubic_metres)) // max )) * si::cubic_metres; } From f77fc2d8284ceede82303f3f2f9c3d531fe9ddeb Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Wed, 13 May 2026 12:11:35 +0200 Subject: [PATCH 14/21] debugging init_wet --- include/libcloudph++/common/kappa_koehler.hpp | 20 +++++-------------- 1 file changed, 5 insertions(+), 15 deletions(-) diff --git a/include/libcloudph++/common/kappa_koehler.hpp b/include/libcloudph++/common/kappa_koehler.hpp index 6f596d05..fd5600a3 100644 --- a/include/libcloudph++/common/kappa_koehler.hpp +++ b/include/libcloudph++/common/kappa_koehler.hpp @@ -52,7 +52,7 @@ namespace libcloudphxx ) { assert(kappa > 0); - assert(rw3 >= (rd3 + rd3_insol)); + 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); } @@ -86,21 +86,11 @@ namespace libcloudphxx using std::pow; using std::cbrt; #endif - // std::cerr << "==============================" << std::endl; - // std::cerr << "toms for min : " << real_t((rd3 + rd3_insol) / si::cubic_metres) * real_t(1.000001) << std::endl; - // std::cerr << "max: " << real_t(rw3_eq_nokelvin(rd3, rd3_insol, kappa, RH) / si::cubic_metres) << std::endl; - // std::cerr << "T: " << real_t(T / si::kelvin) << std::endl; - // std::cerr << "RH: " << real_t(RH) << std::endl; - // std::cerr << "Seq: " << real_t(a_w(rw3 * si::cubic_metres, this->rd3, this->rd3_insol, this->kappa) * kelvin::klvntrm(real_t(cbrt(rw3)) * si::metres, this->T)) << std::endl; - // std::cerr << "minfun: " << real_t(RH) - real_t(a_w(rw3 * si::cubic_metres, this->rd3, this->rd3_insol, this->kappa) * kelvin::klvntrm(real_t(cbrt(rw3)) * si::metres, this->T)) << std::endl; - // std::cerr << "rw3: " << real_t(rw3) << std::endl; - // std::cerr << "rd3_insol: " << real_t(rd3_insol / si::cubic_meters) << std::endl; - // std::cerr << "rd3: " << real_t(rd3 / si::cubic_meters) << std::endl; - // std::cerr << "kappa: " << real_t(kappa) << std::endl; - // std::cerr << "Seq no insol: " << real_t(a_w(rw3 * si::cubic_metres, this->rd3, real_t(0) * si::cubic_meters, this->kappa) * kelvin::klvntrm(real_t(cbrt(rw3)) * si::metres, this->T)) << std::endl; + + real_t rw3_safe = std::max(rw3, real_t((this->rd3 + this->rd3_insol) / si::cubic_metres)); return this->RH - - a_w(rw3 * si::cubic_metres, this->rd3, this->rd3_insol, 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); } }; From 99c867c45ccbc4163d96374e92feb7697c6ae3e9 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Wed, 13 May 2026 12:22:10 +0200 Subject: [PATCH 15/21] code cleanup --- .../libcloudph++/common/ice_deposition.hpp | 136 ------------ .../ice/particles_impl_ice_dep_common.ipp | 200 ------------------ 2 files changed, 336 deletions(-) delete mode 100644 include/libcloudph++/common/ice_deposition.hpp delete mode 100644 src/impl/ice/particles_impl_ice_dep_common.ipp diff --git a/include/libcloudph++/common/ice_deposition.hpp b/include/libcloudph++/common/ice_deposition.hpp deleted file mode 100644 index 549d681d..00000000 --- a/include/libcloudph++/common/ice_deposition.hpp +++ /dev/null @@ -1,136 +0,0 @@ -#pragma once - -#include "units.hpp" -#include "const_cp.hpp" -#include - -#if defined(__NVCC__) -# include -#endif - -namespace libcloudphxx -{ - namespace common - { - namespace ice_deposition - { - - // capacitance of ice crystals - template - BOOST_GPU_ENABLED - quantity ice_capacitance( - const quantity ice_a, // equatorial radius - const quantity ice_c // polar radius - ) - { - if (ice_a < std::numeric_limits::epsilon) throw std::runtime_error("ice_a too small for calculation of capacitance"); - real_t phi = ice_c / ice_a; - real_t e = std::sqrt(real_t(1) - phi * phi); // eccentricity - return phi < real_t(1) ? - ice_a * e / std::asin(e) : ice_a * e / std::log(1 + e) / phi; - } - - - // rate of change of ice mass - template - BOOST_GPU_ENABLED - quantity::type, real_t> dmi_dt( - const quantity D, // D - const quantity K, // K - const quantity rho_v, // ambient water vapour density - const quantity T, // ambient temperature - const quantity p, // ambient pressure - const quantity RH_i, // p_v/p_vsi = relative humidity w.r.t. ice - const quantity ice_a, // ice equatorial radius - const quantity ice_c // ice polar radius - ) - { - using moist_air::R_v; - - quantity::type, real_t> - l_s = const_cp::l_s(T); - - return real_t(4) * pi() * ice_capacitance(ice_a, ice_c) - * (real_t(1) - real_t(1) / RH_i) - / ( - real_t(1) - / D - / rho_v - + - l_s - / K - / RH_i - / T - * (l_s / R_v() / T - real_t(1)) - ) - ; - } - - - // deposition density from Shima et al. (2020) - template - BOOST_GPU_ENABLED - quantity rho_dep( - const quantity D, // D - const quantity K, // K - const quantity rho_v, // ambient water vapour density - const quantity T, // ambient temperature - const quantity RH_i, // p_v/p_vsi = relative humidity w.r.t. ice - const quantity ice_a // ice equatorial radius - ) - { - - return (ice_growth_ratio(T) < real_t(1) && ice_a < real_t(1e-4)*si::meters) ? - moist_air::rho_i<>() : rho_dep_CL94(D, K, rho_v, T, RH_i); - } - - - // deposition density from Chen and Lamb (1994) - template - BOOST_GPU_ENABLED - quantity rho_dep_CL94( - const quantity D, // D - const quantity K, // K - const quantity rho_v, // ambient water vapour density - const quantity T, // ambient temperature - const quantity RH_i // p_v/p_vsi = relative humidity w.r.t. ice - ) - { - using moist_air::R_v; - quantity::type, real_t> - l_s = const_cp::l_s(T); - - using const_cp::p_vs, const_cp::p_vsi; - - quantity d_rhoi = (std::min(RH_i, p_vs/p_vsi) - real_t(1)) - / RH_i - / D - / ( - real_t(1) - / D - / rho_v - + - l_s - / K - / RH_i - / T - * (l_s / R_v() / T - real_t(1)) - ); - - return moist_air::rho_i<>() - * std::exp(- real_t(3) * std::max(d_rhoi/si::kilograms*si::cubic_meters - real_t(5e-5), real_t(0)) - / ice_growth_ratio(T)); - } - - - // inherent ice growth ratio - template - BOOST_GPU_ENABLED - quantity ice_growth_ratio(const quantity T) - { - } - - - }; - }; -}; diff --git a/src/impl/ice/particles_impl_ice_dep_common.ipp b/src/impl/ice/particles_impl_ice_dep_common.ipp deleted file mode 100644 index d51bb322..00000000 --- a/src/impl/ice/particles_impl_ice_dep_common.ipp +++ /dev/null @@ -1,200 +0,0 @@ -// // vim:filetype=cpp -// /** @file -// * @copyright University of Warsaw -// * @section LICENSE -// * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) -// */ -// -// // #include -// #include -// #include -// #include -// #include -// #include -// #include -// #include -// -// namespace libcloudphxx -// { -// namespace lgrngn -// { -// namespace detail -// { -// -// template -// struct advance_ice_a_minfun -// { -// const quantity ice_a_old; -// const quantity dt; -// const quantity rhod; -// const quantity rv; -// const quantity T; -// const quantity p; -// const quantity RH_i; -// const quantity eta; -// const quantity rd3; -// const quantity kpa; -// const quantity vt; -// const quantity RH_max; -// const quantity lambda_D; -// const quantity lambda_K; -// -// // ctor -// BOOST_GPU_ENABLED -// advance_ice_a_minfun( -// const real_t &dt, -// const real_t &ice_a, -// const thrust::tuple, real_t, real_t> &tpl, -// const real_t &RH_max -// ) : -// dt(dt * si::seconds), -// ice_a_old(ice_a * si::square_metres), -// rhod( thrust::get<0>(thrust::get<0>(tpl)) * si::kilograms / si::cubic_metres), -// rv( thrust::get<1>(thrust::get<0>(tpl))), -// T( thrust::get<2>(thrust::get<0>(tpl)) * si::kelvins), -// eta( thrust::get<3>(thrust::get<0>(tpl)) * si::pascals * si::seconds), -// rd3( thrust::get<4>(thrust::get<0>(tpl)) * si::cubic_metres), -// kpa( thrust::get<5>(thrust::get<0>(tpl))), -// vt( thrust::get<6>(thrust::get<0>(tpl)) * si::metres_per_second), -// p( thrust::get<1>(tpl) * si::pascals), -// RH_i( 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), -// RH_max(RH_max) -// {} -// -// BOOST_GPU_ENABLED -// quantity::type, real_t> d_ice_a_dt(const quantity &ice_a) const -// { -// using namespace common::maxwell_mason; -// using namespace common::kappa_koehler; -// using namespace common::kelvin; -// using common::moist_air::D_0; -// using common::moist_air::K_0; -// using common::moist_air::c_pd; -// using common::transition_regime::beta; -// using common::ventil::Sh; -// using common::ventil::Nu; -// #if !defined(__NVCC__) -// using std::sqrt; -// #endif -// -// const quantity -// Re = common::ventil::Re(vt, ice_a, rhod, eta), -// Sc = common::ventil::Sc(eta, rhod, D_0()), // TODO? cache -// Pr = common::ventil::Pr(eta, c_pd(), K_0()); // TODO? cache -// -// const quantity -// D = D_0() * beta(lambda_D / ice_a) * (Sh(Sc, Re) / 2); -// -// const quantity -// K = K_0() * beta(lambda_K / ice_a) * (Nu(Pr, Re) / 2); -// -// return da_dt( -// D, -// K, -// rhod * rv, -// T, -// p, -// RH_i > RH_max ? RH_max : RH_i -// ); -// } -// -// BOOST_GPU_ENABLED -// real_t operator()(const real_t &ice_a_unitless) const -// { -// const quantity ice_a = ice_a_unitless * si::metres; -// return (ice_a_old + dt * d_ice_a_dt(ice_a) - ice_a) / si::metres; -// } -// }; -// -// -// -// template -// struct advance_ice_axis -// { -// const real_t dt, RH_max; -// -// advance_ice_axis(const real_t &dt, const real_t &RH_max) -// : dt(dt), RH_max(RH_max) {} -// -// BOOST_GPU_ENABLED -// thrust::tuple operator()( -// const thrust::tuple &ac_old, -// const thrust::tuple< -// thrust::tuple, -// real_t, real_t> &tpl -// ) const -// { -// #if !defined(__NVCC__) -// using std::max; -// using std::isnan; -// using std::isinf; -// #endif -// const real_t a_old = thrust::get<0>(ac_old); -// const real_t c_old = thrust::get<1>(ac_old); -// -// // Skip liquid droplets -// if (a_old <= real_t(0) || c_old <= real_t(0)) -// return ac_old; -// -// advance_ice_a_minfun f_a(dt, a_old * a_old, tpl, RH_max); -// advance_ice_c_minfun f_c(dt, c_old * c_old, tpl, RH_max); -// -// const real_t da_dt = (f_a.drw2_dt(a_old * a_old * si::square_metres) / (2 * a_old * si::metres)) -// * si::seconds / si::metres; -// const real_t dc_dt = (f_c.drw2_dt(c_old * c_old * si::square_metres) / (2 * c_old * si::metres)) -// * si::seconds / si::metres; -// -// // to store the result -// real_t a_new, c_new; -// -// if (da_dt == 0 && dc_dt == 0) return ac_old; -// -// const real_t rd = cbrt(thrust::get<4>(tpl_in)); -// -// const real_t -// a = max(rd2, rw2_old + min(real_t(0), cond_mlt * drw2)), -// b = rw2_old + max(real_t(0), cond_mlt * drw2); -// -// // numerics (drw2 != 0 but a==b) -// if (a == b) -// { -// if constexpr (apply) -// return rw2_old; -// else -// return real_t(0); -// } -// -// real_t fa, fb; -// -// if (drw2 > 0) -// { -// fa = drw2; // for implicit Euler its equal to min_fun(x_old) -// fb = f(b); -// } -// else -// { -// fa = f(a); -// fb = drw2; // for implicit Euler its equal to min_fun(x_old) -// } -// -// // root-finding ill posed => explicit Euler -// if (fa * fb > 0) rw2_new = rw2_old + drw2; -// // otherwise implicit Euler -// else -// { -// auto _n_iter = n_iter; // we need a copy because toms748_solve expects non-const ref (n_iter is modified by it) -// rw2_new = common::detail::toms748_solve(f, a, b, fa, fb, eps_tolerance, _n_iter); -// } -// -// // check if it doesn't evaporate too much -// if(rw2_new < rd2) rw2_new = rd2; -// -// return thrust::make_tuple(a_new, c_new, rho_new); -// } -// }; -// -// }; -// }; -// }; From 6edfaa705db30d246f150beed715ceadf0728ac3 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Wed, 13 May 2026 13:06:54 +0200 Subject: [PATCH 16/21] cleanup --- models/kinematic_2D/cases/icmw8_case1.hpp | 1 - .../kinematic_2D/src/kin_cloud_2d_lgrngn.hpp | 35 ++----------------- models/kinematic_2D/src/opts_lgrngn.hpp | 25 ++++--------- 3 files changed, 9 insertions(+), 52 deletions(-) diff --git a/models/kinematic_2D/cases/icmw8_case1.hpp b/models/kinematic_2D/cases/icmw8_case1.hpp index 76133a3f..3f44a4f7 100644 --- a/models/kinematic_2D/cases/icmw8_case1.hpp +++ b/models/kinematic_2D/cases/icmw8_case1.hpp @@ -46,7 +46,6 @@ namespace config //aerosol chemical composition parameters (needed for activation) // for lgrngn: quantity kappa; // CCN-derived value from Table 1 in Petters and Kreidenweis 2007 - quantity rd_insol; // insoluble dry radius // for blk_2m: quantity chem_b; //ammonium sulphate //chem_b = 1.33; // sodium chloride // for lagrangian simulations with aq. chemistry diff --git a/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp b/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp index b228fd48..29c96498 100644 --- a/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp +++ b/models/kinematic_2D/src/kin_cloud_2d_lgrngn.hpp @@ -27,7 +27,7 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common // member fields std::unique_ptr> prtcls; - bool coal, sedi, ice_switch, ice_nucl, time_dep_ice_nucl, depo; + bool coal, sedi; // helper methods void diag() @@ -79,32 +79,6 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common } rng_num++; } - // ice_a - int rng_num = 0; - for (auto &rng_moms : params.out_ice) - { - auto &rng(rng_moms.first); - prtcls->diag_ice_a_rng(rng.first / si::metres, rng.second / si::metres); - for (auto &mom : rng_moms.second) - { - prtcls->diag_ice_a_mom(mom); - this->record_aux(aux_name("ice_a", rng_num, mom), prtcls->outbuf()); - } - rng_num++; - } - // ice_c - int rng_num = 0; - for (auto &rng_moms : params.out_ice) - { - auto &rng(rng_moms.first); - prtcls->diag_ice_c_rng(rng.first / si::metres, rng.second / si::metres); - for (auto &mom : rng_moms.second) - { - prtcls->diag_ice_c_mom(mom); - this->record_aux(aux_name("ice_c", rng_num, mom), prtcls->outbuf()); - } - rng_num++; - } } { // rw3(rd) @@ -155,10 +129,6 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common { coal = params.cloudph_opts.coal; sedi = params.cloudph_opts.sedi; - ice_switch = params.cloudph_opts_init.ice_switch; - ice_nucl = params.cloudph_opts.ice_nucl; - time_dep_ice_nucl = params.cloudph_opts_init.time_dep_ice_nucl; - depo = params.cloudph_opts.depo; parent_t::hook_ante_loop(nt); @@ -181,7 +151,6 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common this->record_aux_const("n1_stp", this->setup.n1_stp * si::cubic_metres); this->record_aux_const("n2_stp", this->setup.n2_stp * si::cubic_metres); this->record_aux_const("kappa", this->setup.kappa); - this->record_aux_const("rd_insol", this->setup.rd_insol / si::metres); assert(params.backend != -1); assert(params.dt != 0); @@ -333,7 +302,7 @@ class kin_cloud_2d_lgrngn : public kin_cloud_2d_common bool async = true; libcloudphxx::lgrngn::opts_t cloudph_opts; libcloudphxx::lgrngn::opts_init_t cloudph_opts_init; - outmom_t out_dry, out_wet, out_ice; + outmom_t out_dry, out_wet; }; protected: diff --git a/models/kinematic_2D/src/opts_lgrngn.hpp b/models/kinematic_2D/src/opts_lgrngn.hpp index d6635168..924921d2 100644 --- a/models/kinematic_2D/src/opts_lgrngn.hpp +++ b/models/kinematic_2D/src/opts_lgrngn.hpp @@ -76,8 +76,7 @@ void parse_moms( std::vector> min_maxnum; outmom_t &moms = opt == "out_dry" ? rt_params.out_dry : - opt == "out_wet" ? rt_params.out_wet : - rt_params.out_ice; + rt_params.out_wet; const bool result = qi::phrase_parse(first, last, *( @@ -141,10 +140,9 @@ void parse_moms( outmom_t &moms = opt == "out_dry" ? rt_params.out_dry : - opt == "out_wet" ? rt_params.out_wet : - opt == "out_ice" ? rt_params.out_ice : - opt == "out_chem" ? rt_params.out_chem : - rt_params.out_wet_pH; + opt == "out_wet" ? rt_params.out_wet : + opt == "out_chem" ? rt_params.out_chem : + rt_params.out_wet_pH; const bool result = qi::phrase_parse(first, last, *( @@ -200,7 +198,7 @@ void setopts_micro_chem( {solver_t::ix::th, {"th", "[K]"}}, {solver_t::ix::rv, {"rv", "[kg kg-1]"}} }; - out_set = {"out_dry", "out_wet", "out_ice"}; + out_set = {"out_dry", "out_wet"}; } template @@ -219,7 +217,7 @@ void setopts_micro_chem( {solver_t::ix::NH3g, {"NH3g", "[dimesnionless]"}}, {solver_t::ix::HNO3g, {"HNO3g","[dimensionless]"}} }; - out_set = {"out_dry", "out_wet", "out_ice", "out_chem", "out_wet_pH"}; + out_set = {"out_dry", "out_wet", "out_chem", "out_wet_pH"}; } // simulation and output parameters for micro=lgrngn @@ -248,10 +246,6 @@ void setopts_micro( ("chem_dsc", po::value()->default_value(rt_params.cloudph_opts.chem_dsc) , "dissociation (1=on, 0=off)") ("chem_rct", po::value()->default_value(rt_params.cloudph_opts.chem_rct) , "aqueous chemistry (1=on, 0=off)") ("chem_switch", po::value()->default_value(rt_params.cloudph_opts_init.chem_switch) , "aqueous chemistry (1=on, 0=off)") - ("ice_switch", po::value()->default_value(rt_params.cloudph_opts_init.ice_switch) , "enable ice (1=on, 0=off)") - ("ice_nucl", po::value()->default_value(rt_params.cloudph_opts.ice_nucl) , "ice nucleation (1=on, 0=off)") - ("time_dep_ice_nucl", po::value()->default_value(rt_params.cloudph_opts_init.time_dep_ice_nucl) , "time dependent ice nucleation (1=on, 0=off)") - ("depo", po::value()->default_value(rt_params.cloudph_opts.depo) , "ice depositional growth (1=on, 0=off)") // free parameters ("sstp_cond", po::value()->default_value(rt_params.cloudph_opts_init.sstp_cond), "no. of substeps for condensation") ("sstp_coal", po::value()->default_value(rt_params.cloudph_opts_init.sstp_coal), "no. of substeps for coalescence") @@ -261,7 +255,6 @@ void setopts_micro( // output ("out_dry", po::value()->default_value("0:1|0"), "dry radius ranges and moment numbers (r1:r2|n1,n2...;...)") ("out_wet", po::value()->default_value(".5e-6:25e-6|0,1,2,3;25e-6:1|0,3,6"), "wet radius ranges and moment numbers (r1:r2|n1,n2...;...)") - ("out_ice", po::value()->default_value(".1e-6:1|0,1,2,3"), "ice semi-axis ranges and moment numbers (r1:r2|n1,n2...;...)") ("out_wet_pH", po::value()->default_value("0:1|0"), "wet radius ranges for output of H+ and S_VI)") ("out_chem", po::value()->default_value("0:1|0"), "dry radius ranges for which chem mass is outputted") // collision and sedimentation @@ -304,7 +297,7 @@ void setopts_micro( */ rt_params.cloudph_opts_init.dry_distros.emplace( - libcloudphxx::lgrngn::kappa_rd_insol_t{setup.kappa, setup.rd_insol / si::meters}, // kappa, rd_insol + libcloudphxx::lgrngn::kappa_rd_insol_t{setup.kappa, config::real_t(0.)}, // kappa, rd_insol std::make_shared> (setup) ); @@ -313,10 +306,6 @@ void setopts_micro( rt_params.cloudph_opts.sedi = vm["sedi"].as(); rt_params.cloudph_opts.cond = vm["cond"].as(); rt_params.cloudph_opts.coal = vm["coal"].as(); - rt_params.cloudph_opts_init.ice_switch = vm["ice_switch"].as(); - rt_params.cloudph_opts.ice_nucl = vm["ice_nucl"].as(); - rt_params.cloudph_opts_init.time_dep_ice_nucl = vm["time_dep_ice_nucl"].as(); - rt_params.cloudph_opts.depo = vm["depo"].as(); rt_params.cloudph_opts.rcyc = vm["rcyc"].as(); rt_params.cloudph_opts.chem_dsl = vm["chem_dsl"].as(); From b8c340bf1312551c63b105c71f0b5d1de33fc790 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Wed, 13 May 2026 17:56:20 +0200 Subject: [PATCH 17/21] debugging unit tests --- .../particles_impl_update_incloud_time.ipp | 6 +++--- tests/python/unit/api_common.py | 4 ++-- tests/python/unit/diag_incloud_time.py | 6 +++++- tests/python/unit/source.py | 2 +- 4 files changed, 11 insertions(+), 7 deletions(-) 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 5a3eec8d..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,12 +43,12 @@ 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() - ), - rd3_insol.begin() + ) )), // input - 2nd arg rc2.begin(), // output detail::rw3_cr() // op 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..71c8512c 100644 --- a/tests/python/unit/diag_incloud_time.py +++ b/tests/python/unit/diag_incloud_time.py @@ -19,8 +19,12 @@ 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.RH_max = 0.999 # to comply with the assert(RH<1) at init 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()) From 95f08b02f68fdd6c6bbb4853c68d1680f454f898 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Thu, 14 May 2026 10:46:05 +0200 Subject: [PATCH 18/21] fix rw3_cr arguments --- src/impl/housekeeping/particles_impl_hskpng_rc2.ipp | 4 ++-- src/particles_diag.ipp | 12 ++++++------ tests/python/unit/diag_incloud_time.py | 1 + 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp b/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp index 02a36727..6ad3b1a8 100644 --- a/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp +++ b/src/impl/housekeeping/particles_impl_hskpng_rc2.ipp @@ -21,9 +21,9 @@ 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), - rd3_insol.begin() + thrust::make_constant_iterator(opts_init.rc2_T + 273.15) )), // input - 2nd arg rc2.begin(), // condition argument rc2.begin(), // output diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index 7496aa93..5c95c676 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -43,8 +43,8 @@ namespace libcloudphxx BOOST_GPU_ENABLED real_t operator()(const real_t &rd3, const thrust::tuple &tpl) { - const quantity kpa = thrust::get<0>(tpl); - const quantity rd3_insol = thrust::get<1>(tpl) * si::cubic_meters; + 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__) @@ -100,8 +100,8 @@ namespace libcloudphxx BOOST_GPU_ENABLED real_t operator()(const real_t &rd3, const thrust::tuple &tpl) { - const quantity kpa = thrust::get<0>(tpl); - const quantity rd3_insol = thrust::get<1>(tpl) * si::cubic_meters; + 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); @@ -365,8 +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() @@ -397,8 +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/diag_incloud_time.py b/tests/python/unit/diag_incloud_time.py index 71c8512c..de78274b 100644 --- a/tests/python/unit/diag_incloud_time.py +++ b/tests/python/unit/diag_incloud_time.py @@ -27,6 +27,7 @@ def lognormal(lnr): 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) From 59e7c621b027d4dbed5cf6c90ecdf96f75e98219 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Fri, 15 May 2026 10:46:42 +0200 Subject: [PATCH 19/21] suggestions from copilot --- include/libcloudph++/common/kappa_koehler.hpp | 7 +++++-- .../diagnose_SD_attributes/particles_impl_fill_outbuf.ipp | 2 +- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/include/libcloudph++/common/kappa_koehler.hpp b/include/libcloudph++/common/kappa_koehler.hpp index fd5600a3..4896e096 100644 --- a/include/libcloudph++/common/kappa_koehler.hpp +++ b/include/libcloudph++/common/kappa_koehler.hpp @@ -85,9 +85,12 @@ namespace libcloudphxx #if !defined(__NVCC__) using std::pow; using std::cbrt; + using std::max; +#else + using thrust::max; #endif - real_t rw3_safe = std::max(rw3, real_t((this->rd3 + this->rd3_insol) / si::cubic_metres)); + real_t rw3_safe = max(rw3, real_t((this->rd3 + this->rd3_insol) / si::cubic_metres)); return this->RH - 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); @@ -146,7 +149,7 @@ 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, rd3_insol, kappa, T), // the above-defined functor 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 beade301..ff86b1cd 100644 --- a/src/impl/diagnose_SD_attributes/particles_impl_fill_outbuf.ipp +++ b/src/impl/diagnose_SD_attributes/particles_impl_fill_outbuf.ipp @@ -44,7 +44,7 @@ namespace libcloudphxx 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 == "rd3_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."); } From f8977fbc2320395b230e9ff391b34e5c2f4fb878 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Fri, 15 May 2026 11:10:41 +0200 Subject: [PATCH 20/21] fix for cuda builds --- include/libcloudph++/common/kappa_koehler.hpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/include/libcloudph++/common/kappa_koehler.hpp b/include/libcloudph++/common/kappa_koehler.hpp index 4896e096..376de243 100644 --- a/include/libcloudph++/common/kappa_koehler.hpp +++ b/include/libcloudph++/common/kappa_koehler.hpp @@ -52,7 +52,12 @@ namespace libcloudphxx ) { assert(kappa > 0); - assert(rw3 >= rd3 + rd3_insol - std::numeric_limits::epsilon() * si::cubic_metres); +#if !defined(__NVCC__) + const real_t eps = std::numeric_limits::epsilon(); +#else + const real_t eps = thrust::numeric_limits::epsilon(); +#endif + assert(rw3 >= rd3 + rd3_insol - eps * si::cubic_metres); return (rw3 - rd3 - rd3_insol) / (rw3 - rd3 * (real_t(1) - kappa) - rd3_insol); } From 6d9d52d86e44c2ea30ffaea14fc60d429d3fb9bd Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Fri, 15 May 2026 12:05:18 +0200 Subject: [PATCH 21/21] another fix for cuda builds --- include/libcloudph++/common/kappa_koehler.hpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/include/libcloudph++/common/kappa_koehler.hpp b/include/libcloudph++/common/kappa_koehler.hpp index 376de243..4896e096 100644 --- a/include/libcloudph++/common/kappa_koehler.hpp +++ b/include/libcloudph++/common/kappa_koehler.hpp @@ -52,12 +52,7 @@ namespace libcloudphxx ) { assert(kappa > 0); -#if !defined(__NVCC__) - const real_t eps = std::numeric_limits::epsilon(); -#else - const real_t eps = thrust::numeric_limits::epsilon(); -#endif - assert(rw3 >= rd3 + rd3_insol - eps * si::cubic_metres); + 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); }