diff --git a/include/libcloudph++/lgrngn/opts.hpp b/include/libcloudph++/lgrngn/opts.hpp index 8bf6b59f..fed44209 100644 --- a/include/libcloudph++/lgrngn/opts.hpp +++ b/include/libcloudph++/lgrngn/opts.hpp @@ -20,7 +20,7 @@ namespace libcloudphxx struct opts_t { // process toggling - bool adve, sedi, subs, cond, coal, src, rlx, rcyc, turb_adve, turb_cond, turb_coal, ice_nucl; + bool adve, sedi, subs, cond, coal, src, rlx, rcyc, turb_adve, turb_cond, turb_coal, ice_nucl, depo; // RH limit for drop growth real_t RH_max; @@ -42,7 +42,7 @@ namespace libcloudphxx opts_t() : adve(true), sedi(true), subs(false), cond(true), coal(true), src(false), rlx(false), rcyc(false), chem_dsl(false), chem_dsc(false), chem_rct(false), - turb_adve(false), turb_cond(false), turb_coal(false), ice_nucl(false), + turb_adve(false), turb_cond(false), turb_coal(false), ice_nucl(false), depo(false), RH_max(44), // :) (anything greater than 1.1 would be enough dt(-1) // negative means that we do not override dt in this step { diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index 7acaa132..35d08bc1 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -81,7 +81,7 @@ namespace libcloudphxx // super-droplet advection scheme as_t adve_scheme; - // RH formula + // RH formula, for ice, pv_tet and rv_tet work as pv_cc and rv_cc, respectively, because there is no Tetens formula for ice RH_formula_t RH_formula; // diff --git a/src/impl/common/calc_liq_ice_content_change.ipp b/src/impl/common/calc_liq_ice_content_change.ipp index 0fe73c37..6416bb40 100644 --- a/src/impl/common/calc_liq_ice_content_change.ipp +++ b/src/impl/common/calc_liq_ice_content_change.ipp @@ -28,7 +28,7 @@ namespace libcloudphxx if(opts_init.ice_switch) { - thrust_device::vector &d_ice_mass = d_ice_mass_gp->get(); + thrust_device::vector &d_ice_mass_percell = d_ice_mass_percell_gp->get(); moms_gt0(ice_a.begin()); // choose ice particles (ice_a>0) moms_calc(thrust::make_transform_iterator( @@ -40,8 +40,8 @@ namespace libcloudphxx thrust::transform( count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg - thrust::make_permutation_iterator(d_ice_mass.begin(), count_ijk.begin()), // 2nd arg - thrust::make_permutation_iterator(d_ice_mass.begin(), count_ijk.begin()), // output + thrust::make_permutation_iterator(d_ice_mass_percell.begin(), count_ijk.begin()), // 2nd arg + thrust::make_permutation_iterator(d_ice_mass_percell.begin(), count_ijk.begin()), // output thrust::plus() ); } diff --git a/src/impl/common/particles_impl_update_th_rv.ipp b/src/impl/common/particles_impl_update_th_rv.ipp index 36849b83..c3049f7d 100644 --- a/src/impl/common/particles_impl_update_th_rv.ipp +++ b/src/impl/common/particles_impl_update_th_rv.ipp @@ -103,15 +103,15 @@ namespace libcloudphxx // else if (phase == phase_change::deposition) // if(opts_init.ice_switch) // TODO: call update_th_rv once per cond/depo // { - // thrust_device::vector &d_ice_mass = d_ice_mass_gp->get(); - // nancheck(d_ice_mass, "update_th_rv: input d_ice_mass"); + // thrust_device::vector &d_ice_mass_percell = d_ice_mass_percell_gp->get(); + // nancheck(d_ice_mass_percell, "update_th_rv: input d_ice_mass_percell"); // thrust::transform( - // d_ice_mass.begin(), d_ice_mass.end(), // input - 1st arg + // d_ice_mass_percell.begin(), d_ice_mass_percell.end(), // input - 1st arg // thrust::make_constant_iterator( // input - 2nd arg // - real_t(1) // ), - // d_ice_mass.begin(), // output + // d_ice_mass_percell.begin(), // output // thrust::multiplies() // ); // } @@ -130,12 +130,12 @@ namespace libcloudphxx if(opts_init.ice_switch) { - thrust_device::vector &d_ice_mass = d_ice_mass_gp->get(); - nancheck(d_ice_mass, "update_th_rv: input d_ice_mass"); + thrust_device::vector &d_ice_mass_percell = d_ice_mass_percell_gp->get(); + nancheck(d_ice_mass_percell, "update_th_rv: input d_ice_mass_percell"); thrust::transform( rv.begin(), rv.end(), // input - 1st arg - d_ice_mass.begin(), // input - 2nd arg + d_ice_mass_percell.begin(), // input - 2nd arg rv.begin(), // output thrust::minus() ); @@ -170,13 +170,13 @@ namespace libcloudphxx // else if (phase == phase_change::deposition) if(opts_init.ice_switch) { - thrust_device::vector &d_ice_mass = d_ice_mass_gp->get(); + thrust_device::vector &d_ice_mass_percell = d_ice_mass_percell_gp->get(); thrust::transform( th.begin(), th.end(), // input - 1st arg thrust::make_transform_iterator( zip_it_t(thrust::make_tuple( - d_ice_mass.begin(), // - T.begin(), // dth = d_ice_mass * d_th_d_rv(T, th) + d_ice_mass_percell.begin(), // + T.begin(), // dth = d_ice_mass_percell * d_th_d_rv(T, th) th.begin() // )), detail::dth_dep() @@ -186,7 +186,7 @@ namespace libcloudphxx ); } drw_mom3_gp.reset(); // destroy guard to tmp array that stored change in 3rd moment of rw - d_ice_mass_gp.reset(); // destroy guard to tmp array that stored change in 3rd moment of rw + d_ice_mass_percell_gp.reset(); // destroy guard to tmp array that stored change in 3rd moment of rw nancheck(th, "update_th_rv: th after update"); } diff --git a/src/impl/common/save_liq_ice_content_before_change.ipp b/src/impl/common/save_liq_ice_content_before_change.ipp index af14a02e..3ea76ada 100644 --- a/src/impl/common/save_liq_ice_content_before_change.ipp +++ b/src/impl/common/save_liq_ice_content_before_change.ipp @@ -32,8 +32,8 @@ namespace libcloudphxx if(opts_init.ice_switch) { - reset_guardp(d_ice_mass_gp, tmp_device_real_cell); - thrust_device::vector &d_ice_mass = d_ice_mass_gp->get(); + reset_guardp(d_ice_mass_percell_gp, tmp_device_real_cell); + thrust_device::vector &d_ice_mass_percell = d_ice_mass_percell_gp->get(); moms_gt0(ice_a.begin()); // choose ice particles (ice_a>0) moms_calc(thrust::make_transform_iterator( @@ -46,13 +46,13 @@ namespace libcloudphxx // fill with 0s if not all cells have particles if(count_n!=n_cell) { - thrust::fill(d_ice_mass.begin(), d_ice_mass.end(), real_t(0.)); + thrust::fill(d_ice_mass_percell.begin(), d_ice_mass_percell.end(), real_t(0.)); // thrust::fill(ice_mass.begin(), ice_mass.end(), real_t(0.)); } thrust::transform( count_mom.begin(), count_mom.begin() + count_n, - thrust::make_permutation_iterator(d_ice_mass.begin(), count_ijk.begin()), + thrust::make_permutation_iterator(d_ice_mass_percell.begin(), count_ijk.begin()), thrust::negate() ); } diff --git a/src/impl/condensation/perparticle/apply_perparticle_drw3_to_perparticle_rv_and_th.ipp b/src/impl/condensation/perparticle/apply_perparticle_drw3_to_perparticle_rv_and_th.ipp deleted file mode 100644 index cb1343eb..00000000 --- a/src/impl/condensation/perparticle/apply_perparticle_drw3_to_perparticle_rv_and_th.ipp +++ /dev/null @@ -1,62 +0,0 @@ -namespace libcloudphxx -{ - namespace lgrngn - { - template - void particles_t::impl::apply_perparticle_drw3_to_perparticle_rv_and_th() - { - namespace arg = thrust::placeholders; - thrust_device::vector &drw3 = drwX_gp->get(); - thrust_device::vector &Tp = Tp_gp->get(); - - thrust::transform( - drw3.begin(), drw3.end(), - thrust::make_zip_iterator(thrust::make_tuple( - sstp_tmp_rh.begin(), - n.begin(), - thrust::make_permutation_iterator(dv.begin(), ijk.begin()) - )), - drw3.begin(), - detail::rw3diff2drv( - - common::moist_air::rho_w() / si::kilograms * si::cubic_metres - * real_t(4./3) * pi(), n_dims - ) - ); - - if(opts_init.sstp_cond_mix) - update_pstate(sstp_tmp_rv, drw3); - else - thrust::transform( - drw3.begin(), drw3.end(), - sstp_tmp_rv.begin(), - sstp_tmp_rv.begin(), - thrust::plus() - ); - - thrust::transform( - thrust::make_zip_iterator(thrust::make_tuple( - drw3.begin(), - Tp.begin(), - sstp_tmp_th.begin() - )), - thrust::make_zip_iterator(thrust::make_tuple( - drw3.end(), - Tp.end(), - sstp_tmp_th.end() - )), - drw3.begin(), - detail::dth() - ); - - if(opts_init.sstp_cond_mix) - update_pstate(sstp_tmp_th, drw3); - else - thrust::transform( - drw3.begin(), drw3.end(), - sstp_tmp_th.begin(), - sstp_tmp_th.begin(), - thrust::plus() - ); - } - }; -}; diff --git a/src/impl/condensation/perparticle/cond_perparticle_advance_rw2.ipp b/src/impl/condensation/perparticle/cond_perparticle_advance_rw2.ipp deleted file mode 100644 index cd084e41..00000000 --- a/src/impl/condensation/perparticle/cond_perparticle_advance_rw2.ipp +++ /dev/null @@ -1,128 +0,0 @@ -namespace libcloudphxx -{ - namespace lgrngn - { - namespace detail - { - template - struct RH_sgs - { - RH resolved_RH; - - BOOST_GPU_ENABLED - RH_sgs(RH_formula_t RH_formula): - resolved_RH(RH_formula) - {} - - BOOST_GPU_ENABLED - real_t operator()(const thrust::tuple &tpl) noexcept - { - return resolved_RH(thrust::make_tuple(thrust::get<0>(tpl), thrust::get<1>(tpl), thrust::get<2>(tpl))) + thrust::get<3>(tpl); - } - }; - }; - - - template - void particles_t::impl::cond_perparticle_advance_rw2( - const real_t &RH_max, - const bool turb_cond - ) { - - namespace arg = thrust::placeholders; - - thrust_device::vector &Tp = Tp_gp->get(); - - // calculate perparticle temperature; TODO: skip it and use theta in drw2? - if(opts_init.th_dry) - { - thrust::transform( - sstp_tmp_th.begin(), sstp_tmp_th.end(), - sstp_tmp_rh.begin(), - Tp.begin(), - detail::common__theta_dry__T_rhod() - ); - } - else - { - thrust::transform( - sstp_tmp_th.begin(), sstp_tmp_th.end(), - thrust::make_zip_iterator(thrust::make_tuple( - sstp_tmp_rv.begin(), - sstp_tmp_p.begin() - )), - Tp.begin(), - detail::common__theta_std__T_p() - ); - } - - // advance rw2 - if(!opts_init.const_p) - { - auto pressure_iter = thrust::make_transform_iterator( - thrust::make_zip_iterator(thrust::make_tuple( - sstp_tmp_rh.begin(), - sstp_tmp_rv.begin(), - Tp.begin() - )), - detail::common__theta_dry__p() - ); - - if(turb_cond) - perparticle_advance_rw2(RH_max, Tp, - pressure_iter, - thrust::make_transform_iterator( - thrust::make_zip_iterator(thrust::make_tuple( - pressure_iter, - sstp_tmp_rv.begin(), - Tp.begin(), - ssp.begin() - )), - detail::RH_sgs(opts_init.RH_formula) - ) - ); - else - perparticle_advance_rw2 - (RH_max, Tp, - pressure_iter, - thrust::make_transform_iterator( - thrust::make_zip_iterator(thrust::make_tuple( - pressure_iter, - sstp_tmp_rv.begin(), - Tp.begin() - )), - detail::RH(opts_init.RH_formula) - ) - ); - } - else - { - if(turb_cond) - perparticle_advance_rw2(RH_max, Tp, - sstp_tmp_p.begin(), - thrust::make_transform_iterator( - thrust::make_zip_iterator(thrust::make_tuple( - sstp_tmp_p.begin(), - sstp_tmp_rv.begin(), - Tp.begin(), - ssp.begin() - )), - detail::RH_sgs(opts_init.RH_formula) - ) - ); - else - perparticle_advance_rw2(RH_max, Tp, - sstp_tmp_p.begin(), - thrust::make_transform_iterator( - thrust::make_zip_iterator(thrust::make_tuple( - sstp_tmp_p.begin(), - sstp_tmp_rv.begin(), - Tp.begin() - )), - detail::RH(opts_init.RH_formula) - ) - ); - } - } - }; -}; 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 c1c03db5..00000000 --- a/src/impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp +++ /dev/null @@ -1,331 +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_sgs(RH_formula)( - thrust::make_tuple(sstp_tmp_p, sstp_tmp_rv, Tp, ssp) - ) : - detail::RH(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) - { - 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 < rc2 && (rw2 + sstp_cond * drw2) > rc2 ) || - ( rw2 > rc2 && (rw2 + sstp_cond * drw2) < rc2 ) ) - { - sstp_cond = sstp_cond_act; - first_cond_step_done_in_adaptation = false; - } - } - 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) - } - } - - 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 - ) - ); - } - }; -}; diff --git a/src/impl/condensation/common/apply_perparticle_sgs_supersat.ipp b/src/impl/condensation_deposition/common/apply_perparticle_sgs_supersat.ipp similarity index 100% rename from src/impl/condensation/common/apply_perparticle_sgs_supersat.ipp rename to src/impl/condensation_deposition/common/apply_perparticle_sgs_supersat.ipp diff --git a/src/impl/condensation_deposition/common/cond_depo_common.ipp b/src/impl/condensation_deposition/common/cond_depo_common.ipp new file mode 100644 index 00000000..362bfd96 --- /dev/null +++ b/src/impl/condensation_deposition/common/cond_depo_common.ipp @@ -0,0 +1,80 @@ +// 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 massdiff2drv + { + real_t mlt; + int n_dims; + + BOOST_GPU_ENABLED + massdiff2drv(const real_t &mlt, const int &n_dims): + mlt(mlt), n_dims(n_dims) {} + + BOOST_GPU_ENABLED + real_t operator()(const real_t &massdiff, const thrust::tuple &tpl) noexcept + { + if(n_dims > 0) + return mlt * massdiff * thrust::get<1>(tpl) / thrust::get<0>(tpl) / thrust::get<2>(tpl); + else // for parcel setup use 1/rhod instead of dv, dv will be updated in hskpng_Tpr in async + return mlt * massdiff * thrust::get<1>(tpl); + } + }; + + template + struct rw2torwX + { + BOOST_GPU_ENABLED + real_t operator()(const real_t &rw2) noexcept + { + #if !defined(__NVCC__) + using std::pow; + #endif + return pow(rw2, real_t(power) / real_t(2)); + } + }; + + template + struct rw2torwX + { + BOOST_GPU_ENABLED + real_t operator()(const real_t &rw2) noexcept + { + #if !defined(__NVCC__) + using std::sqrt; + #endif + return rw2 * sqrt(rw2); + } + }; + + template + struct rw2torwX + { + BOOST_GPU_ENABLED + real_t operator()(const real_t &rw2) noexcept + { + return rw2; + } + }; + }; + }; +}; diff --git a/src/impl/condensation/common/particles_impl_cond_common.ipp b/src/impl/condensation_deposition/common/particles_impl_cond_common.ipp similarity index 58% rename from src/impl/condensation/common/particles_impl_cond_common.ipp rename to src/impl/condensation_deposition/common/particles_impl_cond_common.ipp index 49b61bc9..bfd4a4fd 100644 --- a/src/impl/condensation/common/particles_impl_cond_common.ipp +++ b/src/impl/condensation_deposition/common/particles_impl_cond_common.ipp @@ -19,63 +19,7 @@ namespace libcloudphxx namespace lgrngn { namespace detail - { - template - struct rw3diff2drv - { - real_t mlt; - int n_dims; - - BOOST_GPU_ENABLED - rw3diff2drv(const real_t &mlt, const int &n_dims): - mlt(mlt), n_dims(n_dims) {} - - BOOST_GPU_ENABLED - real_t operator()(const real_t &rw3diff, const thrust::tuple &tpl) noexcept - { - if(n_dims > 0) - return mlt * rw3diff * thrust::get<1>(tpl) / thrust::get<0>(tpl) / thrust::get<2>(tpl); - else // for parcel setup use 1/rhod instead of dv, dv will be updated in hskpng_Tpr in async - return mlt * rw3diff * thrust::get<1>(tpl); - } - }; - - template - struct rw2torwX - { - BOOST_GPU_ENABLED - real_t operator()(const real_t &rw2) noexcept - { - #if !defined(__NVCC__) - using std::pow; - #endif - return pow(rw2, real_t(power) / real_t(2)); - } - }; - - template - struct rw2torwX - { - BOOST_GPU_ENABLED - real_t operator()(const real_t &rw2) noexcept - { - #if !defined(__NVCC__) - using std::sqrt; - #endif - return rw2 * sqrt(rw2); - } - }; - - template - struct rw2torwX - { - BOOST_GPU_ENABLED - real_t operator()(const real_t &rw2) noexcept - { - return rw2; - } - }; - + { template struct advance_rw2_minfun { @@ -197,7 +141,13 @@ namespace libcloudphxx #endif // Skip ice particles - if (rw2_old <= 0) return rw2_old; + if (rw2_old <= 0) + { + if constexpr (apply) + return rw2_old; + else + return real_t(0); + } auto& tpl_in = thrust::get<0>(tpl); const advance_rw2_minfun f(dt, rw2_old, tpl, RH_max); @@ -336,142 +286,6 @@ namespace libcloudphxx return rw2_new - rw2_old; } }; - - template - struct advance_rw2_minfun_ice - { - const quantity r2_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_rw2_minfun_ice( - const real_t &dt, - const real_t &rw2, - const thrust::tuple, real_t, real_t> &tpl, - const real_t &RH_max - ) : - dt(dt * si::seconds), - r2_old(rw2 * 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> drw2_dt(const quantity &rw2) 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 rw = sqrt(real_t(rw2 / si::square_metres)) * si::metres; - const quantity rw3 = rw * rw * rw;; - - const quantity - Re = common::ventil::Re(vt, rw, 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 / rw) * (Sh(Sc, Re) / 2); - - const quantity - K = K_0() * beta(lambda_K / rw) * (Nu(Pr, Re) / 2); - - return real_t(2) * rdrdt_i( - D, - K, - rhod * rv, - T, - p, - RH_i > RH_max ? RH_max : RH_i - ); - } - - BOOST_GPU_ENABLED - real_t operator()(const real_t &rw2_unitless) const - { - const quantity rw2 = rw2_unitless * si::square_metres; - return (r2_old + dt * drw2_dt(rw2) - rw2) / si::square_metres; - } - }; - - - template - struct advance_ice_ac - { - const real_t dt, RH_max; - - advance_ice_ac(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 <= 0 || c_old <= 0) - return ac_old; - - advance_rw2_minfun_ice f_a(dt, a_old * a_old, tpl, RH_max); - advance_rw2_minfun_ice 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; - - // forward Euler for simplicity - const real_t a_new = max(a_old + dt * da_dt, real_t(1e-9)); - const real_t c_new = max(c_old + dt * dc_dt, real_t(1e-9)); - - return thrust::make_tuple(a_new, c_new); - } - }; - }; }; }; diff --git a/src/impl/condensation_deposition/common/particles_impl_depo_common.ipp b/src/impl/condensation_deposition/common/particles_impl_depo_common.ipp new file mode 100644 index 00000000..0460484f --- /dev/null +++ b/src/impl/condensation_deposition/common/particles_impl_depo_common.ipp @@ -0,0 +1,169 @@ +// 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_rw2_minfun_ice + { + const quantity r2_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_rw2_minfun_ice( + const real_t &dt, + const real_t &rw2, + const thrust::tuple, real_t, real_t> &tpl, + const real_t &RH_max + ) : + dt(dt * si::seconds), + r2_old(rw2 * 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> drw2_dt(const quantity &rw2) 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 rw = sqrt(real_t(rw2 / si::square_metres)) * si::metres; + const quantity rw3 = rw * rw * rw;; + + const quantity + Re = common::ventil::Re(vt, rw, 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 / rw) * (Sh(Sc, Re) / 2); + + const quantity + K = K_0() * beta(lambda_K / rw) * (Nu(Pr, Re) / 2); + + return real_t(2) * rdrdt_i( + D, + K, + rhod * rv, + T, + p, + RH_i > RH_max ? RH_max : RH_i + ); + } + + BOOST_GPU_ENABLED + real_t operator()(const real_t &rw2_unitless) const + { + const quantity rw2 = rw2_unitless * si::square_metres; + return (r2_old + dt * drw2_dt(rw2) - rw2) / si::square_metres; + } + }; + + + template + struct advance_ice_ac + { + const real_t dt, RH_max; + + advance_ice_ac(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 <= 0 || c_old <= 0) + { + if constexpr (apply) + return ac_old; + else + return thrust::make_tuple(real_t(0), real_t(0)); + } + + // TODO: growth of ice_a and _c as in Shima 2020 + advance_rw2_minfun_ice f_a(dt, a_old * a_old, tpl, RH_max); + advance_rw2_minfun_ice 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; + + // forward Euler for simplicity, TODO: implicit Euler / toms748 as in condensation + const real_t a_new = max(a_old + dt * da_dt, real_t(1e-9)); + const real_t c_new = max(c_old + dt * dc_dt, real_t(1e-9)); + + if constexpr (apply) + return thrust::make_tuple(a_new, c_new); + else + return thrust::make_tuple(a_new - a_old, c_new - c_old); + } + }; + + }; + }; +}; diff --git a/src/impl/condensation/common/sstp_save.ipp b/src/impl/condensation_deposition/common/sstp_save.ipp similarity index 100% rename from src/impl/condensation/common/sstp_save.ipp rename to src/impl/condensation_deposition/common/sstp_save.ipp diff --git a/src/impl/condensation/percell/particles_impl_cond.ipp b/src/impl/condensation_deposition/percell/particles_impl_cond.ipp similarity index 100% rename from src/impl/condensation/percell/particles_impl_cond.ipp rename to src/impl/condensation_deposition/percell/particles_impl_cond.ipp diff --git a/src/impl/condensation_deposition/percell/particles_impl_depo.ipp b/src/impl/condensation_deposition/percell/particles_impl_depo.ipp new file mode 100644 index 00000000..ccbccd3e --- /dev/null +++ b/src/impl/condensation_deposition/percell/particles_impl_depo.ipp @@ -0,0 +1,152 @@ +// vim:filetype=cpp +/** @file + * @copyright University of Warsaw + * @section LICENSE + * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) + */ + + +namespace libcloudphxx +{ + namespace lgrngn + { + template + void particles_t::impl::ice_dep( + const real_t &dt, + const real_t &RH_max, + const bool turb_cond, + const int step + ) { + + namespace arg = thrust::placeholders; + + thrust_device::vector &lambda_D(lambda_D_gp->get()); + thrust_device::vector &lambda_K(lambda_K_gp->get()); + + hskpng_sort(); + + // Vector to store 3rd moment + if(step == 0) + reset_guardp(ice_mass_percell_gp, tmp_device_real_cell); + thrust_device::vector &ice_mass_percell = ice_mass_percell_gp->get(); + + if(step == 0) + { + if(count_n!=n_cell) // NOTE: we rely that moms_calc was called recently (e.g. in save_liq_ice_content_before_change) + thrust::fill(ice_mass_percell.begin(), ice_mass_percell.end(), real_t(0.)); + } + else // copy ice_mass from previous step + { + reset_guardp(d_ice_mass_percell_gp, tmp_device_real_cell); + thrust_device::vector &d_ice_mass_percell = d_ice_mass_percell_gp->get(); + // drv = -ice_mass precond + thrust::transform( + ice_mass_percell.begin(), ice_mass_percell.end(), + d_ice_mass_percell.begin(), + thrust::negate() + ); + } + + auto hlpr_zip_iter = thrust::make_zip_iterator(thrust::make_tuple( + thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), + thrust::make_permutation_iterator(rv.begin(), ijk.begin()), + thrust::make_permutation_iterator(T.begin(), ijk.begin()), + thrust::make_permutation_iterator(eta.begin(), ijk.begin()), + rd3.begin(), + kpa.begin(), + vt.begin(), + thrust::make_permutation_iterator(lambda_D.begin(), ijk.begin()), + thrust::make_permutation_iterator(lambda_K.begin(), ijk.begin()) + )); + + auto ice_ac_iter = thrust::make_zip_iterator(thrust::make_tuple( + ice_a.begin(), + ice_c.begin() + )); + + if(turb_cond) // with SGS supersaturation + { + auto RH_i_plus_ssp_g = tmp_device_real_part.get_guard(); + thrust_device::vector &RH_i_plus_ssp = RH_i_plus_ssp_g.get(); + thrust::transform( + ssp.begin(), ssp.end(), + thrust::make_permutation_iterator(RH_i.begin(), ijk.begin()), + RH_i_plus_ssp.begin(), + arg::_1 + arg::_2 + ); + + thrust::transform( + ice_ac_iter, + ice_ac_iter + n_part, + thrust::make_zip_iterator( + thrust::make_tuple( + hlpr_zip_iter, + thrust::make_permutation_iterator(p.begin(), ijk.begin()), + RH_i_plus_ssp.begin() + ) + ), + ice_ac_iter, + detail::advance_ice_ac(dt / sstp_cond, RH_max) + ); + } + else // no SGS supersaturation + { + thrust::transform( + ice_ac_iter, + ice_ac_iter + n_part, + thrust::make_zip_iterator( + thrust::make_tuple( + hlpr_zip_iter, + thrust::make_permutation_iterator(p.begin(), ijk.begin()), + thrust::make_permutation_iterator(RH_i.begin(), ijk.begin()) + ) + ), + ice_ac_iter, + detail::advance_ice_ac(dt / sstp_cond, RH_max) + ); + } + nancheck(ice_a, "ice_a after deposition (no sub-steps)"); + nancheck(ice_c, "ice_c after deposition (no sub-steps)"); + + // Compute per-cell 3rd moment of ice after deposition. It is stored in count_mom + moms_gt0(ice_a.begin()); // choose ice particles (ice_a>0) + moms_calc(thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple(ice_a.begin(), ice_c.begin(), ice_rho.begin())), + detail::ice_mass() + ), + real_t(1)); + nancheck_range(count_mom.begin(), count_mom.begin() + count_n, "count_mom (3rd ice moment) after deposition"); + + // Adding the ice volume after deposition to d_ice_mass_percell + if(step < sstp_cond - 1) + { + thrust_device::vector &d_ice_mass_percell = d_ice_mass_percell_gp->get(); + thrust::copy( + count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg + thrust::make_permutation_iterator(ice_mass_percell.begin(), count_ijk.begin()) // output + ); + + // adding the third moment after deposition to ice_mass + thrust::transform( + ice_mass_percell.begin(), ice_mass_percell.end(), + d_ice_mass_percell.begin(), + d_ice_mass_percell.begin(), + thrust::plus() + ); + } + else // last step, calculate change in 3rd moment and update th and rv + { + thrust_device::vector &d_ice_mass_percell = d_ice_mass_percell_gp->get(); + + thrust::transform( + count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg + thrust::make_permutation_iterator(d_ice_mass_percell.begin(), count_ijk.begin()), // input - 2nd arg + thrust::make_permutation_iterator(d_ice_mass_percell.begin(), count_ijk.begin()), // output + thrust::plus() + ); + ice_mass_percell_gp.reset(); // destroy guard to tmp array that stored ice_mass + } + + } + }; + } \ No newline at end of file diff --git a/src/impl/condensation/percell/sstp_percell_step.ipp b/src/impl/condensation_deposition/percell/sstp_percell_step.ipp similarity index 100% rename from src/impl/condensation/percell/sstp_percell_step.ipp rename to src/impl/condensation_deposition/percell/sstp_percell_step.ipp diff --git a/src/impl/condensation/perparticle/acquire_arrays_for_perparticle_sstp.ipp b/src/impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp similarity index 59% rename from src/impl/condensation/perparticle/acquire_arrays_for_perparticle_sstp.ipp rename to src/impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp index 1f890168..7e612960 100644 --- a/src/impl/condensation/perparticle/acquire_arrays_for_perparticle_sstp.ipp +++ b/src/impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp @@ -4,7 +4,7 @@ namespace libcloudphxx namespace lgrngn { template - void particles_t::impl::acquire_arrays_for_perparticle_sstp() + void particles_t::impl::acquire_arrays_for_perparticle_sstp(const bool cond, const bool depo) { reset_guardp(sstp_dlt_rv_gp, tmp_device_real_part); reset_guardp(sstp_dlt_th_gp, tmp_device_real_part); @@ -14,14 +14,22 @@ namespace libcloudphxx if(!sstp_cond_exact_nomix_adaptive) { - reset_guardp(rwX_gp, tmp_device_real_part); - reset_guardp(drwX_gp, tmp_device_real_part); reset_guardp(Tp_gp, tmp_device_real_part); + if(cond) + { + reset_guardp(rwX_gp, tmp_device_real_part); + reset_guardp(drwX_gp, tmp_device_real_part); + } + if(depo) + { + reset_guardp(ice_mass_gp, tmp_device_real_part); + reset_guardp(d_ice_mass_gp, tmp_device_real_part); + } } if(opts_init.adaptive_sstp_cond) { - reset_guardp(perparticle_sstp_cond_gp, tmp_device_n_part); + reset_guardp(perparticle_sstp_cond_gp, tmp_device_n_part); } } }; diff --git a/src/impl/condensation_deposition/perparticle/add_perparticle_ice_mass_to_d_ice_mass.ipp b/src/impl/condensation_deposition/perparticle/add_perparticle_ice_mass_to_d_ice_mass.ipp new file mode 100644 index 00000000..4affa91d --- /dev/null +++ b/src/impl/condensation_deposition/perparticle/add_perparticle_ice_mass_to_d_ice_mass.ipp @@ -0,0 +1,46 @@ +namespace libcloudphxx +{ + namespace lgrngn + { + template + void particles_t::impl::add_perparticle_ice_mass_to_d_ice_mass(const bool store_ice_mass) + { + thrust_device::vector &d_ice_mass = d_ice_mass_gp->get(); + + auto ice_mass_it = + thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple(ice_a.begin(), ice_c.begin(), ice_rho.begin())), + detail::ice_mass() + ); + + if(store_ice_mass) + { + thrust_device::vector &ice_mass = ice_mass_gp->get(); + + thrust::transform( + ice_mass_it, ice_mass_it + n_part, + ice_mass.begin(), + thrust::identity() + ); + + thrust::transform( + ice_mass.begin(), ice_mass.end(), + d_ice_mass.begin(), + d_ice_mass.begin(), + thrust::plus() + ); + + } + else + { + thrust::transform( + ice_mass_it, ice_mass_it + n_part, + d_ice_mass.begin(), + d_ice_mass.begin(), + thrust::plus() + ); + } + } + }; +}; + diff --git a/src/impl/condensation/perparticle/add_perparticle_rwX_to_drwX.ipp b/src/impl/condensation_deposition/perparticle/add_perparticle_rwX_to_drwX.ipp similarity index 100% rename from src/impl/condensation/perparticle/add_perparticle_rwX_to_drwX.ipp rename to src/impl/condensation_deposition/perparticle/add_perparticle_rwX_to_drwX.ipp diff --git a/src/impl/condensation/perparticle/apply_noncond_perparticle_sstp_delta.ipp b/src/impl/condensation_deposition/perparticle/apply_noncond_perparticle_sstp_delta.ipp similarity index 100% rename from src/impl/condensation/perparticle/apply_noncond_perparticle_sstp_delta.ipp rename to src/impl/condensation_deposition/perparticle/apply_noncond_perparticle_sstp_delta.ipp diff --git a/src/impl/condensation/perparticle/apply_perparticle_cond_change_to_percell_rv_and_th.ipp b/src/impl/condensation_deposition/perparticle/apply_perparticle_cond_change_to_percell_rv_and_th.ipp similarity index 100% rename from src/impl/condensation/perparticle/apply_perparticle_cond_change_to_percell_rv_and_th.ipp rename to src/impl/condensation_deposition/perparticle/apply_perparticle_cond_change_to_percell_rv_and_th.ipp 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 new file mode 100644 index 00000000..e58378ac --- /dev/null +++ b/src/impl/condensation_deposition/perparticle/apply_perparticle_drw3_or_d_ice_mass_to_perparticle_rv_and_th.ipp @@ -0,0 +1,130 @@ +namespace libcloudphxx +{ + namespace lgrngn + { + template + void particles_t::impl::apply_perparticle_drw3_or_d_ice_mass_to_perparticle_rv_and_th(const bool rw3_changed, const bool ice_mass_changed) + { + assert(rw3_changed || ice_mass_changed); + + namespace arg = thrust::placeholders; + thrust_device::vector &Tp = Tp_gp->get(); + + if(rw3_changed) + { + // get change in liquid water mass + thrust_device::vector &drw3 = drwX_gp->get(); + // ice crystals should have drw3==0 + thrust::transform( + drw3.begin(), drw3.end(), + drw3.begin(), // in-place! + arg::_1 * (common::moist_air::rho_w() / si::kilograms * si::cubic_metres) + * real_t(4./3) * pi() + ); + + // calculate rv change from the change in liquid water mass + thrust::transform( + drw3.begin(), drw3.end(), + thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_rh.begin(), + n.begin(), + thrust::make_permutation_iterator(dv.begin(), ijk.begin()) + )), + drw3.begin(), + detail::massdiff2drv(real_t(-1), n_dims) + ); + + // apply rv change to rv (of this particle only or of all particles, depending on sstp_cond_mix) + if(opts_init.sstp_cond_mix) + update_pstate(sstp_tmp_rv, drw3); + else + thrust::transform( + drw3.begin(), drw3.end(), + sstp_tmp_rv.begin(), + sstp_tmp_rv.begin(), + thrust::plus() + ); + + // calculate dth from drv + thrust::transform( + thrust::make_zip_iterator(thrust::make_tuple( + drw3.begin(), + Tp.begin(), + sstp_tmp_th.begin() + )), + thrust::make_zip_iterator(thrust::make_tuple( + drw3.end(), + Tp.end(), + sstp_tmp_th.end() + )), + drw3.begin(), + detail::dth() + ); + } + + // 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_gp->get(); + // calculate rv change from the change in solid water mass + thrust::transform( + d_ice_mass.begin(), d_ice_mass.end(), + thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_rh.begin(), + n.begin(), + thrust::make_permutation_iterator(dv.begin(), ijk.begin()) + )), + d_ice_mass.begin(), + detail::massdiff2drv(real_t(-1), n_dims) + ); + + // apply rv change to rv (of this particle only or of all particles, depending on sstp_cond_mix) + if(opts_init.sstp_cond_mix) + update_pstate(sstp_tmp_rv, d_ice_mass); + else + thrust::transform( + d_ice_mass.begin(), d_ice_mass.end(), + sstp_tmp_rv.begin(), + sstp_tmp_rv.begin(), + thrust::plus() + ); + + // calculate dth from the change in ice mass + thrust::transform( + thrust::make_zip_iterator(thrust::make_tuple( + d_ice_mass.begin(), + Tp.begin(), + sstp_tmp_th.begin() + )), + thrust::make_zip_iterator(thrust::make_tuple( + d_ice_mass.end(), + Tp.end(), + sstp_tmp_th.end() + )), + d_ice_mass.begin(), + detail::dth_dep() + ); + } + + 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_gp->get().begin(), + dth.begin(), + thrust::plus() + ); + + // apply dth from condensation and deposition + if(opts_init.sstp_cond_mix) + update_pstate(sstp_tmp_th, dth); + else + thrust::transform( + dth.begin(), dth.end(), + sstp_tmp_th.begin(), + sstp_tmp_th.begin(), + thrust::plus() + ); + } + }; +}; diff --git a/src/impl/condensation_deposition/perparticle/calc_perparticle_T.ipp b/src/impl/condensation_deposition/perparticle/calc_perparticle_T.ipp new file mode 100644 index 00000000..d7ed71fe --- /dev/null +++ b/src/impl/condensation_deposition/perparticle/calc_perparticle_T.ipp @@ -0,0 +1,37 @@ +namespace libcloudphxx +{ + namespace lgrngn + { + template + void particles_t::impl::calc_perparticle_T( + ) { + + namespace arg = thrust::placeholders; + + thrust_device::vector &Tp = Tp_gp->get(); + + // calculate perparticle temperature + if(opts_init.th_dry) + { + thrust::transform( + sstp_tmp_th.begin(), sstp_tmp_th.end(), + sstp_tmp_rh.begin(), + Tp.begin(), + detail::common__theta_dry__T_rhod() + ); + } + else + { + thrust::transform( + sstp_tmp_th.begin(), sstp_tmp_th.end(), + thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_rv.begin(), + sstp_tmp_p.begin() + )), + Tp.begin(), + detail::common__theta_std__T_p() + ); + } + } + }; +}; diff --git a/src/impl/condensation/perparticle/calculate_noncond_perparticle_sstp_delta.ipp b/src/impl/condensation_deposition/perparticle/calculate_noncond_perparticle_sstp_delta.ipp similarity index 100% rename from src/impl/condensation/perparticle/calculate_noncond_perparticle_sstp_delta.ipp rename to src/impl/condensation_deposition/perparticle/calculate_noncond_perparticle_sstp_delta.ipp diff --git a/src/impl/condensation/perparticle/init_perparticle_sstp.ipp b/src/impl/condensation_deposition/perparticle/init_perparticle_sstp.ipp similarity index 100% rename from src/impl/condensation/perparticle/init_perparticle_sstp.ipp rename to src/impl/condensation_deposition/perparticle/init_perparticle_sstp.ipp diff --git a/src/impl/condensation/perparticle/perparticle_advance_rw2.ipp b/src/impl/condensation_deposition/perparticle/perparticle_advance_hlpr.ipp similarity index 50% rename from src/impl/condensation/perparticle/perparticle_advance_rw2.ipp rename to src/impl/condensation_deposition/perparticle/perparticle_advance_hlpr.ipp index 8386c1c5..f0484a73 100644 --- a/src/impl/condensation/perparticle/perparticle_advance_rw2.ipp +++ b/src/impl/condensation_deposition/perparticle/perparticle_advance_hlpr.ipp @@ -3,8 +3,8 @@ namespace libcloudphxx namespace lgrngn { template - template - void particles_t::impl::perparticle_advance_rw2( + template + void particles_t::impl::perparticle_advance_hlpr( const real_t &RH_max, const thrust_device::vector &Tp, const pres_iter &pi, @@ -30,16 +30,38 @@ namespace libcloudphxx thrust::make_permutation_iterator(lambda_K.begin(), ijk.begin()) )); - thrust::transform( - rw2.begin(), rw2.end(), - thrust::make_zip_iterator(thrust::make_tuple( - hlpr_zip_iter, - pi, - rhi - )), - rw2.begin(), - detail::advance_rw2(dt / sstp_cond, RH_max, config.eps_tolerance, config.cond_mlt, config.n_iter) - ); + if constexpr(!ice) + { + thrust::transform( + rw2.begin(), rw2.end(), + thrust::make_zip_iterator(thrust::make_tuple( + hlpr_zip_iter, + pi, + rhi + )), + rw2.begin(), + detail::advance_rw2(dt / sstp_cond, RH_max, config.eps_tolerance, config.cond_mlt, config.n_iter) + ); + } + else // ice + { + auto hlpr_zip_iter_ice = thrust::make_zip_iterator(thrust::make_tuple( + ice_a.begin(), + ice_c.begin() + )); + + thrust::transform( + hlpr_zip_iter_ice, + hlpr_zip_iter_ice + n_part, + thrust::make_zip_iterator(thrust::make_tuple( + hlpr_zip_iter, + pi, + rhi + )), + hlpr_zip_iter_ice, + detail::advance_ice_ac(dt / sstp_cond, RH_max) + ); + } } }; }; diff --git a/src/impl/condensation_deposition/perparticle/perparticle_advance_size.ipp b/src/impl/condensation_deposition/perparticle/perparticle_advance_size.ipp new file mode 100644 index 00000000..618087e9 --- /dev/null +++ b/src/impl/condensation_deposition/perparticle/perparticle_advance_size.ipp @@ -0,0 +1,197 @@ +namespace libcloudphxx +{ + namespace lgrngn + { + namespace detail + { + template + struct RH_hlpr + { + std::conditional_t, RH > resolved_RH; + + BOOST_GPU_ENABLED + RH_hlpr(RH_formula_t RH_formula): + resolved_RH(RH_formula) + {} + + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple &tpl) noexcept + { + assert(sgs); + return resolved_RH(thrust::make_tuple(thrust::get<0>(tpl), thrust::get<1>(tpl), thrust::get<2>(tpl))) + thrust::get<3>(tpl); + } + + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple &tpl) noexcept + { + assert(!sgs); + return resolved_RH(thrust::make_tuple(thrust::get<0>(tpl), thrust::get<1>(tpl), thrust::get<2>(tpl))); + } + }; + }; + + + template + template + void particles_t::impl::perparticle_advance_size( + const real_t &RH_max, + const bool turb_cond + ) { + + namespace arg = thrust::placeholders; + + thrust_device::vector &Tp = Tp_gp->get(); + + // std::conditional_t)> advance_func = ice ? advance_ice_a_c> : advance_hlpr>; + + // advance rw2 + if(!opts_init.const_p) + { + auto pressure_iter = thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_rh.begin(), + sstp_tmp_rv.begin(), + Tp.begin() + )), + detail::common__theta_dry__p() + ); + + if(turb_cond) + perparticle_advance_hlpr(RH_max, Tp, + pressure_iter, + thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple( + pressure_iter, + sstp_tmp_rv.begin(), + Tp.begin(), + ssp.begin() + )), + detail::RH_hlpr(opts_init.RH_formula) + ) + ); + else + perparticle_advance_hlpr + (RH_max, Tp, + pressure_iter, + thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple( + pressure_iter, + sstp_tmp_rv.begin(), + Tp.begin() + )), + detail::RH_hlpr(opts_init.RH_formula) + // detail::RH(opts_init.RH_formula) + ) + ); + } + else + { + if(turb_cond) + perparticle_advance_hlpr(RH_max, Tp, + sstp_tmp_p.begin(), + thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_p.begin(), + sstp_tmp_rv.begin(), + Tp.begin(), + ssp.begin() + )), + detail::RH_hlpr(opts_init.RH_formula) + ) + ); + else + perparticle_advance_hlpr(RH_max, Tp, + sstp_tmp_p.begin(), + thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple( + sstp_tmp_p.begin(), + sstp_tmp_rv.begin(), + Tp.begin() + )), + detail::RH_hlpr(opts_init.RH_formula) + ) + ); + } + } + + + // template + // void particles_t::impl::cond_perparticle_advance_hlpr( + // const real_t &RH_max, + // const bool turb_cond + // ) { + + // namespace arg = thrust::placeholders; + + // thrust_device::vector &Tp = Tp_gp->get(); + + // // advance rw2 + // if(!opts_init.const_p) + // { + // auto pressure_iter = thrust::make_transform_iterator( + // thrust::make_zip_iterator(thrust::make_tuple( + // sstp_tmp_rh.begin(), + // sstp_tmp_rv.begin(), + // Tp.begin() + // )), + // detail::common__theta_dry__p() + // ); + + // if(turb_cond) + // perparticle_advance_hlpr(RH_max, Tp, + // pressure_iter, + // thrust::make_transform_iterator( + // thrust::make_zip_iterator(thrust::make_tuple( + // pressure_iter, + // sstp_tmp_rv.begin(), + // Tp.begin(), + // ssp.begin() + // )), + // detail::RH_sgs(opts_init.RH_formula) + // ) + // ); + // else + // perparticle_advance_hlpr + // (RH_max, Tp, + // pressure_iter, + // thrust::make_transform_iterator( + // thrust::make_zip_iterator(thrust::make_tuple( + // pressure_iter, + // sstp_tmp_rv.begin(), + // Tp.begin() + // )), + // detail::RH(opts_init.RH_formula) + // ) + // ); + // } + // else + // { + // if(turb_cond) + // perparticle_advance_hlpr(RH_max, Tp, + // sstp_tmp_p.begin(), + // thrust::make_transform_iterator( + // thrust::make_zip_iterator(thrust::make_tuple( + // sstp_tmp_p.begin(), + // sstp_tmp_rv.begin(), + // Tp.begin(), + // ssp.begin() + // )), + // detail::RH_sgs(opts_init.RH_formula) + // ) + // ); + // else + // perparticle_advance_hlpr(RH_max, Tp, + // sstp_tmp_p.begin(), + // thrust::make_transform_iterator( + // thrust::make_zip_iterator(thrust::make_tuple( + // sstp_tmp_p.begin(), + // sstp_tmp_rv.begin(), + // Tp.begin() + // )), + // detail::RH(opts_init.RH_formula) + // ) + // ); + // } + // } + }; +}; diff --git a/src/impl/condensation_deposition/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp b/src/impl/condensation_deposition/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp new file mode 100644 index 00000000..c5891365 --- /dev/null +++ b/src/impl/condensation_deposition/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp @@ -0,0 +1,477 @@ +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, ice_switch, cond, depo; + 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), + ice_switch(opts_init.ice_switch), + cond(opts.cond), + depo(opts.depo) + {} + + 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)); + thrust::tuple ice_ac; // = thrust::make_tuple(0, 0); // only used if ice_switch is true, but we need to define it here to be able to use it in lambdas + real_t ice_rho; + + if(rw2 > 0 && !cond) return; // skip liquid droplets if condensation is turned off + + if(ice_switch) + { + real_t ice_a = thrust::get<7>(thrust::get<1>(tpl)); + if(ice_a > 0 && !depo) return; // skip ice particles if deposition is turned off + real_t ice_c = thrust::get<8>(thrust::get<1>(tpl)); + ice_ac = thrust::make_tuple(ice_a, ice_c); + ice_rho = thrust::get<9>(thrust::get<1>(tpl)); + } + + const bool ice = rw2<=0; // flag to indicate if the particle is ice or liquid; + + real_t drw2, Tp, RH; + thrust::tuple dice_ac; + + // 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 + { + if(!ice) + 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) + ); + else + 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; + thrust::tuple dice_ac_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 dX convergence for increasing number of substeps + for(int sstp_cond_try = 1; sstp_cond_try <= sstp_cond_max; 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(); + if(!ice) + (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 + ) + ); + else // ice + (sstp_cond_try == 1 ? dice_ac : dice_ac_new) = + detail::advance_ice_ac(dt / sstp_cond_try, RH_max)( + ice_ac, + 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 + { + // TODO: get rid of _max and of actviation adaptation? + if(!ice) + { + 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; + } + else // ice + { + const real_t da_dt_rel = cuda::std::abs(thrust::get<0>(dice_ac_new) * 2 - thrust::get<0>(dice_ac)) / (thrust::get<0>(dice_ac) + 1e-20); + const real_t dc_dt_rel = cuda::std::abs(thrust::get<1>(dice_ac_new) * 2 - thrust::get<1>(dice_ac)) / (thrust::get<1>(dice_ac) + 1e-20); + if(da_dt_rel <= sstp_cond_adapt_drw2_eps && da_dt_rel <= sstp_cond_adapt_drw2_eps + && da_dt_rel <= sstp_cond_adapt_drw2_max && dc_dt_rel <= sstp_cond_adapt_drw2_max) + { + 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; + } + dice_ac = dice_ac_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 < rc2 && (rw2 + sstp_cond * drw2) > rc2 ) || + ( rw2 > rc2 && (rw2 + sstp_cond * drw2) < rc2 ) ) + { + 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; + + // actual condensation substepping + if(!ice) + { + 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(); + + 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::massdiff2drv( + - 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; + } + } + else // actual deposition substepping + { + auto _advance_ice_ac = detail::advance_ice_ac(dt / sstp_cond, RH_max); + real_t &ice_mass = drw2; // if its ice, drw2 is not used at all + real_t dice_mass; + + auto ice_mass_fctr = detail::ice_mass(); + + for(int step = 0; step < sstp_cond; ++step) + { + dice_mass = step > 0 ? -ice_mass : -ice_mass_fctr(thrust::make_tuple(thrust::get<0>(ice_ac), thrust::get<1>(ice_ac), ice_rho)); + + if(first_cond_step_done_in_adaptation && step == 0) + { + thrust::get<0>(ice_ac) += thrust::get<0>(dice_ac); + thrust::get<1>(ice_ac) += thrust::get<1>(dice_ac); + } + else + { + _apply_noncond_perparticle_sstp_delta(delta_fraction_applied); + _calc_Tp(); + _calc_sstp_tmp_p(); + _calc_RH(); + + ice_ac = _advance_ice_ac( + ice_ac, + 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) + { + ice_mass = ice_mass_fctr(thrust::make_tuple(thrust::get<0>(ice_ac), thrust::get<1>(ice_ac), ice_rho)); + dice_mass += ice_mass; + } + else + dice_mass += ice_mass_fctr(thrust::make_tuple(thrust::get<0>(ice_ac), thrust::get<1>(ice_ac), ice_rho)); + + dice_mass = detail::massdiff2drv(real_t(1), n_dims) + (dice_mass, + thrust::make_tuple(sstp_tmp_rh, n, dv) + ); + + sstp_tmp_rv += dice_mass; + + dice_mass = detail::dth_dep()( + thrust::make_tuple(dice_mass, Tp, sstp_tmp_th) + ); + + sstp_tmp_th += dice_mass; + } + } + + // 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; + if(!ice) + thrust::get<1>(thrust::get<1>(tpl)) = rw2; + else // ice + { + thrust::get<7>(thrust::get<1>(tpl)) = thrust::get<0>(ice_ac); + thrust::get<8>(thrust::get<1>(tpl)) = thrust::get<1>(ice_ac); + } + } + }; + }; + + 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(), + ice_a.begin(), + ice_c.begin(), + ice_rho.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 + ) + ); + } + }; +}; diff --git a/src/impl/condensation/perparticle/release_arrays_for_perparticle_sstp.ipp b/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp similarity index 68% rename from src/impl/condensation/perparticle/release_arrays_for_perparticle_sstp.ipp rename to src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp index ea8d4de6..9dbf2401 100644 --- a/src/impl/condensation/perparticle/release_arrays_for_perparticle_sstp.ipp +++ b/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp @@ -4,7 +4,7 @@ namespace libcloudphxx namespace lgrngn { template - void particles_t::impl::release_arrays_for_perparticle_sstp() + void particles_t::impl::release_arrays_for_perparticle_sstp(const bool cond, const bool depo) { sstp_dlt_rv_gp.reset(); sstp_dlt_th_gp.reset(); @@ -14,9 +14,18 @@ namespace libcloudphxx if(!sstp_cond_exact_nomix_adaptive) { - rwX_gp.reset(); - drwX_gp.reset(); Tp_gp.reset(); + + if(cond) + { + rwX_gp.reset(); + drwX_gp.reset(); + } + if(depo) + { + ice_mass_gp.reset(); + d_ice_mass_gp.reset(); + } } if(opts_init.adaptive_sstp_cond) @@ -25,4 +34,4 @@ namespace libcloudphxx } } }; -}; \ No newline at end of file +}; diff --git a/src/impl/condensation_deposition/perparticle/set_perparticle_d_ice_mass_to_minus_ice_mass.ipp b/src/impl/condensation_deposition/perparticle/set_perparticle_d_ice_mass_to_minus_ice_mass.ipp new file mode 100644 index 00000000..59bebe2a --- /dev/null +++ b/src/impl/condensation_deposition/perparticle/set_perparticle_d_ice_mass_to_minus_ice_mass.ipp @@ -0,0 +1,35 @@ +namespace libcloudphxx +{ + namespace lgrngn + { + template + void particles_t::impl::set_perparticle_d_ice_mass_to_minus_ice_mass(const bool use_stored_ice_mass) + { + thrust_device::vector &d_ice_mass = d_ice_mass_gp->get(); + if(!use_stored_ice_mass) + { + auto ice_mass_it = + thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple(ice_a.begin(), ice_c.begin(), ice_rho.begin())), + detail::ice_mass() + ); + + thrust::transform( + ice_mass_it, ice_mass_it + n_part, + d_ice_mass.begin(), + thrust::negate() + ); + } + else + { + thrust_device::vector &ice_mass = ice_mass_gp->get(); + thrust::transform( + ice_mass.begin(), ice_mass.end(), + d_ice_mass.begin(), + thrust::negate() + ); + } + } + }; +}; + diff --git a/src/impl/condensation/perparticle/set_perparticle_drwX_to_minus_rwX.ipp b/src/impl/condensation_deposition/perparticle/set_perparticle_drwX_to_minus_rwX.ipp similarity index 100% rename from src/impl/condensation/perparticle/set_perparticle_drwX_to_minus_rwX.ipp rename to src/impl/condensation_deposition/perparticle/set_perparticle_drwX_to_minus_rwX.ipp diff --git a/src/impl/housekeeping/particles_impl_hskpng_Tpr.ipp b/src/impl/housekeeping/particles_impl_hskpng_Tpr.ipp index 97c735ee..417b6510 100644 --- a/src/impl/housekeeping/particles_impl_hskpng_Tpr.ipp +++ b/src/impl/housekeeping/particles_impl_hskpng_Tpr.ipp @@ -193,8 +193,10 @@ namespace libcloudphxx switch (RH_formula) { case RH_formula_t::pv_cc: + case RH_formula_t::pv_tet: return RHi_pv_cc(thrust::get<0>(tpl), thrust::get<1>(tpl), thrust::get<2>(tpl)); case RH_formula_t::rv_cc: + case RH_formula_t::rv_tet: return RHi_rv_cc(thrust::get<0>(tpl), thrust::get<1>(tpl), thrust::get<2>(tpl)); // NOTE: no Tetens formulas for ice default: @@ -269,18 +271,14 @@ namespace libcloudphxx if(opts_init.ice_switch) { - // RH_i - thrust::transform( - zip_it_t(thrust::make_tuple(p.begin(), rv.begin(), T.begin())), // input - begin - zip_it_t(thrust::make_tuple(p.end(), rv.end(), T.end() )), // input - end - RH_i.begin(), // output - detail::RH_i( - opts_init.RH_formula == RH_formula_t::pv_tet ? RH_formula_t::pv_cc : - opts_init.RH_formula == RH_formula_t::rv_tet ? RH_formula_t::rv_cc : - opts_init.RH_formula - ) - ); - } + // RH_i + thrust::transform( + zip_it_t(thrust::make_tuple(p.begin(), rv.begin(), T.begin())), // input - begin + zip_it_t(thrust::make_tuple(p.end(), rv.end(), T.end() )), // input - end + RH_i.begin(), // output + detail::RH_i(opts_init.RH_formula) + ); + } } // dynamic viscosity diff --git a/src/impl/ice/particles_impl_ice_dep.ipp b/src/impl/ice/particles_impl_ice_dep.ipp deleted file mode 100644 index 36e69615..00000000 --- a/src/impl/ice/particles_impl_ice_dep.ipp +++ /dev/null @@ -1,134 +0,0 @@ -// vim:filetype=cpp -/** @file - * @copyright University of Warsaw - * @section LICENSE - * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) - */ - - -namespace libcloudphxx -{ - namespace lgrngn - { - template - void particles_t::impl::ice_dep( - const real_t &dt, - const real_t &RH_max, - const int step - ) { - - namespace arg = thrust::placeholders; - - thrust_device::vector &lambda_D(lambda_D_gp->get()); - thrust_device::vector &lambda_K(lambda_K_gp->get()); - - hskpng_sort(); - - // Vector to store 3rd moment - if(step == 0) - reset_guardp(ice_mass_gp, tmp_device_real_cell); - thrust_device::vector &ice_mass = ice_mass_gp->get(); - - if(step == 0) - { - if(count_n!=n_cell) // NOTE: we rely that moms_calc was called recently (e.g. in save_liq_ice_content_before_change) - thrust::fill(ice_mass.begin(), ice_mass.end(), real_t(0.)); - } - else // copy ice_mass from previous step - { - reset_guardp(d_ice_mass_gp, tmp_device_real_cell); - thrust_device::vector &d_ice_mass = d_ice_mass_gp->get(); - // drv = -ice_mass precond - thrust::transform( - ice_mass.begin(), ice_mass.end(), - d_ice_mass.begin(), - thrust::negate() - ); - } - - auto hlpr_zip_iter = thrust::make_zip_iterator(thrust::make_tuple( - thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), - thrust::make_permutation_iterator(rv.begin(), ijk.begin()), - thrust::make_permutation_iterator(T.begin(), ijk.begin()), - thrust::make_permutation_iterator(eta.begin(), ijk.begin()), - rd3.begin(), - kpa.begin(), - vt.begin(), - thrust::make_permutation_iterator(lambda_D.begin(), ijk.begin()), - thrust::make_permutation_iterator(lambda_K.begin(), ijk.begin()) - )); - - // deposition for ice crystals - thrust::transform( - thrust::make_zip_iterator( - thrust::make_tuple( - ice_a.begin(), - ice_c.begin() - ) - ), - thrust::make_zip_iterator( - thrust::make_tuple( - ice_a.end(), - ice_c.end() - ) - ), - thrust::make_zip_iterator( - thrust::make_tuple( - hlpr_zip_iter, - thrust::make_permutation_iterator(p.begin(), ijk.begin()), - thrust::make_permutation_iterator(RH_i.begin(), ijk.begin()) - ) - ), - thrust::make_zip_iterator( - thrust::make_tuple( - ice_a.begin(), - ice_c.begin() - ) - ), - detail::advance_ice_ac(dt / sstp_cond, RH_max) - ); - nancheck(ice_a, "ice_a after deposition (no sub-steps"); - nancheck(ice_c, "ice_c after deposition (no sub-steps"); - - // Compute per-cell 3rd moment of ice after deposition. It is stored in count_mom - moms_gt0(ice_a.begin()); // choose ice particles (ice_a>0) - moms_calc(thrust::make_transform_iterator( - thrust::make_zip_iterator(thrust::make_tuple(ice_a.begin(), ice_c.begin(), ice_rho.begin())), - detail::ice_mass() - ), - real_t(1)); - nancheck_range(count_mom.begin(), count_mom.begin() + count_n, "count_mom (3rd ice moment) after deposition"); - - // Adding the ice volume after deposition to d_ice_mass - if(step < sstp_cond - 1) - { - thrust_device::vector &d_ice_mass = d_ice_mass_gp->get(); - thrust::copy( - count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg - thrust::make_permutation_iterator(ice_mass.begin(), count_ijk.begin()) // output - ); - - // adding the third moment after deposition to ice_mass - thrust::transform( - ice_mass.begin(), ice_mass.end(), - d_ice_mass.begin(), - d_ice_mass.begin(), - thrust::plus() - ); - } - else // last step, calculate change in 3rd moment and update th and rv - { - thrust_device::vector &d_ice_mass = d_ice_mass_gp->get(); - - thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg - thrust::make_permutation_iterator(d_ice_mass.begin(), count_ijk.begin()), // input - 2nd arg - thrust::make_permutation_iterator(d_ice_mass.begin(), count_ijk.begin()), // output - thrust::plus() - ); - ice_mass_gp.reset(); // destroy guard to tmp array that stored ice_mass - } - - } - }; - } \ No newline at end of file diff --git a/src/impl/initialization/particles_impl_init_sanity_check.ipp b/src/impl/initialization/particles_impl_init_sanity_check.ipp index ec1768ff..e1976522 100644 --- a/src/impl/initialization/particles_impl_init_sanity_check.ipp +++ b/src/impl/initialization/particles_impl_init_sanity_check.ipp @@ -168,10 +168,10 @@ namespace libcloudphxx throw std::runtime_error("libcloudph++: relaxation does not work with ice."); if(opts_init.src_type==src_t::matching) // because we dont account for ice/water when matching and initializing aerosols from this type of source throw std::runtime_error("libcloudph++: 'matching' source type does not work with ice."); - if(opts_init.turb_cond_switch) // because we dont want to add SGS RH to RH_i + if(opts_init.turb_cond_switch) // current code would simply add SGS RH to RH_i, but this is not correct (?); how are they related? throw std::runtime_error("libcloudph++: SGS condensation does not work with ice."); - if(opts_init.exact_sstp_cond) - throw std::runtime_error("libcloudph++: deposition works only with per-cell substepping"); + // if(opts_init.adaptive_sstp_cond) + // throw std::runtime_error("libcloudph++: deposition does not work with adaptive substepping (opts_init.ice_switch==true implies that deposition may be modeled)"); } } }; diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 681bf689..d2a0d684 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -232,8 +232,10 @@ namespace libcloudphxx // drw3_gp, Tp_gp, // rw3_gp, + ice_mass_gp, + ice_mass_percell_gp, d_ice_mass_gp, - ice_mass_gp; + d_ice_mass_percell_gp; std::unique_ptr< typename tmp_vector_pool>::guard @@ -693,11 +695,14 @@ namespace libcloudphxx void sedi(const real_t &dt); void subs(const real_t &dt); + void calc_perparticle_T(); + // condensation methods void cond(const real_t &dt, const real_t &RH_max, const bool turb_cond, const int step); - void cond_perparticle_advance_rw2(const real_t &RH_max, const bool turb_cond); - template - void perparticle_advance_rw2(const real_t &RH_max, const thrust_device::vector &Tp, const pres_iter &pi, const RH_iter &rhi); + template + void perparticle_advance_size(const real_t &RH_max, const bool turb_cond); + template + void perparticle_advance_hlpr(const real_t &RH_max, const thrust_device::vector &Tp, const pres_iter &pi, const RH_iter &rhi); void perparticle_nomixing_adaptive_sstp_cond(const opts_t &); void save_liq_ice_content_before_change(); void calc_liq_ice_content_change(); @@ -705,8 +710,10 @@ namespace libcloudphxx void set_perparticle_drwX_to_minus_rwX(const bool use_stored_rw3); template void add_perparticle_rwX_to_drwX(const bool store_rw3); + void set_perparticle_d_ice_mass_to_minus_ice_mass(const bool use_stored_ice_mass); + void add_perparticle_ice_mass_to_d_ice_mass(const bool store_ice_mass); - void apply_perparticle_drw3_to_perparticle_rv_and_th(); + void apply_perparticle_drw3_or_d_ice_mass_to_perparticle_rv_and_th(const bool rw3_changed, const bool ice_mass_changed); void apply_perparticle_cond_change_to_percell_rv_and_th(); void update_th_rv(); // ======= @@ -715,7 +722,7 @@ namespace libcloudphxx // void cond_dm3_helper(); // void cond(const real_t &dt, const real_t &RH_max, const bool turb_cond, const int step); // void cond_sstp(const real_t &dt, const real_t &RH_max, const bool turb_cond, const int step); - void ice_dep(const real_t &dt, const real_t &RH_max, const int step); + void ice_dep(const real_t &dt, const real_t &RH_max, const bool turb_cond, const int step); // template // void cond_sstp_hlpr(const real_t &dt, const real_t &RH_max, const thrust_device::vector &Tp, const pres_iter &pi, const RH_iter &rhi); void update_th_freezing(thrust_device::vector &); @@ -750,8 +757,8 @@ namespace libcloudphxx void rlx_dry_distros(const real_t); // substepping methods - void acquire_arrays_for_perparticle_sstp(); - void release_arrays_for_perparticle_sstp(); + void acquire_arrays_for_perparticle_sstp(const bool cond, const bool depo); + void release_arrays_for_perparticle_sstp(const bool cond, const bool depo); void calculate_noncond_perparticle_sstp_delta(); void apply_noncond_perparticle_sstp_delta(); void apply_perparticle_sgs_supersat(); diff --git a/src/particles.tpp b/src/particles.tpp index c7e98d1e..95f42d48 100644 --- a/src/particles.tpp +++ b/src/particles.tpp @@ -122,23 +122,29 @@ #include "impl/advection/particles_impl_adve.ipp" #include "impl/advection/particles_impl_turb_adve.ipp" -#include "impl/condensation/common/apply_perparticle_sgs_supersat.ipp" -#include "impl/condensation/common/particles_impl_cond_common.ipp" -#include "impl/condensation/common/sstp_save.ipp" -#include "impl/condensation/percell/sstp_percell_step.ipp" -#include "impl/condensation/percell/particles_impl_cond.ipp" -#include "impl/condensation/perparticle/init_perparticle_sstp.ipp" -#include "impl/condensation/perparticle/acquire_arrays_for_perparticle_sstp.ipp" -#include "impl/condensation/perparticle/release_arrays_for_perparticle_sstp.ipp" -#include "impl/condensation/perparticle/calculate_noncond_perparticle_sstp_delta.ipp" -#include "impl/condensation/perparticle/apply_noncond_perparticle_sstp_delta.ipp" -#include "impl/condensation/perparticle/perparticle_advance_rw2.ipp" -#include "impl/condensation/perparticle/cond_perparticle_advance_rw2.ipp" -#include "impl/condensation/perparticle/apply_perparticle_drw3_to_perparticle_rv_and_th.ipp" -#include "impl/condensation/perparticle/apply_perparticle_cond_change_to_percell_rv_and_th.ipp" -#include "impl/condensation/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp" -#include "impl/condensation/perparticle/set_perparticle_drwX_to_minus_rwX.ipp" -#include "impl/condensation/perparticle/add_perparticle_rwX_to_drwX.ipp" +#include "impl/condensation_deposition/common/apply_perparticle_sgs_supersat.ipp" +#include "impl/condensation_deposition/common/cond_depo_common.ipp" +#include "impl/condensation_deposition/common/particles_impl_cond_common.ipp" +#include "impl/condensation_deposition/common/particles_impl_depo_common.ipp" +#include "impl/condensation_deposition/common/sstp_save.ipp" +#include "impl/condensation_deposition/percell/sstp_percell_step.ipp" +#include "impl/condensation_deposition/percell/particles_impl_cond.ipp" +#include "impl/condensation_deposition/percell/particles_impl_depo.ipp" +#include "impl/condensation_deposition/perparticle/init_perparticle_sstp.ipp" +#include "impl/condensation_deposition/perparticle/calc_perparticle_T.ipp" +#include "impl/condensation_deposition/perparticle/acquire_arrays_for_perparticle_sstp.ipp" +#include "impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp" +#include "impl/condensation_deposition/perparticle/calculate_noncond_perparticle_sstp_delta.ipp" +#include "impl/condensation_deposition/perparticle/apply_noncond_perparticle_sstp_delta.ipp" +#include "impl/condensation_deposition/perparticle/perparticle_advance_size.ipp" +#include "impl/condensation_deposition/perparticle/perparticle_advance_hlpr.ipp" +#include "impl/condensation_deposition/perparticle/apply_perparticle_drw3_or_d_ice_mass_to_perparticle_rv_and_th.ipp" +#include "impl/condensation_deposition/perparticle/apply_perparticle_cond_change_to_percell_rv_and_th.ipp" +#include "impl/condensation_deposition/perparticle/perparticle_nomixing_adaptive_sstp_cond.ipp" +#include "impl/condensation_deposition/perparticle/set_perparticle_drwX_to_minus_rwX.ipp" +#include "impl/condensation_deposition/perparticle/add_perparticle_rwX_to_drwX.ipp" +#include "impl/condensation_deposition/perparticle/set_perparticle_d_ice_mass_to_minus_ice_mass.ipp" +#include "impl/condensation_deposition/perparticle/add_perparticle_ice_mass_to_d_ice_mass.ipp" #include "impl/sedimentation/particles_impl_sedi.ipp" @@ -163,7 +169,6 @@ #include "impl/sources_and_relaxation_of_SDs/particles_impl_rlx.ipp" #include "impl/sources_and_relaxation_of_SDs/particles_impl_rlx_dry_distros.ipp" -#include "impl/ice/particles_impl_ice_dep.ipp" #include "impl/ice/particles_impl_ice_nucl_melt.ipp" diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 5e5a9f23..0120196b 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -172,6 +172,9 @@ namespace libcloudphxx if(opts.turb_cond && !pimpl->opts_init.turb_cond_switch) throw std::runtime_error("libcloudph++: turb_cond_swtich=False, but turb_cond==True"); + if(opts.depo && !pimpl->opts_init.ice_switch) + throw std::runtime_error("libcloudph++: opts.depo==true, but opts_init.ice_switch=False"); + pimpl->should_now_run_cond = false; // dt defined in opts_init can be overriden by dt in opts @@ -185,7 +188,7 @@ namespace libcloudphxx pimpl->ice_nucl_melt(pimpl->dt); // condensation/evaporation - if (opts.cond) + if (opts.cond || opts.depo) { // prerequisite pimpl->hskpng_sort(); @@ -200,11 +203,11 @@ namespace libcloudphxx { if(!pimpl->opts_init.sstp_cond_mix) { - pimpl->save_liq_ice_content_before_change(); // in drw_mom3_gp and d_ice_mass_gp + pimpl->save_liq_ice_content_before_change(); // in drw_mom3_gp and d_ice_mass_percell_gp pimpl->n_filtered_gp.reset(); // n_filtered, acquired and filled in rw_mom3_ante_change, not needed anymore } - pimpl->acquire_arrays_for_perparticle_sstp(); // sstp_dlt_rv_gp, etc. ; as in sstp_percell_step_exact() + pimpl->acquire_arrays_for_perparticle_sstp(opts.cond, opts.depo); // 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? // adaptive per-particle substepping @@ -223,17 +226,31 @@ namespace libcloudphxx pimpl->apply_noncond_perparticle_sstp_delta(); if(opts.turb_cond) pimpl->apply_perparticle_sgs_supersat(); - - pimpl->template set_perparticle_drwX_to_minus_rwX<3>(/*use_stored_rw3=*/ step>0); - pimpl->cond_perparticle_advance_rw2(opts.RH_max, opts.turb_cond); - pimpl->template add_perparticle_rwX_to_drwX<3>(/*store_rw3=*/ step < pimpl->sstp_cond - 1); - pimpl->apply_perparticle_drw3_to_perparticle_rv_and_th(); + pimpl->calc_perparticle_T(); // TODO: implement it by moving code from cond_perparticle_advance_rw2 + + // condensation + if(opts.cond) + { + pimpl->template set_perparticle_drwX_to_minus_rwX<3>(/*use_stored_rw3=*/ step>0); + pimpl->template perparticle_advance_size(opts.RH_max, opts.turb_cond); + pimpl->template add_perparticle_rwX_to_drwX<3>(/*store_rw3=*/ step < pimpl->sstp_cond - 1); + } + + // deposition + if(opts.depo) + { + pimpl->set_perparticle_d_ice_mass_to_minus_ice_mass(step>0); + pimpl->template perparticle_advance_size(opts.RH_max, opts.turb_cond); + pimpl->add_perparticle_ice_mass_to_d_ice_mass(step < pimpl->sstp_cond - 1); + } + + pimpl->apply_perparticle_drw3_or_d_ice_mass_to_perparticle_rv_and_th(opts.cond, opts.depo); } } - pimpl->release_arrays_for_perparticle_sstp(); + pimpl->release_arrays_for_perparticle_sstp(opts.cond, opts.depo); pimpl->apply_perparticle_cond_change_to_percell_rv_and_th(); } - else // apply per-cell sstp logic, always with mixing (sstp_cond_mix==true required) + else // apply per-cell sstp logic, always with mixing { for (int step = 0; step < pimpl->sstp_cond; ++step) { @@ -242,13 +259,13 @@ namespace libcloudphxx pimpl->apply_perparticle_sgs_supersat(); pimpl->hskpng_Tpr(); if(step == 0) - pimpl->save_liq_ice_content_before_change(); // in drw_mom3_gp and d_ice_mass_gp + pimpl->save_liq_ice_content_before_change(); // in drw_mom3_gp and d_ice_mass_percell_gp pimpl->cond(pimpl->dt, opts.RH_max, opts.turb_cond, step); // pimpl->cond(pimpl->dt / pimpl->sstp_cond, opts.RH_max, opts.turb_cond, step); - if (pimpl->opts_init.ice_switch) + if (opts.depo) //pimpl->opts_init.ice_switch) { // pimpl->ice_dep(pimpl->dt / pimpl->sstp_cond, opts.RH_max, step); - pimpl->ice_dep(pimpl->dt, opts.RH_max, step); + pimpl->ice_dep(pimpl->dt, opts.RH_max, opts.turb_cond, step); } pimpl->update_th_rv(); }