Skip to content

Insoluble dry radius in S_eq formulas #488

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
AgnieszkaMakulska wants to merge 24 commits into
base: master
Choose a base branch
from
Open

Insoluble dry radius in S_eq formulas #488

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
213f53e
ice growth formulas
AgnieszkaMakulska Feb 3, 2026
ae1d70d
dep density formulas
AgnieszkaMakulska Feb 4, 2026
eb99094
Merge branch 'master' into ice
AgnieszkaMakulska Apr 8, 2026
fa2aff0
kinematic_2D ice options
AgnieszkaMakulska Apr 17, 2026
ed2331c
fix u01 for T_freeze
AgnieszkaMakulska Apr 17, 2026
0e3a860
rd_insol in Kohler functions
AgnieszkaMakulska Apr 22, 2026
c1c421d
init wet with rd_insol
AgnieszkaMakulska Apr 23, 2026
6ef77f8
update calls to functions with rd_insol
AgnieszkaMakulska Apr 23, 2026
0b74c79
changing rd2_insol to rd3_insol
AgnieszkaMakulska Apr 24, 2026
c645b35
debugging
AgnieszkaMakulska Apr 24, 2026
a41d59e
more debugging
AgnieszkaMakulska Apr 27, 2026
64c270f
always initialize rd insol
AgnieszkaMakulska Apr 27, 2026
d3ae28a
debugging wet init
AgnieszkaMakulska Apr 29, 2026
0b4e204
correct limits for toms
AgnieszkaMakulska Apr 29, 2026
f77fc2d
debugging init_wet
AgnieszkaMakulska May 13, 2026
99c867c
code cleanup
AgnieszkaMakulska May 13, 2026
6edfaa7
cleanup
AgnieszkaMakulska May 13, 2026
b8c340b
debugging unit tests
AgnieszkaMakulska May 13, 2026
95f08b0
fix rw3_cr arguments
AgnieszkaMakulska May 14, 2026
ad6f8b8
Merge branch 'master' into ice
pdziekan May 14, 2026
59e7c62
suggestions from copilot
AgnieszkaMakulska May 15, 2026
9ddf421
Merge branch 'ice' of github.com:AgnieszkaMakulska/libcloudphxx into ice
AgnieszkaMakulska May 15, 2026
f8977fb
fix for cuda builds
AgnieszkaMakulska May 15, 2026
6d9d52d
another fix for cuda builds
AgnieszkaMakulska May 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions bindings/python/common.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,20 +119,22 @@ namespace libcloudphxx
}

template <typename real_t>
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 <typename real_t>
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
);
Comment thread
AgnieszkaMakulska marked this conversation as resolved.
Expand Down
18 changes: 9 additions & 9 deletions include/libcloudph++/common/ice_nucleation.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ namespace libcloudphxx
BOOST_GPU_ENABLED
quantity<si::temperature, real_t> 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)
Expand All @@ -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))
{
Expand All @@ -54,12 +54,12 @@ namespace libcloudphxx
BOOST_GPU_ENABLED
real_t operator()(const thrust::tuple<real_t, real_t> &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<real_t>(
INP_type,
rd2_insol,
rd3_insol,
rand
) / si::kelvin;
}
Expand All @@ -71,21 +71,21 @@ 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__)
* pi<real_t>()
#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<real_t>(T * si::kelvin)/ const_cp::p_vs<real_t>(T * si::kelvin); // water activity
if (INP_type == INP_t::mineral)
{
Expand Down Expand Up @@ -127,13 +127,13 @@ namespace libcloudphxx
BOOST_GPU_ENABLED
real_t operator()(const thrust::tuple<real_t, real_t, real_t> &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<real_t>(
INP_type,
rd2_insol,
rd3_insol,
rw2,
T,
dt
Expand Down
60 changes: 38 additions & 22 deletions include/libcloudph++/common/kappa_koehler.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,15 +29,16 @@ namespace libcloudphxx
template <typename real_t>
BOOST_GPU_ENABLED
quantity<si::volume, real_t> rw3_eq_nokelvin(
quantity<si::volume, real_t> rd3,
quantity<si::volume, real_t> rd3,
quantity<si::volume, real_t> rd3_insol,
quantity<si::dimensionless, real_t> kappa,
quantity<si::dimensionless, real_t> 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
Expand All @@ -46,11 +47,13 @@ namespace libcloudphxx
quantity<si::dimensionless, real_t> a_w(
quantity<si::volume, real_t> rw3,
quantity<si::volume, real_t> rd3,
quantity<si::volume, real_t> rd3_insol,
quantity<si::dimensionless, real_t> kappa
)
{
assert(kappa > 0);
return (rw3 - rd3) / (rw3 - rd3 * (real_t(1) - kappa));
assert(rw3 >= rd3 + rd3_insol - std::numeric_limits<real_t>::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
Expand All @@ -62,6 +65,7 @@ namespace libcloudphxx
{
const quantity<si::dimensionless, real_t> RH;
const quantity<si::volume, real_t> rd3;
const quantity<si::volume, real_t> rd3_insol;
const quantity<si::dimensionless, real_t> kappa;
const quantity<si::temperature, real_t> T;

Expand All @@ -70,20 +74,26 @@ namespace libcloudphxx
rw3_eq_minfun(
const quantity<si::dimensionless, real_t> &RH,
const quantity<si::volume, real_t> &rd3,
const quantity<si::volume, real_t> &rd3_insol,
const quantity<si::dimensionless, real_t> &kappa,
const quantity<si::temperature, real_t> &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
{
#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);
Comment thread
AgnieszkaMakulska marked this conversation as resolved.
}
};

Expand All @@ -92,16 +102,18 @@ namespace libcloudphxx
struct rw3_cr_minfun
{
const quantity<si::volume, real_t> rd3;
const quantity<si::volume, real_t> rd3_insol;
const quantity<si::dimensionless, real_t> kappa;
const quantity<si::temperature, real_t> T;

// ctor
BOOST_GPU_ENABLED
rw3_cr_minfun(
const quantity<si::volume, real_t> &rd3,
const quantity<si::volume, real_t> &rd3_insol,
const quantity<si::dimensionless, real_t> &kappa,
const quantity<si::temperature, real_t> &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
Expand All @@ -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;
}
Expand All @@ -127,7 +139,8 @@ namespace libcloudphxx
template <typename real_t>
BOOST_GPU_ENABLED
quantity<si::volume, real_t> rw3_eq(
quantity<si::volume, real_t> rd3,
quantity<si::volume, real_t> rd3,
quantity<si::volume, real_t> rd3_insol,
quantity<si::dimensionless, real_t> kappa,
quantity<si::dimensionless, real_t> RH,
quantity<si::temperature, real_t> T
Expand All @@ -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<real_t>(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<real_t>(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;
}

Expand All @@ -151,20 +164,22 @@ namespace libcloudphxx
template <typename real_t>
BOOST_GPU_ENABLED
quantity<si::volume, real_t> rw3_cr(
quantity<si::volume, real_t> rd3,
quantity<si::volume, real_t> rd3,
quantity<si::volume, real_t> rd3_insol,
quantity<si::dimensionless, real_t> kappa,
quantity<si::temperature, real_t> T
)
{
assert(kappa > 0);

quantity<si::volume, double> rd3_dbl = static_cast<quantity<si::volume, double> >(rd3);
quantity<si::volume, double> rd3_insol_dbl = static_cast<quantity<si::volume, double> >(rd3_insol);
quantity<si::temperature, double> T_dbl = static_cast<quantity<si::temperature, double> >(T);

return real_t(common::detail::toms748_solve(
detail::rw3_cr_minfun<double>(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<double>(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;
}

Expand All @@ -173,6 +188,7 @@ namespace libcloudphxx
BOOST_GPU_ENABLED
quantity<si::dimensionless, real_t> S_cr(
quantity<si::volume, real_t> rd3,
quantity<si::volume, real_t> rd3_insol,
quantity<si::dimensionless, real_t> kappa,
quantity<si::temperature, real_t> T
)
Expand All @@ -182,9 +198,9 @@ namespace libcloudphxx
using std::pow;
using std::cbrt;
#endif
quantity<si::volume, real_t> rw3 = rw3_cr(rd3, kappa, T);
quantity<si::volume, real_t> 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
);
Expand Down
19 changes: 12 additions & 7 deletions src/impl/condensation/common/particles_impl_cond_common.ipp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,13 +93,14 @@ namespace libcloudphxx
const quantity<si::dimensionless, real_t> RH_max;
const quantity<si::length, real_t> lambda_D;
const quantity<si::length, real_t> lambda_K;
const quantity<si::volume, real_t> rd3_insol;

// ctor
BOOST_GPU_ENABLED
advance_rw2_minfun(
const real_t &dt,
const real_t &rw2,
const thrust::tuple<thrust::tuple<real_t, real_t, real_t, real_t, real_t, real_t, real_t, real_t, real_t>, real_t, real_t> &tpl,
const thrust::tuple<thrust::tuple<real_t, real_t, real_t, real_t, real_t, real_t, real_t, real_t, real_t, real_t>, real_t, real_t> &tpl,
const real_t &RH_max
) :
dt(dt * si::seconds),
Expand All @@ -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)
{}

Expand Down Expand Up @@ -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)
);
}
Expand Down Expand Up @@ -186,7 +188,7 @@ namespace libcloudphxx
BOOST_GPU_ENABLED
real_t operator()(
const real_t &rw2_old,
const thrust::tuple<thrust::tuple<real_t, real_t, real_t, real_t, real_t, real_t, real_t, real_t, real_t>, real_t, real_t> &tpl
const thrust::tuple<thrust::tuple<real_t, real_t, real_t, real_t, real_t, real_t, real_t, real_t, real_t, real_t>, real_t, real_t> &tpl
) const {
#if !defined(__NVCC__)
using std::min;
Expand Down Expand Up @@ -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
Expand All @@ -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);
}
Expand Down Expand Up @@ -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);
}
Expand Down
3 changes: 2 additions & 1 deletion src/impl/condensation/percell/particles_impl_cond.ipp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion src/impl/condensation/perparticle/perparticle_advance_rw2.ipp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
Loading
Loading