From fec305b547653c94db3e4449b4c5e2b036899abc Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Tue, 25 Aug 2026 15:19:25 +0200 Subject: [PATCH 1/4] T_freeze cant be lower than homogeneous --- include/libcloudph++/common/ice_nucleation.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/libcloudph++/common/ice_nucleation.hpp b/include/libcloudph++/common/ice_nucleation.hpp index 58340c84..76cab019 100644 --- a/include/libcloudph++/common/ice_nucleation.hpp +++ b/include/libcloudph++/common/ice_nucleation.hpp @@ -2,6 +2,7 @@ #include "units.hpp" #include "const_cp.hpp" +#include #include #if defined(__NVCC__) @@ -27,14 +28,14 @@ namespace libcloudphxx real_t A = real_t(4) #if !defined(__NVCC__) * pi() - #else + #else * CUDART_PI #endif * rd2_insol; // surface area of the insoluble particle if (INP_type == INP_t::mineral && A > real_t(1e-20)) { - return real_t(real_t(273.15) + (real_t(8.934) - log(- log(real_t(1.) - rand) / A) ) / real_t(0.517)) * si::kelvin; + return std::max(real_t(235.15), real_t(real_t(273.15) + (real_t(8.934) - log(- log(real_t(1.) - rand) / A) ) / real_t(0.517))) * si::kelvin; } else { From 4e824a5b7276194f23c43ad396a5dc69e34012b6 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 26 Aug 2026 16:24:02 +0200 Subject: [PATCH 2/4] implement T_freeze for silver iodide --- .../libcloudph++/common/ice_nucleation.hpp | 68 ++++++++++++++----- 1 file changed, 51 insertions(+), 17 deletions(-) diff --git a/include/libcloudph++/common/ice_nucleation.hpp b/include/libcloudph++/common/ice_nucleation.hpp index 76cab019..544120bd 100644 --- a/include/libcloudph++/common/ice_nucleation.hpp +++ b/include/libcloudph++/common/ice_nucleation.hpp @@ -15,31 +15,57 @@ namespace libcloudphxx { namespace ice_nucleation { - enum class INP_t {mineral}; // types of ice nucleating particles, TODO: add more types + enum class INP_t {mineral, AgI}; // types of ice nucleating particles, TODO: add more types // Inverse CDF for singular freezing temperature as defined in eq. 1 in Shima et al., 2020 + // frozen_fraction(T) = p(T_f > T) + // CDF of T_f = 1 - frozen_fraction(T) template BOOST_GPU_ENABLED quantity T_freeze_CDF_inv( - const INP_t& INP_type, // type of ice nucleating particle + 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 rand // random number between [0, 1] ) { - real_t A = real_t(4) - #if !defined(__NVCC__) - * pi() - #else - * CUDART_PI - #endif - * rd2_insol; // surface area of the insoluble particle - - if (INP_type == INP_t::mineral && A > real_t(1e-20)) - { - return std::max(real_t(235.15), real_t(real_t(273.15) + (real_t(8.934) - log(- log(real_t(1.) - rand) / A) ) / real_t(0.517))) * si::kelvin; - } - else + static constexpr quantity T_homo = real_t(235.15) * si::kelvin; + // static constexpr quantity Niemand_T_min = real_t(273.15 - 36) * si::kelvin; + // static constexpr quantity Niemand_T_max = real_t(273.15 - 12) * si::kelvin; + static constexpr quantity Omanovic_b = 0.97; + static constexpr quantity Omanovic_k = 0.88; + static constexpr quantity Omanovic_T0 = real_t(263.95) * si::kelvin; + + if(rd2_insol < 1e-20) return T_homo; + + switch(INP_type) { - return real_t(235.15) * si::kelvin; // the default freezing temperature is -38 C + case INP_t::mineral: { + // active site ns(T) parameterization from Niemand et al. 2012 for mineral dust + // Shima et al. 2020 (and many others): p(T_f > T) = 1 - exp(A * ns(T)) + // Niemand et al. 2012 for mineral dust: ns(T) = exp(-0.517(T - 273.15) + 8.934) [m^-2] + // NOTE (not used by us): Shima et al. 2020: use Niemand only in the range -36 C to -12 C, above ns(T)=0 and below ns(T)=ns(-36) + + const real_t A = real_t(4) + #if !defined(__NVCC__) + * pi() + #else + * CUDART_PI + #endif + * rd2_insol; // surface area of the insoluble particle + + const real_t Niemand_T_freeze = real_t(real_t(273.15) + (real_t(8.934) - log(- log(real_t(1.) - rand) / (real_t(4) * pi() * rd2_insol)) ) / real_t(0.517)); + return std::max(real_t(T_homo / si::kelvin), Niemand_T_freeze) * si::kelvin; + break; + } + case INP_t::AgI: { + // Omanovic et al. 2024 for AgI: frozen_fraction = b [ 1 - 1 / (1 + exp(-k(T-T0)))]; b=0.97, k=0.88, T0=263.95K ; applicable for rd_insol>20nm + // leads to T_f = T0 - 1/k ln((1-R)/(b-1+R)) , where R is uniformly distributed [0,1]; as T->0 FF=b, so if R>b, set Tf=0? + if(rand > Omanovic_b) return T_homo; + const real_t Omanovic_T_freeze = real_t(Omanovic_T0/ si::kelvin) - (real_t(1.) / Omanovic_k) * log((real_t(1.) - rand) / (Omanovic_b - real_t(1.) + rand)); + return std::max(real_t(T_homo / si::kelvin), Omanovic_T_freeze) * si::kelvin; + break; + } + default: + throw std::runtime_error("Unrecognized INP type"); } } @@ -93,8 +119,16 @@ namespace libcloudphxx real_t J_het = pow(real_t(10), real_t(-1.35) + real_t(22.62) * d_aw) * real_t(1e4); // nucleation rate return 1 - exp(- J_het * A * dt); } + else if (INP_type == INP_t::AgI) + { + throw std::runtime_error("AgI time-dependent freezing not implemented yet"); + return 0; + } else - return real_t(0.); // TODO: other INP types + { + throw std::runtime_error("Unrecognized INP type"); + return 0; + } } else { From 505bd0cc20b8b64b3f9cae0c0ff414c0c71ce322 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 26 Aug 2026 16:30:43 +0200 Subject: [PATCH 3/4] AgI freezing: parameterization works only for large enough particles --- include/libcloudph++/common/ice_nucleation.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/libcloudph++/common/ice_nucleation.hpp b/include/libcloudph++/common/ice_nucleation.hpp index 544120bd..0022dad2 100644 --- a/include/libcloudph++/common/ice_nucleation.hpp +++ b/include/libcloudph++/common/ice_nucleation.hpp @@ -57,9 +57,10 @@ namespace libcloudphxx break; } case INP_t::AgI: { - // Omanovic et al. 2024 for AgI: frozen_fraction = b [ 1 - 1 / (1 + exp(-k(T-T0)))]; b=0.97, k=0.88, T0=263.95K ; applicable for rd_insol>20nm + // Omanovic et al. 2024 for AgI: frozen_fraction = b [ 1 - 1 / (1 + exp(-k(T-T0)))]; b=0.97, k=0.88, T0=263.95K // leads to T_f = T0 - 1/k ln((1-R)/(b-1+R)) , where R is uniformly distributed [0,1]; as T->0 FF=b, so if R>b, set Tf=0? - if(rand > Omanovic_b) return T_homo; + assert(rd2_insol > 4e-16); // this size-agnostic parameterisation is applicable for rd_insol>20nm + if(rand > Omanovic_b || rd2_insol <= 4e-16) return T_homo; // in non-debug runs, ignore small AgI INPs const real_t Omanovic_T_freeze = real_t(Omanovic_T0/ si::kelvin) - (real_t(1.) / Omanovic_k) * log((real_t(1.) - rand) / (Omanovic_b - real_t(1.) + rand)); return std::max(real_t(T_homo / si::kelvin), Omanovic_T_freeze) * si::kelvin; break; From 275f52769935b5fcd9588770253a5ca1a9aee154 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 27 Aug 2026 18:20:35 +0200 Subject: [PATCH 4/4] T_freeze fixes for CUDA --- .../libcloudph++/common/ice_nucleation.hpp | 48 +++++++++++++------ .../particles_impl_fill_outbuf.ipp | 2 +- 2 files changed, 34 insertions(+), 16 deletions(-) diff --git a/include/libcloudph++/common/ice_nucleation.hpp b/include/libcloudph++/common/ice_nucleation.hpp index 0022dad2..ca826115 100644 --- a/include/libcloudph++/common/ice_nucleation.hpp +++ b/include/libcloudph++/common/ice_nucleation.hpp @@ -3,6 +3,7 @@ #include "units.hpp" #include "const_cp.hpp" #include +#include #include #if defined(__NVCC__) @@ -25,16 +26,25 @@ namespace libcloudphxx 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 rand // random number between [0, 1] + const real_t rand // random number between [0, 1] (edge values are guarded inside the function) ) { - static constexpr quantity T_homo = real_t(235.15) * si::kelvin; - // static constexpr quantity Niemand_T_min = real_t(273.15 - 36) * si::kelvin; - // static constexpr quantity Niemand_T_max = real_t(273.15 - 12) * si::kelvin; - static constexpr quantity Omanovic_b = 0.97; - static constexpr quantity Omanovic_k = 0.88; - static constexpr quantity Omanovic_T0 = real_t(263.95) * si::kelvin; - if(rd2_insol < 1e-20) return T_homo; + #if !defined(__NVCC__) + using std::min; + using std::max; + #endif + + static constexpr real_t T_homo = 235.15; + static constexpr real_t Niemand_T_max = 261.15; + static constexpr real_t Omanovic_b = 0.97; + static constexpr real_t Omanovic_k = 0.88; + static constexpr real_t Omanovic_T0 = 263.95; + + if(rd2_insol < 1e-20) return T_homo * si::kelvin; + + // to avoid rand==0 (breaks mineral dust) and rand==1 (breaks AgI) + const real_t eps = 1e-10;//thrust::numeric_limits::epsilon(); + const real_t rand_guarded = max(eps, min(rand, real_t(1) - eps)); switch(INP_type) { @@ -42,7 +52,7 @@ namespace libcloudphxx // active site ns(T) parameterization from Niemand et al. 2012 for mineral dust // Shima et al. 2020 (and many others): p(T_f > T) = 1 - exp(A * ns(T)) // Niemand et al. 2012 for mineral dust: ns(T) = exp(-0.517(T - 273.15) + 8.934) [m^-2] - // NOTE (not used by us): Shima et al. 2020: use Niemand only in the range -36 C to -12 C, above ns(T)=0 and below ns(T)=ns(-36) + // following Shima et al. 2020: use Niemand only in the range -36 C to -12 C, above ns(T)=0 and below ns(T)=ns(-36) const real_t A = real_t(4) #if !defined(__NVCC__) @@ -52,21 +62,25 @@ namespace libcloudphxx #endif * rd2_insol; // surface area of the insoluble particle - const real_t Niemand_T_freeze = real_t(real_t(273.15) + (real_t(8.934) - log(- log(real_t(1.) - rand) / (real_t(4) * pi() * rd2_insol)) ) / real_t(0.517)); - return std::max(real_t(T_homo / si::kelvin), Niemand_T_freeze) * si::kelvin; + const real_t Niemand_T_freeze = real_t(real_t(273.15) + (real_t(8.934) - log(- log(real_t(1.) - rand_guarded) / (real_t(4) * pi() * rd2_insol)) ) / real_t(0.517)); + return min(Niemand_T_max, max(T_homo, Niemand_T_freeze)) * si::kelvin; break; } case INP_t::AgI: { // Omanovic et al. 2024 for AgI: frozen_fraction = b [ 1 - 1 / (1 + exp(-k(T-T0)))]; b=0.97, k=0.88, T0=263.95K // leads to T_f = T0 - 1/k ln((1-R)/(b-1+R)) , where R is uniformly distributed [0,1]; as T->0 FF=b, so if R>b, set Tf=0? assert(rd2_insol > 4e-16); // this size-agnostic parameterisation is applicable for rd_insol>20nm - if(rand > Omanovic_b || rd2_insol <= 4e-16) return T_homo; // in non-debug runs, ignore small AgI INPs - const real_t Omanovic_T_freeze = real_t(Omanovic_T0/ si::kelvin) - (real_t(1.) / Omanovic_k) * log((real_t(1.) - rand) / (Omanovic_b - real_t(1.) + rand)); - return std::max(real_t(T_homo / si::kelvin), Omanovic_T_freeze) * si::kelvin; + if(rand_guarded < (1-Omanovic_b) || rd2_insol <= 4e-16) return T_homo * si::kelvin; // b=0.97, so 3% dont freeze heterogeneously; in non-debug runs, ignore small AgI INPs + const real_t Omanovic_T_freeze = real_t(Omanovic_T0) - (real_t(1.) / Omanovic_k) * log((real_t(1.) - rand_guarded) / (Omanovic_b - real_t(1.) + rand_guarded)); + return max(T_homo, Omanovic_T_freeze) * si::kelvin; break; } - default: + default: { +#if !defined(__NVCC__) throw std::runtime_error("Unrecognized INP type"); +#endif + return T_homo * si::kelvin; + } } } @@ -122,12 +136,16 @@ namespace libcloudphxx } else if (INP_type == INP_t::AgI) { +#if !defined(__NVCC__) throw std::runtime_error("AgI time-dependent freezing not implemented yet"); +#endif return 0; } else { +#if !defined(__NVCC__) throw std::runtime_error("Unrecognized INP type"); +#endif return 0; } } 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 0908c4fb..1132af92 100644 --- a/src/impl/diagnose_SD_attributes/particles_impl_fill_outbuf.ipp +++ b/src/impl/diagnose_SD_attributes/particles_impl_fill_outbuf.ipp @@ -70,7 +70,7 @@ namespace libcloudphxx // however, if output was done concurrently, values in the diagnosed vector might change after the call to fill_attr_outbuf. std::vector out(n_part); thrust::copy( - dv.begin(), dv.end(), + dv.begin(), dv.begin() + n_part, out.begin() ); return out;