From 01ce76ffefac73414b9a5d27312d086386569c01 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Tue, 17 Feb 2026 14:20:46 +0100 Subject: [PATCH 01/13] diag_sstp_cond_mom + fix of activation adaptation (but it is to be removed anyway?) --- bindings/python/lib.cpp | 1 + include/libcloudph++/lgrngn/particles.hpp | 3 + ...lculate_noncond_perparticle_sstp_delta.ipp | 4 +- ...erparticle_nomixing_adaptive_sstp_cond.ipp | 55 ++++++++++--------- .../release_arrays_for_perparticle_sstp.ipp | 2 +- .../particles_impl_moms.ipp | 52 ++++++++++++++++++ src/impl/particles_impl.ipp | 25 ++++++--- src/particles_diag.ipp | 11 ++++ src/particles_step.ipp | 2 +- 9 files changed, 119 insertions(+), 36 deletions(-) diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index 6719841fd..eccd1c885 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -427,6 +427,7 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def("diag_ice_a_mom", &lgr::particles_proto_t::diag_ice_a_mom) .def("diag_ice_c_mom", &lgr::particles_proto_t::diag_ice_c_mom) .def("diag_ice_mix_ratio", &lgr::particles_proto_t::diag_ice_mix_ratio) + .def("diag_sstp_cond_mom", &lgr::particles_proto_t::diag_sstp_cond_mom) .def("outbuf", &lgrngn::outbuf) .def("get_attr", &lgr::particles_proto_t::get_attr) ; diff --git a/include/libcloudph++/lgrngn/particles.hpp b/include/libcloudph++/lgrngn/particles.hpp index c0546e54d..6157e067c 100644 --- a/include/libcloudph++/lgrngn/particles.hpp +++ b/include/libcloudph++/lgrngn/particles.hpp @@ -118,6 +118,7 @@ namespace libcloudphxx virtual void diag_vp_mom(const int&) { assert(false); } virtual void diag_wp_mom(const int&) { assert(false); } virtual void diag_incloud_time_mom(const int&) { assert(false); } // requires opts_init.diag_incloud_time==true + virtual void diag_sstp_cond_mom(const int &k) { assert(false); } virtual void diag_max_rw() { assert(false); } virtual void diag_vel_div() { assert(false); } virtual std::map diag_puddle() { assert(false); return std::map(); } @@ -212,6 +213,7 @@ namespace libcloudphxx void diag_vp_mom(const int&); void diag_wp_mom(const int&); void diag_incloud_time_mom(const int &k); + void diag_sstp_cond_mom(const int &k); void diag_wet_mass_dens(const real_t&, const real_t&); void diag_chem(const enum common::chem::chem_species_t&); @@ -321,6 +323,7 @@ namespace libcloudphxx void diag_vp_mom(const int&); void diag_wp_mom(const int&); void diag_incloud_time_mom(const int&); + void diag_sstp_cond_mom(const int &k); void diag_wet_mass_dens(const real_t&, const real_t&); std::vector get_attr(const std::string &); real_t *outbuf(); diff --git a/src/impl/condensation/perparticle/calculate_noncond_perparticle_sstp_delta.ipp b/src/impl/condensation/perparticle/calculate_noncond_perparticle_sstp_delta.ipp index 9cf10e1c1..864862f08 100644 --- a/src/impl/condensation/perparticle/calculate_noncond_perparticle_sstp_delta.ipp +++ b/src/impl/condensation/perparticle/calculate_noncond_perparticle_sstp_delta.ipp @@ -22,7 +22,7 @@ namespace libcloudphxx for (int ix = 0; ix < (opts_init.const_p ? n : n-1); ++ix) { - const real_t sstp = sstp_cond; +// const real_t sstp = sstp_cond; thrust::transform( thrust::make_permutation_iterator(scl[ix]->begin(), ijk.begin()), thrust::make_permutation_iterator(scl[ix]->begin(), ijk.end()), @@ -33,4 +33,4 @@ namespace libcloudphxx } } }; -}; \ No newline at end of file +}; 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 c1c03db58..fd64b238b 100644 --- a/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp +++ b/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp @@ -129,34 +129,34 @@ namespace libcloudphxx // check drw convergence for increasing number of substeps for(int sstp_cond_try = 1; sstp_cond_try <= sstp_cond_max; sstp_cond_try*=2) + //for(int sstp_cond_try = 1; sstp_cond_try <= 1; sstp_cond_try*=2) { delta_fraction_applied = sstp_cond_try == 1 ? 1 : -real_t(1) / sstp_cond_try; _apply_noncond_perparticle_sstp_delta(delta_fraction_applied); // _cond_perparticle_drw2(sstp_cond_try, sstp_cond_try == 1 ? drw2 : drw2_new); // also updates Tp! - - _calc_Tp(); - _calc_sstp_tmp_p(); - _calc_RH(); - (sstp_cond_try == 1 ? drw2 : drw2_new) = - detail::advance_rw2(dt / sstp_cond_try, RH_max, eps_tolerance, cond_mlt, n_iter)( - rw2, - thrust::make_tuple( + _calc_Tp(); + _calc_sstp_tmp_p(); + _calc_RH(); + (sstp_cond_try == 1 ? drw2 : drw2_new) = + detail::advance_rw2(dt / sstp_cond_try, RH_max, eps_tolerance, cond_mlt, n_iter)( + rw2, thrust::make_tuple( - sstp_tmp_rh, - sstp_tmp_rv, - Tp, - detail::common__vterm__visc()(Tp), - rd3, - kpa, - vt, - lambda_D, - lambda_K - ), - sstp_tmp_p, - RH - ) - ); + thrust::make_tuple( + sstp_tmp_rh, + sstp_tmp_rv, + Tp, + detail::common__vterm__visc()(Tp), + rd3, + kpa, + vt, + lambda_D, + lambda_K + ), + sstp_tmp_p, + RH + ) + ); if(sstp_cond_try > 1) // check for convergence { @@ -178,8 +178,13 @@ namespace libcloudphxx { const real_t rc2 = thrust::get<2>(thrust::get<2>(tpl)); - if ( ( rw2 < rc2 && (rw2 + sstp_cond * drw2) > rc2 ) || - ( rw2 > rc2 && (rw2 + sstp_cond * drw2) < rc2 ) ) + // if(rw2 > 10. * rc2 && rw2 < 100. * rc2) + // printf("rw2=%e, rc2=%e, drw2=%e\n", rw2, rc2, drw2); + + // if ( ( rw2 < 10. * rc2 && (rw2 + sstp_cond * drw2) * 10 > rc2 ) || + // ( rw2 > rc2 && (rw2 + sstp_cond * drw2) < rc2 ) ) + + if ( rw2 < real_t(144) * rc2 && rw2 > rc2 / real_t(4)) { sstp_cond = sstp_cond_act; first_cond_step_done_in_adaptation = false; @@ -187,7 +192,7 @@ namespace libcloudphxx } if(!first_cond_step_done_in_adaptation) { - _apply_noncond_perparticle_sstp_delta(delta_fraction_applied); // revert to state before adaptation loop (beacause sstp_cond == sstp_cond_max and sstp_cond_max may not be a power of 2) + _apply_noncond_perparticle_sstp_delta(sstp_cond_max == 1 ? -delta_fraction_applied : delta_fraction_applied); // revert to state before adaptation loop (beacause sstp_cond == sstp_cond_max and sstp_cond_max may not be a power of 2); If only one step was tried, whole change was applied; If more steps were tried, we are moving back from the entire step } } diff --git a/src/impl/condensation/perparticle/release_arrays_for_perparticle_sstp.ipp b/src/impl/condensation/perparticle/release_arrays_for_perparticle_sstp.ipp index ea8d4de6c..7e39adc36 100644 --- a/src/impl/condensation/perparticle/release_arrays_for_perparticle_sstp.ipp +++ b/src/impl/condensation/perparticle/release_arrays_for_perparticle_sstp.ipp @@ -21,7 +21,7 @@ namespace libcloudphxx if(opts_init.adaptive_sstp_cond) { - perparticle_sstp_cond_gp.reset(); + // perparticle_sstp_cond_gp.reset(); } } }; diff --git a/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp b/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp index 1895149e0..e29ca23ff 100644 --- a/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp +++ b/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp @@ -396,5 +396,57 @@ namespace libcloudphxx { moms_calc(vec_bgn, n_part, power, specific); } + + // like moms_calc, but not accounting for multiplicity (n) and not dividing by volume (or air mass) + template + template // iterator type + void particles_t::impl::SD_moms_calc( + const it_t &vec_bgn, + const thrust_size_t npart, + const real_t power + ) + { + thrust::pair< + thrust_device::vector::iterator, + typename thrust_device::vector::iterator + > it_pair = thrust::reduce_by_key( + // input - keys + sorted_ijk.begin(), sorted_ijk.begin()+npart, + // input - values + thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple( + thrust::make_constant_iterator(1), + thrust::make_permutation_iterator(vec_bgn, sorted_id.begin()) + )), + detail::moment_counter(power) + ), + // output - keys + count_ijk.begin(), + // output - values + count_mom.begin() + ); + + count_n = it_pair.first - count_ijk.begin(); +#if !defined(NDEBUG) + { + int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n, isnaninf(), 0, thrust::plus()); + if(nan_count>0) + { + std::cout << nan_count << " nan/inf numbers detected in count_mom after reduce_by_key " << std::endl; + } + } +#endif + assert(count_n <= n_cell); + } + + template + template // iterator type + void particles_t::impl::SD_moms_calc( + const it_t &vec_bgn, + const real_t power + ) + { + moms_calc(vec_bgn, n_part, power); + } }; }; diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 681bf6891..b0a39f1d3 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -494,13 +494,13 @@ namespace libcloudphxx distmem_n_vctrs.insert(&n); // number of required temporary real vectors of size npart - int tmp_drp_no = 1; - if(n_dims == 2) - tmp_drp_no = std::max(tmp_drp_no, 2); - if(allow_sstp_cond) - tmp_drp_no = std::max(tmp_drp_no, 2); - if(opts_init.ice_switch && opts_init.time_dep_ice_nucl) - tmp_drp_no = std::max(tmp_drp_no, 2); + int tmp_drp_no = 2; + // if(n_dims == 2) + // tmp_drp_no = std::max(tmp_drp_no, 2); + // if(allow_sstp_cond) + // tmp_drp_no = std::max(tmp_drp_no, 2); + // if(opts_init.ice_switch && opts_init.time_dep_ice_nucl) + // tmp_drp_no = std::max(tmp_drp_no, 2); if(opts_init.chem_switch) tmp_drp_no = std::max(tmp_drp_no, 3); if(n_dims == 3) @@ -670,6 +670,17 @@ namespace libcloudphxx const real_t power, const bool specific = true ); + template // iterator type + void SD_moms_calc( + const it_t &vec_bgn, + const thrust_size_t npart, + const real_t power + ); + template // iterator type + void SD_moms_calc( + const it_t &vec_bgn, + const real_t power + ); void mass_dens_estim( const typename thrust_device::vector::iterator &vec_bgn, diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index 720a39355..0e5459824 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -654,5 +654,16 @@ namespace libcloudphxx { return pimpl->output_puddle; } + + template + void particles_t::diag_sstp_cond_mom(const int &n) + { + if(pimpl->opts_init.exact_sstp_cond && (pimpl->sstp_cond > 1 || pimpl->sstp_cond_act > 1) && pimpl->opts_init.adaptive_sstp_cond) + { + pimpl->SD_moms_calc(pimpl->perparticle_sstp_cond_gp->get().begin(), n); + } + else + assert(0 && "diag_sstp_cond_mom called, but adaptive substepping is off (opts_init.exact_sstp_cond && (opts_ini.sstp_cond > 1 || opts_ini.sstp_cond_act > 1) && opts_init.adaptive_sstp_cond) == False. Therefore number of substeps is defined by opts_init.sstp_cond."); + } }; }; diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 5e5a9f233..11ee839a4 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -205,7 +205,7 @@ namespace libcloudphxx } pimpl->acquire_arrays_for_perparticle_sstp(); // sstp_dlt_rv_gp, etc. ; as in sstp_percell_step_exact() - pimpl->calculate_noncond_perparticle_sstp_delta(); // sstp_dlt_rv_gp, etc. ; as in sstp_percell_step_exact(); returns change / sstp_count; make it just change and multiply afterwards? + pimpl->calculate_noncond_perparticle_sstp_delta(); // sstp_dlt_rv_gp, etc. ; as in sstp_percell_step_exact(); // adaptive per-particle substepping if(pimpl->opts_init.adaptive_sstp_cond) From 6fbd5891dc28305cac67f8b05bdb25f216e1db74 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Fri, 20 Feb 2026 14:19:19 +0100 Subject: [PATCH 02/13] remove merge leftobers --- ...erparticle_nomixing_adaptive_sstp_cond.ipp | 336 ------------------ 1 file changed, 336 deletions(-) delete mode 100644 src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp diff --git a/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp b/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp deleted file mode 100644 index 7fcdfaf3b..000000000 --- a/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp +++ /dev/null @@ -1,336 +0,0 @@ -namespace libcloudphxx -{ - namespace lgrngn - { - namespace detail - { - template - struct perparticle_nomixing_adaptive_sstp_cond_loop - { - const bool th_dry, const_p, turb_cond, adaptive_sstp_cond; - const real_t dt, RH_max, cond_mlt, sstp_cond_adapt_drw2_eps, sstp_cond_adapt_drw2_max; - const common::detail::eps_tolerance eps_tolerance; - const int n_dims, sstp_cond_max, sstp_cond_act; - const RH_formula_t RH_formula; - uintmax_t n_iter; - - perparticle_nomixing_adaptive_sstp_cond_loop( - const opts_init_t &opts_init, - const opts_t &opts, - const int n_dims, - const real_t &dt, - const int &sstp_cond_max, - const int &sstp_cond_act, - const common::detail::eps_tolerance &eps_tolerance, - const real_t &cond_mlt, - const uintmax_t &n_iter - ) : th_dry(opts_init.th_dry), - const_p(opts_init.const_p), - turb_cond(opts.turb_cond), - dt(dt), - RH_max(opts.RH_max), - n_dims(n_dims), - RH_formula(opts_init.RH_formula), - eps_tolerance(eps_tolerance), - cond_mlt(cond_mlt), - n_iter(n_iter), - adaptive_sstp_cond(opts_init.adaptive_sstp_cond), - sstp_cond_act(sstp_cond_act), - sstp_cond_max(sstp_cond_max), - sstp_cond_adapt_drw2_eps(opts_init.sstp_cond_adapt_drw2_eps), - sstp_cond_adapt_drw2_max(opts_init.sstp_cond_adapt_drw2_max) - {} - - template - BOOST_GPU_ENABLED void operator()( - tpl_t tpl - ) //noexcept - { - // copy values into local variables - // variables that are not modified - const real_t sstp_dlt_rv = thrust::get<5>(thrust::get<0>(tpl)); - const real_t sstp_dlt_th = thrust::get<6>(thrust::get<0>(tpl)); - const real_t sstp_dlt_rhod = thrust::get<7>(thrust::get<0>(tpl)); - const real_t sstp_dlt_p = thrust::get<8>(thrust::get<0>(tpl)); - const auto n = thrust::get<2>(thrust::get<1>(tpl)); - const real_t dv = thrust::get<3>(thrust::get<1>(tpl)); - 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 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; - - // variables that are modified, we make local copies regardless and copy back at the end - unsigned int sstp_cond; // its set in this function, old value not important - real_t sstp_tmp_rv = thrust::get<1>(thrust::get<0>(tpl)); - real_t sstp_tmp_th = thrust::get<2>(thrust::get<0>(tpl)); - real_t sstp_tmp_rh = thrust::get<3>(thrust::get<0>(tpl)); - real_t sstp_tmp_p = const_p ? thrust::get<4>(thrust::get<0>(tpl)) : 0; - real_t ssp = turb_cond ? thrust::get<9>(thrust::get<0>(tpl)) : 0; - real_t rw2 = thrust::get<1>(thrust::get<1>(tpl)); - - real_t drw2, Tp, RH; - - // helper functions - auto _apply_noncond_perparticle_sstp_delta = [&] (const real_t &multiplier) -> void - { - sstp_tmp_rv += sstp_dlt_rv * multiplier; - sstp_tmp_th += sstp_dlt_th * multiplier; - sstp_tmp_rh += sstp_dlt_rhod * multiplier; - if(const_p) - sstp_tmp_p += sstp_dlt_p * multiplier; - if(turb_cond) - ssp += dot_ssp * dt * multiplier; - }; - - auto _calc_Tp = [&] () -> void - { - if(th_dry) - Tp = detail::common__theta_dry__T_rhod()(sstp_tmp_th, sstp_tmp_rh); - else - Tp = detail::common__theta_std__T_p()(sstp_tmp_th, - thrust::make_tuple(sstp_tmp_rv, sstp_tmp_p) - ); - }; - - auto _calc_sstp_tmp_p = [&] () -> void - { - if(!const_p) // sstp_tmp_p needs to be allocated even without const_p! - sstp_tmp_p = detail::common__theta_dry__p()( - thrust::make_tuple(sstp_tmp_rh, sstp_tmp_rv, Tp) - ); - }; - - auto _calc_RH = [&] () -> void - { - RH = turb_cond ? - detail::RH_hlpr(RH_formula)( - thrust::make_tuple(sstp_tmp_p, sstp_tmp_rv, Tp, ssp) - ) : - detail::RH_hlpr(RH_formula)( - thrust::make_tuple(sstp_tmp_p, sstp_tmp_rv, Tp) - ); - }; - - // bool converged = false; - real_t delta_fraction_applied; - bool first_cond_step_done_in_adaptation = sstp_cond_max == 1 ? true : false; // actually its true if sstp_cond_max is a power of 2 (?) - // bool activates = false; - - // look for correct number of substeps - // NOTE: this function is actually only called when adaptive_sstp_cond == true, so we skip the check below - // if(adaptive_sstp_cond) - { - real_t drw2_new; - // real_t Tp; // temperature - - sstp_cond = sstp_cond_max; // start with max number of substeps, may be changed due to convergence or if droplets activate in this step - - // check drw convergence for increasing number of substeps - for(int sstp_cond_try = 1; sstp_cond_try <= sstp_cond_max; sstp_cond_try*=2) - //for(int sstp_cond_try = 1; sstp_cond_try <= 1; sstp_cond_try*=2) - { - delta_fraction_applied = sstp_cond_try == 1 ? 1 : -real_t(1) / sstp_cond_try; - _apply_noncond_perparticle_sstp_delta(delta_fraction_applied); - // _cond_perparticle_drw2(sstp_cond_try, sstp_cond_try == 1 ? drw2 : drw2_new); // also updates Tp! - - _calc_Tp(); - _calc_sstp_tmp_p(); - _calc_RH(); - (sstp_cond_try == 1 ? drw2 : drw2_new) = - detail::advance_rw2(dt / sstp_cond_try, RH_max, eps_tolerance, cond_mlt, n_iter)( - rw2, - thrust::make_tuple( - thrust::make_tuple( - sstp_tmp_rh, - sstp_tmp_rv, - Tp, - detail::common__vterm__visc()(Tp), - rd3, - kpa, - vt, - lambda_D, - lambda_K - ), - sstp_tmp_p, - RH - ) - ); - - if(sstp_cond_try > 1) // check for convergence - { - if((cuda::std::abs(drw2_new * 2 - drw2) <= sstp_cond_adapt_drw2_eps * rw2) // drw2 relative to rw2 converged - && cuda::std::abs(drw2 < sstp_cond_adapt_drw2_max * rw2)) // otherwise for small droplets (near activation?) drw2_new == 2*drw already for 2 substeps, but we ativate too many droplets - // if(cuda::std::abs(drw2_new * 2 - drw2) <= tol * drw2) // drw2 converged - { - sstp_cond = sstp_cond_try / 2; - _apply_noncond_perparticle_sstp_delta(-delta_fraction_applied); // revert last addition to get to a state after one step of converged number - first_cond_step_done_in_adaptation = true; - break; - } - drw2 = drw2_new; - } - } - - // override number of substeps for SDs that de/activate in this timestep; - if(sstp_cond_act > 1) - { - const real_t rc2 = thrust::get<2>(thrust::get<2>(tpl)); - - // if(rw2 > 10. * rc2 && rw2 < 100. * rc2) - // printf("rw2=%e, rc2=%e, drw2=%e\n", rw2, rc2, drw2); - - // if ( ( rw2 < 10. * rc2 && (rw2 + sstp_cond * drw2) * 10 > rc2 ) || - // ( rw2 > rc2 && (rw2 + sstp_cond * drw2) < rc2 ) ) - - if ( rw2 < real_t(144) * rc2 && rw2 > rc2 / real_t(4)) - { - sstp_cond = sstp_cond_act; - first_cond_step_done_in_adaptation = false; - } - } - if(!first_cond_step_done_in_adaptation) - { - _apply_noncond_perparticle_sstp_delta(sstp_cond_max == 1 ? -delta_fraction_applied : delta_fraction_applied); // revert to state before adaptation loop (beacause sstp_cond == sstp_cond_max and sstp_cond_max may not be a power of 2); If only one step was tried, whole change was applied; If more steps were tried, we are moving back from the entire step - } - } - - delta_fraction_applied = real_t(1) / sstp_cond; - auto _advance_rw2 = detail::advance_rw2(dt / sstp_cond, RH_max, eps_tolerance, cond_mlt, n_iter); - real_t &rw3 = drw2; // drw2 needed only at the start of the first step - real_t drw3; - - auto rw2torw3 = detail::rw2torwX(); - - // actual condensation substepping - for(int step = 0; step < sstp_cond; ++step) - { - drw3 = step > 0 ? -rw3 : -rw2torw3(rw2); - - if(first_cond_step_done_in_adaptation && step == 0) - { - rw2 += drw2; - } - else - { - _apply_noncond_perparticle_sstp_delta(delta_fraction_applied); - _calc_Tp(); - _calc_sstp_tmp_p(); - _calc_RH(); - - rw2 = _advance_rw2( - rw2, - thrust::make_tuple( - thrust::make_tuple( - sstp_tmp_rh, - sstp_tmp_rv, - Tp, - detail::common__vterm__visc()(Tp), - rd3, - kpa, - vt, - lambda_D, - lambda_K - ), - sstp_tmp_p, - RH - ) - ); - } - - if (step < sstp_cond - 1) - { - rw3 = rw2torw3(rw2); - drw3 += rw3; - } - else - drw3 += rw2torw3(rw2); - - drw3 = detail::rw3diff2drv( - - common::moist_air::rho_w() / si::kilograms * si::cubic_metres - * real_t(4./3) * real_t(3.14159265358979323846264338), n_dims - ) (drw3, - thrust::make_tuple(sstp_tmp_rh, n, dv) - ); - - sstp_tmp_rv += drw3; - - drw3 = detail::dth()( - thrust::make_tuple(drw3, Tp, sstp_tmp_th) - ); - - sstp_tmp_th += drw3; - } - - // copy back modified variables - thrust::get<0>(thrust::get<0>(tpl)) = sstp_cond; - thrust::get<1>(thrust::get<0>(tpl)) = sstp_tmp_rv; - thrust::get<2>(thrust::get<0>(tpl)) = sstp_tmp_th; - thrust::get<3>(thrust::get<0>(tpl)) = sstp_tmp_rh; - if(const_p) - thrust::get<4>(thrust::get<0>(tpl)) = sstp_tmp_p; - if(turb_cond) - thrust::get<9>(thrust::get<0>(tpl)) = ssp; - thrust::get<1>(thrust::get<1>(tpl)) = rw2; - } - }; - }; - - template - void particles_t::impl::perparticle_nomixing_adaptive_sstp_cond(const opts_t &opts) { - - auto &perparticle_sstp_cond = perparticle_sstp_cond_gp->get(); - auto &sstp_dlt_rv = sstp_dlt_rv_gp->get(); - auto &sstp_dlt_th = sstp_dlt_th_gp->get(); - auto &sstp_dlt_rhod = sstp_dlt_rhod_gp->get(); - auto &sstp_dlt_p = sstp_dlt_p_gp->get(); - // auto &Tp = Tp_gp->get(); - // auto &drwX = drwX_gp->get(); - // auto &rwX = rwX_gp->get(); - const auto &lambda_D = lambda_D_gp->get(); - const auto &lambda_K = lambda_K_gp->get(); - - auto pptcl_nomix_sstp_cond_args_zip = - thrust::make_zip_iterator(thrust::make_tuple( - thrust::make_zip_iterator(thrust::make_tuple( - perparticle_sstp_cond.begin(), - sstp_tmp_rv.begin(), - sstp_tmp_th.begin(), - sstp_tmp_rh.begin(), - sstp_tmp_p.begin(), - sstp_dlt_rv.begin(), - sstp_dlt_th.begin(), - sstp_dlt_rhod.begin(), - sstp_dlt_p.begin(), - ssp.begin() - )), - thrust::make_zip_iterator(thrust::make_tuple( - dot_ssp.begin(), - // Tp.begin(), - // drwX.begin(), - // rwX.begin(), - rw2.begin(), - n.begin(), - 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() - )), - thrust::make_zip_iterator(thrust::make_tuple( - kpa.begin(), - vt.begin(), - rc2.begin() - )) - )); - - thrust::for_each( - pptcl_nomix_sstp_cond_args_zip, - pptcl_nomix_sstp_cond_args_zip + n_part, - detail::perparticle_nomixing_adaptive_sstp_cond_loop( - opts_init, opts, n_dims, dt, sstp_cond, sstp_cond_act, config.eps_tolerance, config.cond_mlt, config.n_iter - ) - ); - } - }; -}; From fd5d91d6b1c141cf8c6f5b2a63952f830181ec42 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 23 Feb 2026 09:59:00 +0100 Subject: [PATCH 03/13] incease tmp_drp_no, because sstp_cond is not released (due to diag_sstp_cond) --- src/impl/particles_impl.ipp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 20ab390c1..224ad1873 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -516,7 +516,7 @@ namespace libcloudphxx if(allow_sstp_cond && opts_init.exact_sstp_cond && sstp_cond_exact_nomix_adaptive) tmp_drp_no = std::max(tmp_drp_no, 4); // why 5? not 4? if(allow_sstp_cond && opts_init.exact_sstp_cond && !sstp_cond_exact_nomix_adaptive) - tmp_drp_no = std::max(tmp_drp_no, 7); // why 8? not 7? + tmp_drp_no = std::max(tmp_drp_no, 8); // 8, because sstp_cond is not released, for diag_sstp_cond to work // if(allow_sstp_cond && opts_init.exact_sstp_cond && opts_init.const_p) // tmp_drp_no = std::max(tmp_drp_no, 7); tmp_device_real_part.add_vectors(tmp_drp_no-1); // -1 because 1 is already created in the ctor From ff57222099bdcc4bffd5c574560739bb8afc11e4 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Thu, 26 Feb 2026 10:28:50 +0100 Subject: [PATCH 04/13] diag ice rng in python --- bindings/python/lib.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index eccd1c885..446a64812 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -426,6 +426,10 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def("diag_water_cons", &lgr::particles_proto_t::diag_water_cons) .def("diag_ice_a_mom", &lgr::particles_proto_t::diag_ice_a_mom) .def("diag_ice_c_mom", &lgr::particles_proto_t::diag_ice_c_mom) + .def("diag_ice_a_rng", &lgr::particles_proto_t::diag_ice_a_rng) + .def("diag_ice_c_rng", &lgr::particles_proto_t::diag_ice_c_rng) + .def("diag_ice_a_rng_cons", &lgr::particles_proto_t::diag_ice_a_rng_cons) + .def("diag_ice_c_rng_cons", &lgr::particles_proto_t::diag_ice_c_rng_cons) .def("diag_ice_mix_ratio", &lgr::particles_proto_t::diag_ice_mix_ratio) .def("diag_sstp_cond_mom", &lgr::particles_proto_t::diag_sstp_cond_mom) .def("outbuf", &lgrngn::outbuf) From 81bfd89a38adce8f2f3957b5658caae1124242f8 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Thu, 26 Feb 2026 12:11:43 +0100 Subject: [PATCH 05/13] bindigns missing functions --- bindings/python/lib.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index 446a64812..b0802b488 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -131,10 +131,13 @@ BOOST_PYTHON_MODULE(libcloudphxx) bp::def("exner", &common::exner); bp::def("p_v", &common::p_v); bp::def("p_vs", &common::p_vs); + bp::def("p_vsi", &common::p_vsi); bp::def("r_vs", &common::r_vs); bp::def("r_vsi", &common::r_vsi); bp::def("p_vs_tet", &common::p_vs_tet); bp::def("l_v", &common::l_v); + bp::def("l_s", &common::l_s); + bp::def("l_f", &common::l_f); bp::def("T", &common::T); bp::def("p", &common::p); bp::def("visc", &common::visc); From 8e8b6705a05ff384a41694e1be4a27470935a4ac Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Thu, 26 Feb 2026 13:58:19 +0100 Subject: [PATCH 06/13] debug --- bindings/python/lib.cpp | 1 + src/impl/particles_impl.ipp | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index b0802b488..677d1677b 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -284,6 +284,7 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("turb_cond", &lgr::opts_t::turb_cond) .def_readwrite("turb_coal", &lgr::opts_t::turb_coal) .def_readwrite("ice_nucl", &lgr::opts_t::ice_nucl) + .def_readwrite("depo", &lgr::opts_t::depo) .def_readwrite("dt", &lgr::opts_t::dt) .add_property("src_dry_distros", &lgrngn::get_sdd, &lgrngn::set_sdd) .add_property("src_dry_sizes", &lgrngn::get_ds, &lgrngn::set_sds) diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 368952876..ea50e97ae 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -516,7 +516,7 @@ namespace libcloudphxx if(allow_sstp_cond && opts_init.exact_sstp_cond && sstp_cond_exact_nomix_adaptive) tmp_drp_no = std::max(tmp_drp_no, 4); // why 5? not 4? if(allow_sstp_cond && opts_init.exact_sstp_cond && !sstp_cond_exact_nomix_adaptive) - tmp_drp_no = std::max(tmp_drp_no, 8); // 8, because sstp_cond is not released, for diag_sstp_cond to work + tmp_drp_no = std::max(tmp_drp_no, 10); // for some reason it fails for less than 10 // if(allow_sstp_cond && opts_init.exact_sstp_cond && opts_init.const_p) // tmp_drp_no = std::max(tmp_drp_no, 7); tmp_device_real_part.add_vectors(tmp_drp_no-1); // -1 because 1 is already created in the ctor From ddfd80be446581de347dc3eeb9fcb433488896a2 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Tue, 3 Mar 2026 20:14:54 +0100 Subject: [PATCH 07/13] debugging segfault --- .../perparticle/acquire_arrays_for_perparticle_sstp.ipp | 1 + .../perparticle/release_arrays_for_perparticle_sstp.ipp | 1 + src/impl/particles_impl.ipp | 2 +- 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp b/src/impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp index 7e612960a..422bee3c2 100644 --- a/src/impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp +++ b/src/impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp @@ -24,6 +24,7 @@ namespace libcloudphxx { reset_guardp(ice_mass_gp, tmp_device_real_part); reset_guardp(d_ice_mass_gp, tmp_device_real_part); + reset_guardp(d_ice_mass_percell_gp, tmp_device_real_part); } } diff --git a/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp b/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp index 097ddb50e..c073fc947 100644 --- a/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp +++ b/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp @@ -25,6 +25,7 @@ namespace libcloudphxx { ice_mass_gp.reset(); d_ice_mass_gp.reset(); + d_ice_mass_percell_gp.reset(); } } diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index ea50e97ae..092a1409b 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -516,7 +516,7 @@ namespace libcloudphxx if(allow_sstp_cond && opts_init.exact_sstp_cond && sstp_cond_exact_nomix_adaptive) tmp_drp_no = std::max(tmp_drp_no, 4); // why 5? not 4? if(allow_sstp_cond && opts_init.exact_sstp_cond && !sstp_cond_exact_nomix_adaptive) - tmp_drp_no = std::max(tmp_drp_no, 10); // for some reason it fails for less than 10 + tmp_drp_no = std::max(tmp_drp_no, 11); // for some reason it fails for less than 11 // if(allow_sstp_cond && opts_init.exact_sstp_cond && opts_init.const_p) // tmp_drp_no = std::max(tmp_drp_no, 7); tmp_device_real_part.add_vectors(tmp_drp_no-1); // -1 because 1 is already created in the ctor From 7b0bc9264259d50044146f90896722684665076e Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Wed, 4 Mar 2026 11:41:30 +0100 Subject: [PATCH 08/13] d_ice_mass per particle --- .../perparticle/acquire_arrays_for_perparticle_sstp.ipp | 1 - ...particle_drw3_or_d_ice_mass_to_perparticle_rv_and_th.ipp | 6 +++--- .../perparticle/release_arrays_for_perparticle_sstp.ipp | 1 - src/impl/particles_impl.ipp | 2 +- 4 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp b/src/impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp index 422bee3c2..7e612960a 100644 --- a/src/impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp +++ b/src/impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp @@ -24,7 +24,6 @@ namespace libcloudphxx { reset_guardp(ice_mass_gp, tmp_device_real_part); reset_guardp(d_ice_mass_gp, tmp_device_real_part); - reset_guardp(d_ice_mass_percell_gp, tmp_device_real_part); } } diff --git a/src/impl/condensation_deposition/perparticle/apply_perparticle_drw3_or_d_ice_mass_to_perparticle_rv_and_th.ipp b/src/impl/condensation_deposition/perparticle/apply_perparticle_drw3_or_d_ice_mass_to_perparticle_rv_and_th.ipp index 066d1e757..e58378ac4 100644 --- a/src/impl/condensation_deposition/perparticle/apply_perparticle_drw3_or_d_ice_mass_to_perparticle_rv_and_th.ipp +++ b/src/impl/condensation_deposition/perparticle/apply_perparticle_drw3_or_d_ice_mass_to_perparticle_rv_and_th.ipp @@ -65,7 +65,7 @@ namespace libcloudphxx // add the change in ice mass // TODO: very similar to the part for liquid water... if(ice_mass_changed) { - thrust_device::vector &d_ice_mass = d_ice_mass_percell_gp->get(); + thrust_device::vector &d_ice_mass = d_ice_mass_gp->get(); // calculate rv change from the change in solid water mass thrust::transform( d_ice_mass.begin(), d_ice_mass.end(), @@ -106,11 +106,11 @@ namespace libcloudphxx ); } - thrust_device::vector &dth = rw3_changed ? drwX_gp->get() : d_ice_mass_percell_gp->get(); + thrust_device::vector &dth = rw3_changed ? drwX_gp->get() : d_ice_mass_gp->get(); if(rw3_changed && ice_mass_changed) thrust::transform( dth.begin(), dth.end(), - d_ice_mass_percell_gp->get().begin(), + d_ice_mass_gp->get().begin(), dth.begin(), thrust::plus() ); diff --git a/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp b/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp index c073fc947..097ddb50e 100644 --- a/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp +++ b/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp @@ -25,7 +25,6 @@ namespace libcloudphxx { ice_mass_gp.reset(); d_ice_mass_gp.reset(); - d_ice_mass_percell_gp.reset(); } } diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 092a1409b..ea50e97ae 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -516,7 +516,7 @@ namespace libcloudphxx if(allow_sstp_cond && opts_init.exact_sstp_cond && sstp_cond_exact_nomix_adaptive) tmp_drp_no = std::max(tmp_drp_no, 4); // why 5? not 4? if(allow_sstp_cond && opts_init.exact_sstp_cond && !sstp_cond_exact_nomix_adaptive) - tmp_drp_no = std::max(tmp_drp_no, 11); // for some reason it fails for less than 11 + tmp_drp_no = std::max(tmp_drp_no, 10); // for some reason it fails for less than 10 // if(allow_sstp_cond && opts_init.exact_sstp_cond && opts_init.const_p) // tmp_drp_no = std::max(tmp_drp_no, 7); tmp_device_real_part.add_vectors(tmp_drp_no-1); // -1 because 1 is already created in the ctor From f6aa6891accab1fe41b545474feea11540a826e1 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Thu, 19 Mar 2026 12:11:59 +0100 Subject: [PATCH 09/13] diag_sstp_mom only if sstp vector exists --- src/particles_diag.ipp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index 0e5459824..1c7c1e36d 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -660,7 +660,10 @@ namespace libcloudphxx { if(pimpl->opts_init.exact_sstp_cond && (pimpl->sstp_cond > 1 || pimpl->sstp_cond_act > 1) && pimpl->opts_init.adaptive_sstp_cond) { - pimpl->SD_moms_calc(pimpl->perparticle_sstp_cond_gp->get().begin(), n); + if (pimpl->perparticle_sstp_cond_gp) + { + pimpl->SD_moms_calc(pimpl->perparticle_sstp_cond_gp->get().begin(), n); + } } else assert(0 && "diag_sstp_cond_mom called, but adaptive substepping is off (opts_init.exact_sstp_cond && (opts_ini.sstp_cond > 1 || opts_ini.sstp_cond_act > 1) && opts_init.adaptive_sstp_cond) == False. Therefore number of substeps is defined by opts_init.sstp_cond."); From 2bc5462210d58303675e98fc0d0de52986db8405 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Thu, 9 Apr 2026 20:42:41 +0200 Subject: [PATCH 10/13] diag_sstp_cond_mom on GPU --- src/particles_multi_gpu_diag.ipp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/particles_multi_gpu_diag.ipp b/src/particles_multi_gpu_diag.ipp index 797b134f3..2f12a33fc 100644 --- a/src/particles_multi_gpu_diag.ipp +++ b/src/particles_multi_gpu_diag.ipp @@ -213,6 +213,12 @@ namespace libcloudphxx pimpl->mcuda_run(&particles_t::diag_incloud_time_mom, k); } + template + void particles_t::diag_sstp_cond_mom(const int &k) + { + pimpl->mcuda_run(&particles_t::diag_sstp_cond_mom, k); + } + template void particles_t::diag_wet_mass_dens(const real_t &a, const real_t &b) { From 3d5cb7c7413d446a00a77f2fc2121e407d82c222 Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Fri, 17 Apr 2026 15:24:52 +0200 Subject: [PATCH 11/13] Revert "Merge branch 'ice' into diag_sstp_cond" This reverts commit c905d9198d8a9b1716bf6af9e7f3827be54acb7f, reversing changes made to 2bc5462210d58303675e98fc0d0de52986db8405. --- .../libcloudph++/common/ice_deposition.hpp | 136 ------------ 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 +-- .../ice/particles_impl_ice_dep_common.ipp | 200 ------------------ .../particles_impl_init_T_freeze.ipp | 4 +- 6 files changed, 11 insertions(+), 390 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 549d681db..000000000 --- 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/models/kinematic_2D/cases/icmw8_case1.hpp b/models/kinematic_2D/cases/icmw8_case1.hpp index 76133a3fb..3f44a4f79 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 b228fd48b..29c964985 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 d6635168b..924921d24 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(); 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 d51bb3227..000000000 --- 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); -// } -// }; -// -// }; -// }; -// }; diff --git a/src/impl/initialization/particles_impl_init_T_freeze.ipp b/src/impl/initialization/particles_impl_init_T_freeze.ipp index ce2097db6..ee0f5da55 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(rd2_insol.begin() + n_part_old, u01.begin() + n_part_old)), + thrust::make_zip_iterator(thrust::make_tuple(rd2_insol.end(), u01.end())), T_freeze.begin() + n_part_old, T_freeze_CDF_inv_functor(opts_init.inp_type) ); From b2cb994d0b18214635c107e31c7980d8bd8b270a Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Fri, 17 Apr 2026 15:04:53 +0200 Subject: [PATCH 12/13] 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 ee0f5da55..ce2097db6 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 7ddcf00f0d2e3438d3bf30ce3065c70e3bfb8a6b Mon Sep 17 00:00:00 2001 From: AgnieszkaMakulska Date: Wed, 22 Apr 2026 15:22:45 +0200 Subject: [PATCH 13/13] parcel moms fix --- src/impl/diagnose_SD_attributes/particles_impl_moms.ipp | 2 +- src/particles_step.ipp | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp b/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp index e29ca23ff..b811501a2 100644 --- a/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp +++ b/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp @@ -316,7 +316,7 @@ namespace libcloudphxx } #endif assert(count_n <= n_cell); - if(specific) + if(specific && n_dims>0) // parcel (n_dims==0) implicitly assumes 1kg of air { // dividing by dv thrust::transform( diff --git a/src/particles_step.ipp b/src/particles_step.ipp index ecc78b669..67206f9e5 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -249,6 +249,7 @@ namespace libcloudphxx } pimpl->release_arrays_for_perparticle_sstp(opts.cond, opts.depo); pimpl->apply_perparticle_cond_change_to_percell_rv_and_th(); + pimpl->n_filtered_gp.reset(); // used in apply_perparticle_cond_change_to_percell_rv_and_th, not needed anymore } else // apply per-cell sstp logic, always with mixing {