diff --git a/CMakeLists.txt b/CMakeLists.txt index eb8bde863..1a2b6504a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -277,7 +277,7 @@ if (CMAKE_CUDA_COMPILER) target_compile_options(cloudphxx_lgrngn PRIVATE $<$: -O3 -use_fast_math -Xcompiler=-Ofast,-march=native,-fopenmp>) # Debug mode cuda flags elseif(CMAKE_BUILD_TYPE STREQUAL "Debug") - target_compile_options(cloudphxx_lgrngn PRIVATE $<$: -g -DTHRUST_DEBUG -lineinfo -Xcompiler=-fopenmp,-g,-Og,-rdynamic>) + target_compile_options(cloudphxx_lgrngn PRIVATE $<$: -g -DTHRUST_DEBUG -lineinfo -Xcompiler=-fopenmp,-g,-Og,-rdynamic,-Wfatal-errors>) # Release mode cuda flags else() target_compile_options(cloudphxx_lgrngn PRIVATE $<$: -DNDEBUG -O3 -use_fast_math -Xcompiler=-Ofast,-march=native,-DNDEBUG,-fopenmp>) @@ -400,7 +400,7 @@ add_subdirectory(include) if(NOT LIBCLOUDPHXX_DISABLE_BINDINGS) # some tests require bindings, so we disable them too enable_testing() - add_subdirectory(bindings) + add_subdirectory(bindings) add_subdirectory(tests) endif() diff --git a/bindings/python/common.hpp b/bindings/python/common.hpp index 58c174c18..c361d6f46 100644 --- a/bindings/python/common.hpp +++ b/bindings/python/common.hpp @@ -58,6 +58,12 @@ namespace libcloudphxx return cmn::const_cp::p_vs(T * si::kelvins) / si::pascals; } + template + real_t p_vsi(const real_t &T) + { + return cmn::const_cp::p_vsi(T * si::kelvins) / si::pascals; + } + template real_t p_vs_tet(const real_t &T) { @@ -69,6 +75,12 @@ namespace libcloudphxx { return cmn::const_cp::r_vs(T * si::kelvins, p * si::pascals); } + + template + real_t r_vsi(const real_t &T, const real_t &p) + { + return cmn::const_cp::r_vsi(T * si::kelvins, p * si::pascals); + } template real_t l_v(const real_t &T) @@ -76,6 +88,18 @@ namespace libcloudphxx return cmn::const_cp::l_v(T * si::kelvins) / si::joules * si::kilograms; } + template + real_t l_s(const real_t &T) + { + return cmn::const_cp::l_s(T * si::kelvins) / si::joules * si::kilograms; + } + + template + real_t l_f(const real_t &T) + { + return cmn::const_cp::l_f(T * si::kelvins) / si::joules * si::kilograms; + } + template real_t T(const real_t &th, const real_t &rhod) { diff --git a/bindings/python/lgrngn.hpp b/bindings/python/lgrngn.hpp index 693d13663..9223cd6e3 100644 --- a/bindings/python/lgrngn.hpp +++ b/bindings/python/lgrngn.hpp @@ -222,7 +222,7 @@ namespace libcloudphxx map_t map = arg->diag_puddle(); bp::dict dict; for(auto& x : map) - dict[static_cast(x.first)] = x.second; + dict[cmn::output_names.at(x.first)] = x.second; return dict; } @@ -234,45 +234,71 @@ namespace libcloudphxx return *arg->opts_init; } + // set dry distros from a dict with (kappa, rd_insol) as key template - void set_dd( // dry_distro + void set_dd( lgr::opts_init_t *arg, - const bp::dict &kappa_func) + const bp::dict &kappa_func) // a dict keyed by (kappa, rd_insol) { arg->dry_distros.clear(); for (int i = 0; i < len(kappa_func.keys()); ++i) + { + bp::tuple key = bp::extract(kappa_func.keys()[i]); + const real_t kappa = bp::extract(key[0]); + const real_t rd_insol = bp::extract(key[1]); arg->dry_distros.emplace( - bp::extract(kappa_func.keys()[i]), + libcloudphxx::lgrngn::kappa_rd_insol_t{kappa, rd_insol}, std::make_shared>(kappa_func.values()[i]) ); + } } + // src_dry_distros moved from opts_init to opts template void set_sdd( // src_dry_distro - lgr::opts_init_t *arg, - const bp::dict &kappa_func) + lgr::opts_t *arg, + const bp::dict &kappa_func) // a dict keyed by (kappa, rd_insol) { arg->src_dry_distros.clear(); for (int i = 0; i < len(kappa_func.keys()); ++i) + { + bp::tuple key = bp::extract(kappa_func.keys()[i]); + bp::tuple val = bp::extract(kappa_func.values()[i]); + const real_t kappa = bp::extract(key[0]); + const real_t rd_insol = bp::extract(key[1]); + + const int sd_conc = bp::extract(val[1]); + const int supstp = bp::extract(val[2]); arg->src_dry_distros.emplace( - bp::extract(kappa_func.keys()[i]), - std::make_shared>(kappa_func.values()[i]) + libcloudphxx::lgrngn::kappa_rd_insol_t{kappa, rd_insol}, + std::make_tuple(std::static_pointer_cast>( + std::make_shared>(val[0])), + sd_conc, supstp) ); + } } + // set dry sizes from a dict with (kappa, rd_insol) as key template - void set_ds( // dry_sizes + void set_ds( lgr::opts_init_t *arg, - const bp::dict &kappa_func + const bp::dict &kappa_func // a dict keyed by (kappa, rd_insol) ) { arg->dry_sizes.clear(); if(len(kappa_func.keys()) == 0) return; - // loop over kappas + // loop over kappas and rd_insol for (int j = 0; j < len(kappa_func.keys()); ++j) { + + // extract the key tuple (kappa, rd_insol) + const bp::tuple key = bp::extract(kappa_func.keys()[j]); + const real_t kappa = bp::extract(key[0]); + const real_t rd_insol = bp::extract(key[1]); + + // extract size : {conc, count} dict for this (kappa, rd_insol) const bp::dict size_conc = bp::extract(kappa_func.values()[j]); std::map> size_conc_map; @@ -283,43 +309,50 @@ namespace libcloudphxx assert(len(conc_count_list) == 2); const real_t conc = bp::extract(conc_count_list[0]); const int count = bp::extract (conc_count_list[1]); - size_conc_map[bp::extract(size_conc.keys()[i])] = std::make_pair(conc, count); + size_conc_map[bp::extract(size_conc.keys()[i])] = std::make_pair(conc, count); } - const real_t kappa = bp::extract(kappa_func.keys()[j]); - arg->dry_sizes[kappa] = size_conc_map; + arg->dry_sizes[libcloudphxx::lgrngn::kappa_rd_insol_t{kappa, rd_insol}] = size_conc_map; } } + // src_dry_sizes moved from opts_init to opts template void set_sds( // src_dry_sizes - lgr::opts_init_t *arg, - const bp::dict &kappa_func + lgr::opts_t *arg, + const bp::dict &kappa_func // a dict keyed by (kappa, rd_insol) ) { arg->src_dry_sizes.clear(); if(len(kappa_func.keys()) == 0) return; - // loop over kappas + // loop over kappas and rd_insol for (int j = 0; j < len(kappa_func.keys()); ++j) { + // extract the key tuple (kappa, rd_insol) + const bp::tuple key = bp::extract(kappa_func.keys()[j]); + const real_t kappa = bp::extract(key[0]); + const real_t rd_insol = bp::extract(key[1]); + + // extract size : {conc, count, supstp_src} dict for this (kappa, rd_insol) const bp::dict size_conc = bp::extract(kappa_func.values()[j]); - std::map> size_conc_map; + std::map> size_conc_map; - // turn the size : {conc, count} dict into a size : {conc, count} map + // turn the size : {conc, count, supstp_src} dict into a size : {conc, count, supstp_src} map for (int i = 0; i < len(size_conc.keys()); ++i) { const bp::list conc_count_list = bp::extract(size_conc.values()[i]); - assert(len(conc_count_list) == 2); + assert(len(conc_count_list) == 3); const real_t conc = bp::extract(conc_count_list[0]); const int count = bp::extract (conc_count_list[1]); - size_conc_map[bp::extract(size_conc.keys()[i])] = std::make_pair(conc, count); + const int supstp = bp::extract (conc_count_list[2]); + size_conc_map[bp::extract(size_conc.keys()[i])] = std::make_tuple(conc, count, supstp); } - const real_t kappa = bp::extract(kappa_func.keys()[j]); - arg->src_dry_sizes[kappa] = size_conc_map; + arg->src_dry_sizes[libcloudphxx::lgrngn::kappa_rd_insol_t{kappa, rd_insol}] = size_conc_map; } } + template void set_rdd( // rlx_dry_distros lgr::opts_init_t *arg, diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index edb763121..fca30884e 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -63,6 +63,7 @@ BOOST_PYTHON_MODULE(libcloudphxx) bp::scope().attr("eps") = (real_t) (cmn::moist_air::eps()); bp::scope().attr("rho_stp") = (real_t) (cmn::earth::rho_stp() * si::cubic_metres / si::kilograms); bp::scope().attr("rho_w") = (real_t) (cmn::moist_air::rho_w() * si::cubic_metres / si::kilograms); + bp::scope().attr("rho_i") = (real_t) (cmn::moist_air::rho_i() * si::cubic_metres / si::kilograms); //molar mass of trace gases bp::scope().attr("M_SO2") = (real_t) (cmn::molar_mass::M_SO2() * si::moles / si::kilograms); bp::scope().attr("M_H2O2") = (real_t) (cmn::molar_mass::M_H2O2() * si::moles / si::kilograms); @@ -131,6 +132,7 @@ BOOST_PYTHON_MODULE(libcloudphxx) bp::def("p_v", &common::p_v); bp::def("p_vs", &common::p_vs); 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("T", &common::T); @@ -278,13 +280,14 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("turb_adve", &lgr::opts_t::turb_adve) .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("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) ; bp::class_>("opts_init_t") .add_property("dry_distros", &lgrngn::get_dd, &lgrngn::set_dd) .add_property("dry_sizes", &lgrngn::get_ds, &lgrngn::set_ds) - .add_property("src_dry_distros", &lgrngn::get_sdd, &lgrngn::set_sdd) - .add_property("src_dry_sizes", &lgrngn::get_ds, &lgrngn::set_sds) .add_property("rlx_dry_distros", &lgrngn::get_rdd, &lgrngn::set_rdd) .def_readwrite("nx", &lgr::opts_init_t::nx) .def_readwrite("ny", &lgr::opts_init_t::ny) @@ -321,7 +324,6 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("sstp_cond", &lgr::opts_init_t::sstp_cond) .def_readwrite("sstp_coal", &lgr::opts_init_t::sstp_coal) .def_readwrite("sstp_chem", &lgr::opts_init_t::sstp_chem) - .def_readwrite("supstp_src", &lgr::opts_init_t::supstp_src) .def_readwrite("supstp_rlx", &lgr::opts_init_t::supstp_rlx) .def_readwrite("kernel", &lgr::opts_init_t::kernel) .def_readwrite("adve_scheme", &lgr::opts_init_t::adve_scheme) @@ -329,7 +331,6 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("sd_conc_large_tail", &lgr::opts_init_t::sd_conc_large_tail) .def_readwrite("aerosol_independent_of_rhod", &lgr::opts_init_t::aerosol_independent_of_rhod) .def_readwrite("sd_const_multi", &lgr::opts_init_t::sd_const_multi) - .def_readwrite("src_sd_conc", &lgr::opts_init_t::src_sd_conc) .def_readwrite("rlx_bins", &lgr::opts_init_t::rlx_bins) .def_readwrite("rlx_sd_per_bin", &lgr::opts_init_t::rlx_sd_per_bin) .def_readwrite("rlx_timescale", &lgr::opts_init_t::rlx_timescale) @@ -347,6 +348,8 @@ BOOST_PYTHON_MODULE(libcloudphxx) .add_property("SGS_mix_len", &lgrngn::get_SGS_mix_len, &lgrngn::set_SGS_mix_len) .add_property("kernel_parameters", &lgrngn::get_kp, &lgrngn::set_kp) .def_readwrite("variable_dt_switch", &lgr::opts_init_t::variable_dt_switch) + .def_readwrite("ice_switch", &lgr::opts_init_t::ice_switch) + .def_readwrite("time_dep_ice_nucl", &lgr::opts_init_t::time_dep_ice_nucl) .def_readwrite("const_p", &lgr::opts_init_t::const_p) .def_readwrite("th_dry", &lgr::opts_init_t::th_dry) ; @@ -409,6 +412,13 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def("diag_chem", &lgr::particles_proto_t::diag_chem) .def("diag_precip_rate", &lgr::particles_proto_t::diag_precip_rate) .def("diag_puddle", &lgrngn::diag_puddle) + .def("diag_ice", &lgr::particles_proto_t::diag_ice) + .def("diag_water", &lgr::particles_proto_t::diag_water) + .def("diag_ice_cons", &lgr::particles_proto_t::diag_ice_cons) + .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_mix_ratio", &lgr::particles_proto_t::diag_ice_mix_ratio) .def("outbuf", &lgrngn::outbuf) .def("get_attr", &lgr::particles_proto_t::get_attr) ; diff --git a/include/libcloudph++/common/ice_nucleation.hpp b/include/libcloudph++/common/ice_nucleation.hpp new file mode 100644 index 000000000..fec85b0e7 --- /dev/null +++ b/include/libcloudph++/common/ice_nucleation.hpp @@ -0,0 +1,146 @@ +#pragma once + +#include "units.hpp" +#include "const_cp.hpp" +#include + +#if defined(__NVCC__) +# include +#endif + +namespace libcloudphxx +{ + namespace common + { + namespace ice_nucleation + { + enum class INP_t {mineral}; // types of ice nucleating particles, TODO: add more types + + // Inverse CDF for singular freezing temperature as defined in eq. 1 in Shima et al., 2020 + template + BOOST_GPU_ENABLED + quantity T_freeze_CDF_inv( + const INP_t& INP_type, // type of ice nucleating particle + const real_t rd2_insol, // radius squared of insoluble particle in m^2 + const real_t rand // random number between [0, 1] + ) { + real_t A = real_t(4) + #if !defined(__NVCC__) + * pi() + #else + * CUDART_PI + #endif + * rd2_insol; // surface area of the insoluble particle + + if (INP_type == INP_t::mineral && A > real_t(1e-20)) + { + return (real_t(273.15) + (real_t(8.934) - log(- log(real_t(1.) - rand) / A) ) / real_t(0.517)) * si::kelvin; + } + else + { + return real_t(235.15) * si::kelvin; // the default freezing temperature is -38 C + } + } + + + template + struct T_freeze_CDF_inv_functor + { + INP_t INP_type; + + T_freeze_CDF_inv_functor(INP_t INP_type) + : INP_type(INP_type) {} + + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple &tpl) const + { + const real_t &rd2_insol = thrust::get<0>(tpl); // from rd2 vector + const real_t &rand = thrust::get<1>(tpl); // from rand vector + + return ice_nucleation::template T_freeze_CDF_inv( + INP_type, + rd2_insol, + rand + ) / si::kelvin; + } + }; + + // Probability of time-dependent freezing, + // heterogeneous as in Arabas et al., 2025 and homogeneous as in Koop & Murray, 2016 + template + BOOST_GPU_ENABLED + real_t p_freeze( + const INP_t& INP_type, // type of ice nucleating particle + const real_t rd2_insol, // radius squared of insoluble particle in m^2 + const real_t rw2, // wet radius squared in m^2 + const real_t T, // temperature in kelvin + const real_t dt // time step in seconds + ) + { + if (rd2_insol > real_t(0)) + { + real_t A = real_t(4) + #if !defined(__NVCC__) + * pi() + #else + * CUDART_PI + #endif + * rd2_insol; // surface area of the insoluble particle + real_t d_aw = real_t(1) - const_cp::p_vsi(T * si::kelvin)/ const_cp::p_vs(T * si::kelvin); // water activity + if (INP_type == INP_t::mineral) + { + real_t J_het = pow(real_t(10), real_t(-1.35) + real_t(22.62) * d_aw) * real_t(1e4); // nucleation rate + return 1 - exp(- J_het * A * dt); + } + else + return real_t(0.); // TODO: other INP types + } + else + { + real_t V = real_t(4./3.) + #if !defined(__NVCC__) + * pi() + #else + * CUDART_PI + #endif + * pow(rw2, real_t(3./2.)); // droplet volume + + real_t dT = T - real_t(273.15); + real_t x = - real_t(3020.684) - real_t(425.921)*pow(dT,real_t(1)) - real_t(25.9779)*pow(dT,real_t(2)) + - real_t(0.868451)*pow(dT,real_t(3)) - real_t(0.0166203)*pow(dT,real_t(4)) + - real_t(0.000171736)*pow(dT,real_t(5)) - real_t(0.000000746953)*pow(dT,real_t(6)); + real_t J_hom = pow(real_t(10), x) * real_t(1e6); // nucleation rate + return 1 - exp(- J_hom * V * dt); + } + } + + + template + struct p_freeze_functor + { + INP_t INP_type; + real_t dt; + + p_freeze_functor(INP_t INP_type, real_t dt) + : INP_type(INP_type), dt(dt) {} + + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple &tpl) const + { + const real_t &rd2_insol = thrust::get<0>(tpl); // radius squared of insoluble particle + const real_t &rw2 = thrust::get<1>(tpl); // wet radius squared + const real_t &T = thrust::get<2>(tpl); // temperature in kelvin + + return ice_nucleation::p_freeze( + INP_type, + rd2_insol, + rw2, + T, + dt + ); + } + }; + + }; + }; +}; diff --git a/include/libcloudph++/common/maxwell-mason.hpp b/include/libcloudph++/common/maxwell-mason.hpp index c86de2371..4145570c0 100644 --- a/include/libcloudph++/common/maxwell-mason.hpp +++ b/include/libcloudph++/common/maxwell-mason.hpp @@ -15,34 +15,71 @@ namespace libcloudphxx quantity::type, real_t> rdrdt( 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, // p_v/p_vs = relative humidity - const quantity a_w, // water activity - const quantity klvntrm // the Kelvin term + const quantity rho_v, // ambient water vapour density + const quantity T, // ambient temperature + const quantity p, // ambient pressure + const quantity RH, // p_v/p_vs = relative humidity + const quantity a_w, // water activity + const quantity klvntrm // the Kelvin term ) { using moist_air::rho_w; + using moist_air::rho_i; using moist_air::R_v; - quantity::type, real_t> + quantity::type, real_t> l_v = const_cp::l_v(T); - return (real_t(1) - a_w * klvntrm / RH) - / rho_w() - / ( - real_t(1) + return (real_t(1) - a_w * klvntrm / RH) + / rho_w() + / ( + real_t(1) / D / rho_v - + - l_v + + + l_v / K / RH / T * (l_v / R_v() / T - real_t(1)) - ) - ; + ) + ; + } + + // for ice + // mass accommodation coeff = 1 + // no solute nor curvature effects + template + BOOST_GPU_ENABLED + quantity::type, real_t> rdrdt_i( + 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 + ) + { + using moist_air::rho_i; + using moist_air::R_v; + + quantity::type, real_t> + l_s = const_cp::l_s(T); + + return (real_t(1) - real_t(1) / RH_i) + / rho_i() + / ( + real_t(1) + / D + / rho_v + + + l_s + / K + / RH_i + / T + * (l_s / R_v() / T - real_t(1)) + ) + ; } }; }; diff --git a/include/libcloudph++/common/output.hpp b/include/libcloudph++/common/output.hpp index 90ef25f15..a15253017 100644 --- a/include/libcloudph++/common/output.hpp +++ b/include/libcloudph++/common/output.hpp @@ -18,7 +18,9 @@ namespace libcloudphxx outliq_vol, outdry_vol, outprtcl_num, - outice_vol + outice_mass, + outliq_num, + outice_num }; const std::map output_names @@ -34,7 +36,9 @@ namespace libcloudphxx {outliq_vol, "liquid_volume"}, {outdry_vol, "dry_volume"}, {outprtcl_num, "particle_number"}, - {outice_vol, "ice_volume"} + {outice_mass, "ice_mass"}, + {outliq_num, "liquid_number"}, + {outice_num, "ice_number"} }; inline output_t get_output_enum(const std::string& name) diff --git a/include/libcloudph++/common/theta_dry.hpp b/include/libcloudph++/common/theta_dry.hpp index cf1d49c8e..49f77c07c 100644 --- a/include/libcloudph++/common/theta_dry.hpp +++ b/include/libcloudph++/common/theta_dry.hpp @@ -54,6 +54,7 @@ namespace libcloudphxx // rho/(1+r) R(r)*(1+r) } + // heat of condensation template BOOST_GPU_ENABLED quantity d_th_d_rv( @@ -63,6 +64,26 @@ namespace libcloudphxx return - th / T * const_cp::l_v(T) / c_pd(); } + // heat of deposition + template + BOOST_GPU_ENABLED + quantity d_th_d_rv_dep( + const quantity &T, + const quantity &th // theta dry!!! + ) { + return - th / T * const_cp::l_s(T) / c_pd(); + } + + // heat of freezing + template + BOOST_GPU_ENABLED + quantity d_th_d_rw_freeze( + const quantity &T, + const quantity &th // theta dry!!! + ) { + return - th / T * const_cp::l_f(T) / c_pd(); + } + template BOOST_GPU_ENABLED quantity std2dry( diff --git a/include/libcloudph++/lgrngn/distro_t.hpp b/include/libcloudph++/lgrngn/distro_t.hpp new file mode 100644 index 000000000..9389b5846 --- /dev/null +++ b/include/libcloudph++/lgrngn/distro_t.hpp @@ -0,0 +1,59 @@ +#pragma once + +namespace libcloudphxx +{ + namespace lgrngn + { + using common::unary_function; + + template + struct kappa_rd_insol_t { + real_t kappa; + real_t rd_insol; + + kappa_rd_insol_t(real_t kappa_, real_t rd_insol_) + : kappa(kappa_), rd_insol(rd_insol_) + {} + + bool operator<(const kappa_rd_insol_t &other) const + { + if (kappa != other.kappa) return kappa < other.kappa; + return rd_insol < other.rd_insol; + } + }; + + // initial dry sizes of aerosol + // defined with a distribution + // uses shared_ptr to make opts_init copyable + template + using dry_distros_t = std::map< + kappa_rd_insol_t, // (kappa, rd_insol) + std::shared_ptr> // n(ln(rd)) @ STP; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true + >; + + // defined with a size-number pair + template + using dry_sizes_t = std::map< + kappa_rd_insol_t, // (kappa, rd_insol) + std::map // STP_concentration [1/m^3], number of SD that represent this radius kappa and concentration + > + >; + + // similar, but for sources of aerosols after initialization + template + using src_dry_distros_t = std::map< + kappa_rd_insol_t, // (kappa, rd_insol) + std::tuple>, int, int> // 1st: n(ln(rd)) @ STP created per second; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true; 2nd: sd_conc for this distribution ; 3rd: supstp for this aerosol (interval in timesteps beween addition of these aerosols) + >; + + // defined with a size-number pair + template + using src_dry_sizes_t = std::map< + kappa_rd_insol_t, // (kappa, rd_insol) + std::map // STP_concentration [1/m^3] created per second, number of SD that represent this radius kappa and concentration, supstp + > + >; + }; +}; diff --git a/include/libcloudph++/lgrngn/opts.hpp b/include/libcloudph++/lgrngn/opts.hpp index 8a32a7e25..8bf6b59f9 100644 --- a/include/libcloudph++/lgrngn/opts.hpp +++ b/include/libcloudph++/lgrngn/opts.hpp @@ -8,6 +8,7 @@ #pragma once #include "extincl.hpp" +#include "distro_t.hpp" #include "../common/chem.hpp" namespace libcloudphxx @@ -19,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; + bool adve, sedi, subs, cond, coal, src, rlx, rcyc, turb_adve, turb_cond, turb_coal, ice_nucl; // RH limit for drop growth real_t RH_max; @@ -31,11 +32,17 @@ namespace libcloudphxx // overriding dt from opts_init real_t dt; + // aerosol source dry sizes defined with a distribution, number of SD and supstp + src_dry_distros_t src_dry_distros; + + // aerosol source dry sizes defined with a size-number pair, number of SD and supstp + src_dry_sizes_t src_dry_sizes; + // ctor with defaults (C++03 compliant) ... 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), + turb_adve(false), turb_cond(false), turb_coal(false), ice_nucl(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 6a8d2232b..ca958b93f 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -13,35 +13,26 @@ #include "advection_scheme.hpp" #include "RH_formula.hpp" #include "ccn_source.hpp" +#include "distro_t.hpp" #include "../common/chem.hpp" +#include "../common/ice_nucleation.hpp" namespace libcloudphxx { namespace lgrngn { using common::unary_function; + using common::ice_nucleation::INP_t; // template struct opts_init_t { - // initial dry sizes of aerosol - // defined with a distribution - // uses shared_ptr to make opts_init copyable - typedef std::unordered_map< - real_t, // kappa - std::shared_ptr> // n(ln(rd)) @ STP; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true - > dry_distros_t; - dry_distros_t dry_distros; - - // defined with a size-number pair - typedef std::map< - real_t, // kappa - std::map // STP_concentration [1/m^3], number of SD that represent this radius kappa and concentration - > - > dry_sizes_t; - dry_sizes_t dry_sizes; + // initial dry sizes of aerosol defined with a distribution + dry_distros_t dry_distros; + + // initial dry sizes of aerosol defined with a size-number pair + dry_sizes_t dry_sizes; // Eulerian component parameters int nx, ny, nz; @@ -97,7 +88,12 @@ namespace libcloudphxx turb_adve_switch, // if true, turbulent motion of SDs is modeled turb_cond_switch, // if true, turbulent condensation of SDs is modeled turb_coal_switch, // if true, turbulent coalescence kernels can be used - exact_sstp_cond; // if true, use per-particle sstp_cond logic, if false, use per-cell + ice_switch, // if true, ice is allowed + exact_sstp_cond, // if true, use per-particle sstp_cond logic, if false, use per-cell + time_dep_ice_nucl; // it true, time-dependent freezing, if false, singular freezing + + // ice nucleating particles type + INP_t inp_type; int sstp_chem; real_t chem_rho; @@ -144,22 +140,9 @@ namespace libcloudphxx // type of CCN source src_t src_type; - // source distro per unit time - dry_distros_t src_dry_distros; - - // number of SDs created from src_dry_distros per cell per source iteration - unsigned long long src_sd_conc; - - // dry sizes of droplets added from the source, STP_concentration created per unit time instead of the STP_concentration - dry_sizes_t src_dry_sizes; - // box in which aerosol from source will be created // will be rounded to cell number - cells are supposed to be uniform real_t src_x0, src_y0, src_z0, src_x1, src_y1, src_z1; - - // timestep interval at which source will be applied - int supstp_src; - // --- aerosol relaxation stuff --- // initial dry sizes of aerosol @@ -212,6 +195,9 @@ namespace libcloudphxx coal_switch(true), // coalescence turned on by default src_type(src_t::off), // source turned off by default rlx_switch(false), + ice_switch(false), + time_dep_ice_nucl(false), + inp_type(INP_t::mineral), exact_sstp_cond(false), turb_cond_switch(false), turb_adve_switch(false), @@ -227,14 +213,12 @@ namespace libcloudphxx dev_count(0), dev_id(-1), n_sd_max(0), - src_sd_conc(0), src_x0(0), src_x1(0), src_y0(0), src_y1(0), src_z0(0), src_z1(0), - supstp_src(1), rlx_bins(0), rlx_timescale(1), rlx_sd_per_bin(0), diff --git a/include/libcloudph++/lgrngn/particles.hpp b/include/libcloudph++/lgrngn/particles.hpp index 41c0e3da4..c0546e54d 100644 --- a/include/libcloudph++/lgrngn/particles.hpp +++ b/include/libcloudph++/lgrngn/particles.hpp @@ -84,8 +84,12 @@ namespace libcloudphxx virtual void diag_RH_ge_Sc() { assert(false); } virtual void diag_dry_rng(const real_t&, const real_t&) { assert(false); } virtual void diag_wet_rng(const real_t&, const real_t&) { assert(false); } + virtual void diag_ice_a_rng(const real_t&, const real_t&) { assert(false); } + virtual void diag_ice_c_rng(const real_t&, const real_t&) { assert(false); } virtual void diag_kappa_rng(const real_t&, const real_t&) { assert(false); } - // The 3 following functions are for consecutive selection of SDs. + virtual void diag_ice() { assert(false); } + virtual void diag_water() { assert(false); } + // The following functions are for consecutive selection of SDs. // It allows the user to select SDs based on multiple characteristics, e.g. wet radius (0.5, 1) and kappa (0.1, 0.2): // diag_wet_rng(0.5, 1); diag_kappa_rng_cons(0.1, 0.2); // NOTE: the call with "cons" needs to be right after the previous call to diag_X_rng! @@ -94,13 +98,21 @@ namespace libcloudphxx // and Boost Python does not work well with virtual member functions that have default arguments virtual void diag_dry_rng_cons(const real_t&, const real_t&) { assert(false); } virtual void diag_wet_rng_cons(const real_t&, const real_t&) { assert(false); } + virtual void diag_ice_a_rng_cons(const real_t&, const real_t&) { assert(false); } + virtual void diag_ice_c_rng_cons(const real_t&, const real_t&) { assert(false); } virtual void diag_kappa_rng_cons(const real_t&, const real_t&) { assert(false); } + virtual void diag_ice_cons() { assert(false); } + virtual void diag_water_cons() { assert(false); } virtual void diag_dry_mom(const int&) { assert(false); } virtual void diag_wet_mom(const int&) { assert(false); } + virtual void diag_ice_a_mom(const int&) { assert(false); } + virtual void diag_ice_c_mom(const int&) { assert(false); } + virtual void diag_ice_mix_ratio() { assert(false); } virtual void diag_wet_mass_dens(const real_t&, const real_t&) { assert(false); } virtual void diag_chem(const enum common::chem::chem_species_t&) { assert(false); } virtual void diag_precip_rate() { assert(false); } + virtual void diag_precip_rate_ice_mass() { assert(false); } virtual void diag_kappa_mom(const int&) { assert(false); } virtual void diag_up_mom(const int&) { assert(false); } virtual void diag_vp_mom(const int&) { assert(false); } @@ -178,12 +190,23 @@ namespace libcloudphxx void diag_RH(); void diag_dry_rng(const real_t &r_mi, const real_t &r_mx); void diag_wet_rng(const real_t &r_mi, const real_t &r_mx); + void diag_ice_a_rng(const real_t &a_mi, const real_t &a_mx); + void diag_ice_c_rng(const real_t &c_mi, const real_t &c_mx); void diag_kappa_rng(const real_t &r_mi, const real_t &r_mx); + void diag_ice(); + void diag_water(); void diag_dry_rng_cons(const real_t &r_mi, const real_t &r_mx); void diag_wet_rng_cons(const real_t &r_mi, const real_t &r_mx); + void diag_ice_a_rng_cons(const real_t &a_mi, const real_t &a_mx); + void diag_ice_c_rng_cons(const real_t &c_mi, const real_t &c_mx); void diag_kappa_rng_cons(const real_t &r_mi, const real_t &r_mx); + void diag_ice_cons(); + void diag_water_cons(); void diag_dry_mom(const int &k); void diag_wet_mom(const int &k); + void diag_ice_a_mom(const int &k); + void diag_ice_c_mom(const int &k); + void diag_ice_mix_ratio(); void diag_kappa_mom(const int &k); void diag_up_mom(const int&); void diag_vp_mom(const int&); @@ -196,6 +219,7 @@ namespace libcloudphxx void diag_RH_ge_Sc(); void diag_all(); void diag_precip_rate(); + void diag_precip_rate_ice_mass(); void diag_max_rw(); void diag_vel_div(); std::map diag_puddle(); @@ -275,12 +299,23 @@ namespace libcloudphxx void diag_RH(); void diag_dry_rng(const real_t &r_mi, const real_t &r_mx); void diag_wet_rng(const real_t &r_mi, const real_t &r_mx); + void diag_ice_a_rng(const real_t &a_mi, const real_t &a_mx); + void diag_ice_c_rng(const real_t &c_mi, const real_t &c_mx); void diag_kappa_rng(const real_t &r_mi, const real_t &r_mx); + void diag_ice(); + void diag_water(); void diag_dry_rng_cons(const real_t &r_mi, const real_t &r_mx); void diag_wet_rng_cons(const real_t &r_mi, const real_t &r_mx); + void diag_ice_a_rng_cons(const real_t &a_mi, const real_t &a_mx); + void diag_ice_c_rng_cons(const real_t &c_mi, const real_t &c_mx); void diag_kappa_rng_cons(const real_t &r_mi, const real_t &r_mx); + void diag_ice_cons(); + void diag_water_cons(); void diag_dry_mom(const int &k); void diag_wet_mom(const int &k); + void diag_ice_a_mom(const int &k); + void diag_ice_c_mom(const int &k); + void diag_ice_mix_ratio(); void diag_kappa_mom(const int&); void diag_up_mom(const int&); void diag_vp_mom(const int&); @@ -295,6 +330,7 @@ namespace libcloudphxx void diag_RH_ge_Sc(); void diag_all(); void diag_precip_rate(); + void diag_precip_rate_ice_mass(); void diag_max_rw(); void diag_vel_div(); std::map diag_puddle(); diff --git a/models/kinematic_2D/src/opts_lgrngn.hpp b/models/kinematic_2D/src/opts_lgrngn.hpp index ddae0d1aa..924921d24 100644 --- a/models/kinematic_2D/src/opts_lgrngn.hpp +++ b/models/kinematic_2D/src/opts_lgrngn.hpp @@ -297,7 +297,7 @@ void setopts_micro( */ rt_params.cloudph_opts_init.dry_distros.emplace( - setup.kappa, // key + libcloudphxx::lgrngn::kappa_rd_insol_t{setup.kappa, config::real_t(0.)}, // kappa, rd_insol std::make_shared> (setup) ); diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 4f7fcbfa6..ef57407d5 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -83,7 +83,12 @@ namespace libcloudphxx sstp_tmp_th, // ditto for theta sstp_tmp_rh, // ditto for rho sstp_tmp_p, // ditto for pressure - incloud_time; // time this SD has been within a cloud + incloud_time, // time this SD has been within a cloud + rd2_insol, // dry radii squared of insoluble aerosol + T_freeze, // freezing temperature + ice_a, // equatorial radius of ice + ice_c, // polar radius of ice + ice_rho; // ice apparent density // dry radii distribution characteristics real_t log_rd_min, // logarithm of the lower bound of the distr @@ -139,7 +144,8 @@ namespace libcloudphxx thrust_device::vector T, // temperature [K] p, // pressure [Pa] - RH, // relative humisity + RH, // relative humisity + RH_i, // relative humisity w.r.t. ice eta,// dynamic viscosity diss_rate; // turbulent kinetic energy dissipation rate @@ -196,14 +202,14 @@ namespace libcloudphxx > chem_stepper; // temporary data - tmp_vector_pool> + tmp_vector_pool> tmp_host_real_part, tmp_host_real_grid, tmp_host_real_cell; - tmp_vector_pool> + tmp_vector_pool> tmp_host_size_part, tmp_host_size_cell; - tmp_vector_pool> + tmp_vector_pool> tmp_device_real_part, tmp_device_real_cell; tmp_vector_pool> @@ -226,7 +232,8 @@ namespace libcloudphxx lambda_D_gp, lambda_K_gp, rw_mom3_gp, - rw3_gp; + rw3_gp, + ice_mass_gp; std::unique_ptr< typename tmp_vector_pool>::guard @@ -311,7 +318,9 @@ namespace libcloudphxx // max(1, n) int m1(int n) { return n == 0 ? 1 : n; } - // ctor + enum class phase_change { condensation, deposition }; // enum for choosing between phase change types + + // ctor impl(const opts_init_t &_opts_init, const std::pair &bcond, const int &mpi_rank, const int &mpi_size, const int &n_x_tot) : init_called(false), should_now_run_async(false), @@ -331,7 +340,7 @@ namespace libcloudphxx ), zero(0), n_part(0), - sorted(false), + sorted(false), n_user_params(_opts_init.kernel_parameters.size()), rng(_opts_init.rng_seed), src_stp_ctr(0), @@ -446,6 +455,16 @@ namespace libcloudphxx if(opts_init.diag_incloud_time) distmem_real_vctrs.insert({&incloud_time, detail::no_initial_value}); + + if(opts_init.ice_switch) + { + distmem_real_vctrs.insert({&rd2_insol, detail::no_initial_value}); + distmem_real_vctrs.insert({&ice_a, detail::no_initial_value}); + distmem_real_vctrs.insert({&ice_c, detail::no_initial_value}); + distmem_real_vctrs.insert({&ice_rho, detail::no_initial_value}); + if (! opts_init.time_dep_ice_nucl) + {distmem_real_vctrs.insert({&T_freeze, detail::no_initial_value});} + } // initializing distmem_n_vctrs - list of n_t vectors with properties of SDs that have to be copied/removed/recycled when a SD is copied/removed/recycled distmem_n_vctrs.insert(&n); @@ -463,6 +482,10 @@ namespace libcloudphxx if(opts_init.const_p) tmp_device_real_part.add_vector(); } + if (opts_init.ice_switch && opts_init.time_dep_ice_nucl) + { + tmp_device_real_part.add_vector(); + } // init number of temporary size_t vctrs // 1 needed if n_dims == 1 (i) @@ -486,7 +509,7 @@ namespace libcloudphxx void init_SD_with_distros_sd_conc(const common::unary_function &, const real_t &); void init_SD_with_distros_tail(const common::unary_function &, const real_t); void init_SD_with_distros_const_multi(const common::unary_function &); - void init_SD_with_distros_finalize(const real_t &, const bool unravel_ijk = true); + void init_SD_with_distros_finalize(const kappa_rd_insol_t &, const bool unravel_ijk = true); void init_SD_with_sizes(); void init_sanity_check( const arrinfo_t, const arrinfo_t, const arrinfo_t, @@ -519,6 +542,9 @@ namespace libcloudphxx void init_ijk(); void init_xyz(); void init_kappa(const real_t &); + void init_insol_dry_sizes(real_t); + void init_T_freeze(); + void init_a_c_rho_ice(); void init_incloud_time(); void init_count_num_sd_conc(const real_t & = 1); void init_count_num_const_multi(const common::unary_function &); @@ -572,6 +598,14 @@ namespace libcloudphxx void moms_ge0( const typename thrust_device::vector::iterator &vec_bgn ); + void moms_gt0( + const typename thrust_device::vector::iterator &vec_bgn, + const bool cons = false + ); + void moms_eq0( + const typename thrust_device::vector::iterator &vec_bgn, + const bool cons = false + ); void moms_rng( const real_t &min, const real_t &max, const typename thrust_device::vector::iterator &vec_bgn, @@ -582,15 +616,17 @@ namespace libcloudphxx const real_t &min, const real_t &max, const typename thrust_device::vector::iterator &vec_bgn, const bool cons - ); + ); + template // iterator type void moms_calc( - const typename thrust_device::vector::iterator &vec_bgn, + const it_t &vec_bgn, const thrust_size_t npart, const real_t power, const bool specific = true ); + template // iterator type void moms_calc( - const typename thrust_device::vector::iterator &vec_bgn, + const it_t &vec_bgn, const real_t power, const bool specific = true ); @@ -617,12 +653,16 @@ namespace libcloudphxx void sedi(const real_t &dt); void subs(const real_t &dt); + void ice_nucl_melt(const real_t &dt); + 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); 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_rv(thrust_device::vector &); + void update_th_rv(thrust_device::vector &, phase_change = phase_change::condensation); + void update_th_freezing(thrust_device::vector &); void update_state(thrust_device::vector &, thrust_device::vector &); void update_pstate(thrust_device::vector &, thrust_device::vector &); void update_incloud_time(const real_t &dt); @@ -636,15 +676,15 @@ namespace libcloudphxx void chem_react(const real_t &dt); void chem_cleanup(); void chem_post_step(); - + thrust_size_t rcyc(); void bcnd(); - void src(const real_t &dt); - void src_dry_distros_simple(const real_t &dt); - void src_dry_distros_matching(const real_t &dt); - void src_dry_distros(const real_t &dt); - void src_dry_sizes(const real_t &dt); + void src(const src_dry_distros_t &, const src_dry_sizes_t &); + void src_dry_distros_simple(const src_dry_distros_t &); + void src_dry_distros_matching(const src_dry_distros_t &); + void src_dry_distros(const src_dry_distros_t &); + void src_dry_sizes( const src_dry_sizes_t &); void rlx(const real_t); void rlx_dry_distros(const real_t); diff --git a/src/impl/particles_impl_bcnd.ipp b/src/impl/particles_impl_bcnd.ipp index 1148699c5..57556a3e5 100644 --- a/src/impl/particles_impl_bcnd.ipp +++ b/src/impl/particles_impl_bcnd.ipp @@ -23,17 +23,17 @@ namespace libcloudphxx template struct count_vol - { + { real_t exponent; count_vol(real_t exponent) : exponent(exponent){}; template BOOST_GPU_ENABLED - real_t operator()(const tuple &tup) + real_t operator()(const tuple &tup) // tup is a tuple (n, radius) { #if !defined(__NVCC__) using std::pow; #endif - return 4./3. + return 4./3. #if !defined(__NVCC__) * pi() #else @@ -44,6 +44,43 @@ namespace libcloudphxx thrust::get<1>(tup), // radius at some power exponent); } + }; + + template + struct count_ice_mass + { + count_ice_mass() {} + + template + BOOST_GPU_ENABLED + real_t operator()(const tuple &tup) // tup is a tuple (n, ice_a, ice_c, ice_rho) + { + return 4./3. + #if !defined(__NVCC__) + * pi() + #else + * CUDART_PI + #endif + * thrust::get<0>(tup) // n + * thrust::get<1>(tup) * thrust::get<1>(tup) // a^2 + * thrust::get<2>(tup) //c + * thrust::get<3>(tup); //rho_i + } + }; + + template + struct count_num + { + int ice; // 0 - water only, 1 - ice only, 2 - both + count_num(int ice) : ice(ice){assert(ice == 0 || ice == 1 || ice == 2);}; // TODO: ice enum + + template + BOOST_GPU_ENABLED + real_t operator()(const tuple &tup) // tup is a tuple (n_filtered, rw2) + { + if((ice==0 && thrust::get<1>(tup) == real_t(0)) || (ice==1 && thrust::get<1>(tup) > real_t(0))) return 0.; + return thrust::get<0>(tup); + } }; template @@ -223,7 +260,7 @@ namespace libcloudphxx output_puddle[common::outliq_vol] += thrust::transform_reduce( thrust::make_zip_iterator(thrust::make_tuple( - n_filtered.begin(), rw2.begin())), // input start + n_filtered.begin(), rw2.begin())), // input start thrust::make_zip_iterator(thrust::make_tuple( n_filtered.begin(), rw2.begin())) + n_part, // input end detail::count_vol(3./2.), // operation @@ -243,13 +280,53 @@ namespace libcloudphxx thrust::plus() ); + // add total number of water droplets that fell out in this step + output_puddle[common::outliq_num] += + thrust::transform_reduce( + thrust::make_zip_iterator(thrust::make_tuple( + n_filtered.begin(), rw2.begin())), // input start + thrust::make_zip_iterator(thrust::make_tuple( + n_filtered.begin(), rw2.begin())) + n_part, // input end + detail::count_num(0), // operation + real_t(0), // init val + thrust::plus() + ); + // add total number of particles that fell out in this step - output_puddle[common::outprtcl_num] += + output_puddle[common::outprtcl_num] += thrust::reduce( n_filtered.begin(), // input start n_filtered.begin() + n_part // input end ); + if (opts_init.ice_switch) + { + // add total ice mass that fell out in this step + output_puddle[common::outice_mass] += + thrust::transform_reduce( + thrust::make_zip_iterator(thrust::make_tuple( + n_filtered.begin(), ice_a.begin(), ice_c.begin(), ice_rho.begin())), // input start + thrust::make_zip_iterator(thrust::make_tuple( + n_filtered.begin(), ice_a.begin(), ice_c.begin(), ice_rho.begin())) + n_part, // input end + detail::count_ice_mass(), // operation + real_t(0), // init val + thrust::plus() + ); + + // add total number of ice droplets that fell out in this step + output_puddle[common::outice_num] += + thrust::transform_reduce( + thrust::make_zip_iterator(thrust::make_tuple( + n_filtered.begin(), rw2.begin())), // input start + thrust::make_zip_iterator(thrust::make_tuple( + n_filtered.begin(), rw2.begin())) + n_part, // input end + detail::count_num(1), // operation + real_t(0), // init val + thrust::plus() + ); + } + + if(opts_init.chem_switch) { for (int i = 0; i < chem_all; ++i) diff --git a/src/impl/particles_impl_cond.ipp b/src/impl/particles_impl_cond.ipp index d8c06960d..8df6b7e5b 100644 --- a/src/impl/particles_impl_cond.ipp +++ b/src/impl/particles_impl_cond.ipp @@ -19,35 +19,34 @@ namespace libcloudphxx namespace arg = thrust::placeholders; - thrust_device::vector &lambda_D(lambda_D_gp->get()); - thrust_device::vector &lambda_K(lambda_K_gp->get()); - - // --- calc liquid water content before cond --- - hskpng_sort(); - auto drv_g = tmp_device_real_cell.get_guard(); - thrust_device::vector &drv = drv_g.get(); + 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 + auto drv_liq_g = tmp_device_real_cell.get_guard(); + thrust_device::vector &drv_liq = drv_liq_g.get(); if(step == 0) reset_guardp(rw_mom3_gp, tmp_device_real_cell); - thrust_device::vector &rw_mom3 = rw_mom3_gp->get(); + thrust_device::vector &rw_mom3 = rw_mom3_gp->get(); - // calculating the 3rd wet moment before condensation + // Compute per-cell 3rd moment of liquid droplets before condensation. It is stored in count_mom if(step == 0) { - moms_all(); + moms_gt0(rw2.begin()); // choose liquid particles (rw2>0) moms_calc(rw2.begin(), real_t(3./2.)); nancheck_range(count_mom.begin(), count_mom.begin() + count_n, "count_mom (3rd wet moment) before condensation"); // fill with 0s if not all cells have particles - if(count_n!=n_cell) - { - thrust::fill(drv.begin(), drv.end(), real_t(0.)); + if(count_n!=n_cell) { + thrust::fill(drv_liq.begin(), drv_liq.end(), real_t(0.)); thrust::fill(rw_mom3.begin(), rw_mom3.end(), real_t(0.)); } - - // drv = -rw_mom3 precond + thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output + count_mom.begin(), count_mom.begin() + count_n, + thrust::make_permutation_iterator(drv_liq.begin(), count_ijk.begin()), thrust::negate() ); } @@ -55,8 +54,8 @@ namespace libcloudphxx { // drv = -rw_mom3 precond thrust::transform( - rw_mom3.begin(), rw_mom3.end(), - drv.begin(), + rw_mom3.begin(), rw_mom3.end(), + drv_liq.begin(), thrust::negate() ); } @@ -73,7 +72,7 @@ namespace libcloudphxx thrust::make_permutation_iterator(lambda_K.begin(), ijk.begin()) )); - // calculating drop growth in a timestep using backward Euler + // calculating drop growth in a timestep using backward Euler // TODO: both calls almost identical, use std::bind or sth? if(turb_cond) { @@ -86,6 +85,7 @@ namespace libcloudphxx arg::_1 + arg::_2 ); + // condensation for liquid droplets thrust::transform( rw2.begin(), rw2.end(), // input - 1st arg (zip not as 1st arg not to write zip.end() thrust::make_zip_iterator( // input - 2nd arg @@ -94,12 +94,14 @@ namespace libcloudphxx thrust::make_permutation_iterator(p.begin(), ijk.begin()), RH_plus_ssp.begin() ) - ), + ), rw2.begin(), // output detail::advance_rw2(dt, RH_max) ); } else + { + // condensation for liquid droplets thrust::transform( rw2.begin(), rw2.end(), // input - 1st arg (zip not as 1st arg not to write zip.end() thrust::make_zip_iterator( // input - 2nd arg @@ -108,45 +110,47 @@ namespace libcloudphxx thrust::make_permutation_iterator(p.begin(), ijk.begin()), thrust::make_permutation_iterator(RH.begin(), ijk.begin()) ) - ), + ), rw2.begin(), // output detail::advance_rw2(dt, RH_max) ); - nancheck(rw2, "rw2 after condensation (no sub-steps"); - - // calculating the 3rd wet moment after condensation - moms_calc(rw2.begin(), real_t(3./2.)); - nancheck_range(count_mom.begin(), count_mom.begin() + count_n, "count_mom (3rd wet moment) after condensation"); + nancheck(rw2, "rw2 after condensation (no sub-steps"); + } - if(step < sstp_cond - 1) - { - thrust::copy( - count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg - thrust::make_permutation_iterator(rw_mom3.begin(), count_ijk.begin()) // output - ); + // Compute per-cell 3rd moment of liquid droplets after condensation. It is stored in count_mom + moms_gt0(rw2.begin()); // choose liquid particles (rw2>0) + moms_calc(rw2.begin(), real_t(3./2.)); + nancheck_range(count_mom.begin(), count_mom.begin() + count_n, "count_mom (3rd wet moment) after condensation"); - // adding the third moment after condensation to dm_3 - thrust::transform( - rw_mom3.begin(), rw_mom3.end(), - drv.begin(), - drv.begin(), - thrust::plus() - ); - } - else // last step, calculate change in 3rd moment and update th and rv - { - thrust::transform( - count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // input - 2nd arg - thrust::make_permutation_iterator(drv.begin(), count_ijk.begin()), // output - thrust::plus() - ); - rw_mom3_gp.reset(); // destroy guard to tmp array that stored 3rd moment of rw - } + // Adding the third liquid moment after condensation to drv_liq + if(step < sstp_cond - 1) + { + thrust::copy( + count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg + thrust::make_permutation_iterator(rw_mom3.begin(), count_ijk.begin()) // output + ); - // update th and rv according to changes in third specific wet moment - update_th_rv(drv); + // adding the third moment after condensation to dm_3 + thrust::transform( + rw_mom3.begin(), rw_mom3.end(), + drv_liq.begin(), + drv_liq.begin(), + thrust::plus() + ); + } + else // last step, calculate change in 3rd moment and update th and rv + { + thrust::transform( + count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg + thrust::make_permutation_iterator(drv_liq.begin(), count_ijk.begin()), // input - 2nd arg + thrust::make_permutation_iterator(drv_liq.begin(), count_ijk.begin()), // output + thrust::plus() + ); + rw_mom3_gp.reset(); // destroy guard to tmp array that stored 3rd moment of rw + } - } - }; -}; + // update th and rv according to change in third specific wet moment + update_th_rv(drv_liq, impl::phase_change::condensation); + } + }; + } \ No newline at end of file diff --git a/src/impl/particles_impl_cond_common.ipp b/src/impl/particles_impl_cond_common.ipp index ca2a0f189..d24b4861f 100644 --- a/src/impl/particles_impl_cond_common.ipp +++ b/src/impl/particles_impl_cond_common.ipp @@ -126,13 +126,13 @@ namespace libcloudphxx const quantity K = K_0() * beta(lambda_K / rw) * (Nu(Pr, Re) / 2); - return real_t(2) * rdrdt( + return real_t(2) * rdrdt( D, K, - rhod * rv, - T, - p, - RH > RH_max ? RH_max : RH, + rhod * rv, + T, + p, + RH > RH_max ? RH_max : RH, a_w(rw3, rd3, kpa), klvntrm(rw, T) ); @@ -170,8 +170,11 @@ namespace libcloudphxx using std::isinf; #endif + // Skip ice particles + if (rw2_old <= 0) return rw2_old; + auto& tpl_in = thrust::get<0>(tpl); - const advance_rw2_minfun f(dt, rw2_old, tpl, RH_max); + const advance_rw2_minfun f(dt, rw2_old, tpl, RH_max); const real_t drw2 = dt * f.drw2_dt(rw2_old * si::square_metres) * si::seconds / si::square_metres; #if !defined(NDEBUG) @@ -193,7 +196,7 @@ namespace libcloudphxx "vt: %g " "lambda_D: %g " "lambda_K: %g\n", - drw2, rw2_old, dt, RH_max, + drw2, rw2_old, dt, RH_max, thrust::get<0>(tpl_in), // rhod thrust::get<1>(tpl_in), // rv thrust::get<2>(tpl_in), // T @@ -291,6 +294,7 @@ namespace libcloudphxx return rw2_new; } }; + }; }; }; diff --git a/src/impl/particles_impl_cond_sstp.ipp b/src/impl/particles_impl_cond_sstp.ipp index 62dc8f4d4..a0a9ff306 100644 --- a/src/impl/particles_impl_cond_sstp.ipp +++ b/src/impl/particles_impl_cond_sstp.ipp @@ -20,8 +20,8 @@ namespace libcloudphxx resolved_RH(RH_formula) {} - BOOST_GPU_ENABLED - real_t operator()(const thrust::tuple &tpl) + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple &tpl) { return resolved_RH(thrust::make_tuple(thrust::get<0>(tpl), thrust::get<1>(tpl), thrust::get<2>(tpl))) + thrust::get<3>(tpl); } @@ -36,9 +36,9 @@ namespace libcloudphxx const thrust_device::vector &Tp, const pres_iter &pi, const RH_iter &rhi - ) { - thrust_device::vector &lambda_D(lambda_D_gp->get()); - thrust_device::vector &lambda_K(lambda_K_gp->get()); + ) { + thrust_device::vector &lambda_D(lambda_D_gp->get()); + thrust_device::vector &lambda_K(lambda_K_gp->get()); auto hlpr_zip_iter = thrust::make_zip_iterator(thrust::make_tuple( sstp_tmp_rh.begin(), @@ -56,7 +56,7 @@ namespace libcloudphxx thrust::make_permutation_iterator(lambda_K.begin(), ijk.begin()) )); - // calculating drop growth in a timestep using backward Euler + // calculating drop growth in a timestep using backward Euler thrust::transform( rw2.begin(), rw2.end(), // input - 1st arg (zip not as 1st arg not to write zip.end() thrust::make_zip_iterator(thrust::make_tuple( // input - 2nd arg @@ -65,7 +65,7 @@ namespace libcloudphxx pi, // particle-specific RH rhi - )), + )), rw2.begin(), // output detail::advance_rw2(dt, RH_max) ); @@ -78,12 +78,12 @@ namespace libcloudphxx const real_t &RH_max, const bool turb_cond, const int step - ) { + ) { namespace arg = thrust::placeholders; // prerequisite - hskpng_sort(); + hskpng_sort(); // particle's local change in rv auto pdrv_g = tmp_device_real_part.get_guard(); thrust_device::vector &pdrv = pdrv_g.get(); @@ -121,8 +121,8 @@ namespace libcloudphxx sstp_tmp_th.begin(), sstp_tmp_th.end(), // input - first arg sstp_tmp_rh.begin(), // input - second arg Tp.begin(), // output - detail::common__theta_dry__T_rhod() - ); + detail::common__theta_dry__T_rhod() + ); } else // th_std { @@ -130,16 +130,16 @@ namespace libcloudphxx thrust::transform( sstp_tmp_th.begin(), sstp_tmp_th.end(), // input - first arg thrust::make_zip_iterator(thrust::make_tuple( - sstp_tmp_rv.begin(), // input - second arg + sstp_tmp_rv.begin(), // input - second arg sstp_tmp_p.begin() // input - third arg )), Tp.begin(), // output - detail::common__theta_std__T_p() + detail::common__theta_std__T_p() ); } - // calculating drop growth in a timestep using backward Euler + // calculating drop growth in a timestep using backward Euler if(!opts_init.const_p) { // particle-specific pressure iterator, used twice @@ -154,7 +154,7 @@ namespace libcloudphxx ); if(turb_cond) - cond_sstp_hlpr(dt, RH_max, Tp, + cond_sstp_hlpr(dt, RH_max, Tp, // particle-specific p pressure_iter, // particle-specific RH, resolved + SGS @@ -166,8 +166,8 @@ namespace libcloudphxx ssp.begin() )), detail::RH_sgs(opts_init.RH_formula) - ) - ); + ) + ); else // no RH SGS cond_sstp_hlpr(dt, RH_max, Tp, // particle-specific p @@ -180,8 +180,8 @@ namespace libcloudphxx Tp.begin() )), detail::RH(opts_init.RH_formula) - ) - ); + ) + ); } else // opts_init.const_p { @@ -198,8 +198,8 @@ namespace libcloudphxx ssp.begin() )), detail::RH_sgs(opts_init.RH_formula) - ) - ); + ) + ); else // no RH SGS cond_sstp_hlpr(dt, RH_max, Tp, // particle-specific p @@ -212,8 +212,8 @@ namespace libcloudphxx Tp.begin() )), detail::RH(opts_init.RH_formula) - ) - ); + ) + ); } // rw3_new - rw3_old @@ -258,23 +258,23 @@ namespace libcloudphxx - common::moist_air::rho_w() / si::kilograms * si::cubic_metres * real_t(4./3) * pi(), n_dims ) - ); + ); // apply change in rv to sstp_tmp_rv update_pstate(sstp_tmp_rv, pdrv); // calc particle-specific change in th based on pdrv thrust::transform( - thrust::make_zip_iterator(thrust::make_tuple( - pdrv.begin(), // + thrust::make_zip_iterator(thrust::make_tuple( + pdrv.begin(), // Tp.begin(), // dth = drv * d_th_d_rv(T, th) - sstp_tmp_th.begin() // - )), - thrust::make_zip_iterator(thrust::make_tuple( - pdrv.end(), // + sstp_tmp_th.begin() // + )), + thrust::make_zip_iterator(thrust::make_tuple( + pdrv.end(), // Tp.end(), // dth = drv * d_th_d_rv(T, th) - sstp_tmp_th.end() // - )), + sstp_tmp_th.end() // + )), pdrv.begin(), // in-place detail::dth() ); @@ -282,5 +282,5 @@ namespace libcloudphxx // apply change in th to sstp_tmp_th update_pstate(sstp_tmp_th, pdrv); } - }; -}; + }; +}; \ No newline at end of file diff --git a/src/impl/particles_impl_fill_outbuf.ipp b/src/impl/particles_impl_fill_outbuf.ipp index abf937cb3..82f8e3db9 100644 --- a/src/impl/particles_impl_fill_outbuf.ipp +++ b/src/impl/particles_impl_fill_outbuf.ipp @@ -39,16 +39,31 @@ namespace libcloudphxx template std::vector particles_t::impl::fill_attr_outbuf(const std::string &name) { - const std::set attr_names = {"rw2", "rd3", "kappa", "x", "y", "z"}; // TODO implement "n" - it is n_t type and others are real_t + const std::set attr_names = {"rw2", "rd3", "kappa", "rd2_insol", "T_freeze", "ice_a", "ice_c", "ice_rho", "x", "y", "z"}; // TODO implement "n" - it is n_t type and others are real_t if (std::find(std::begin(attr_names), std::end(attr_names), name) == std::end(attr_names)) throw std::runtime_error("Unknown attribute name passed to get_attr."); + if (!opts_init.ice_switch && + (name == "ice_a" || name == "ice_c" || name == "ice_rho" || name == "rd2_insol" || name == "T_freeze")) + { + throw std::runtime_error("Requested ice attribute '" + name + "' but ice_switch is off."); + } + if (opts_init.time_dep_ice_nucl && name == "T_freeze") + { + throw std::runtime_error("Requested T_freeze but singular ice nucleation is off."); + } + const thrust_device::vector &dv( name == "rw2" ? rw2 : name == "rd3" ? rd3 : - name == "kappa" ? kpa : - name == "x" ? x : - name == "y" ? y : + name == "kappa" ? kpa : + name == "rd2_insol" ? rd2_insol : + name == "T_freeze" ? T_freeze : + name == "ice_a" ? ice_a : + name == "ice_c" ? ice_c : + name == "ice_rho" ? ice_rho : + name == "x" ? x : + name == "y" ? y : z); // NOTE: for host backends (i.e. undefined __NVCC__) we could return the vector directly, without a copy; diff --git a/src/impl/particles_impl_hskpng_Tpr.ipp b/src/impl/particles_impl_hskpng_Tpr.ipp index 0ce9e70d8..3cbe9fe6e 100644 --- a/src/impl/particles_impl_hskpng_Tpr.ipp +++ b/src/impl/particles_impl_hskpng_Tpr.ipp @@ -102,6 +102,24 @@ namespace libcloudphxx / common::tetens::r_vs(T * si::kelvins, p * si::pascals)); } + template + BOOST_GPU_ENABLED + real_t RHi_pv_cc(const real_t &p, const real_t &rv, const real_t &T) + { + return real_t( + common::moist_air::p_v(p * si::pascals, quantity(rv)) + / common::const_cp::p_vsi(T * si::kelvins)); + } + + template + BOOST_GPU_ENABLED + real_t RHi_rv_cc(const real_t &p, const real_t &rv, const real_t &T) + { + return real_t( + rv + / common::const_cp::r_vsi(T * si::kelvins, p * si::pascals)); + } + template struct RH //: thrust::unary_function&, real_t> { @@ -158,7 +176,33 @@ namespace libcloudphxx } } }; - + + template + struct RH_i //: thrust::unary_function&, real_t> + { + const RH_formula_t RH_formula; + // the type of formula to be used for RH + RH_i(RH_formula_t RH_formula): + RH_formula(RH_formula) + {} + + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple &tpl) // p, rv, T + { + switch (RH_formula) + { + case RH_formula_t::pv_cc: + return RHi_pv_cc(thrust::get<0>(tpl), thrust::get<1>(tpl), thrust::get<2>(tpl)); + case RH_formula_t::rv_cc: + return RHi_rv_cc(thrust::get<0>(tpl), thrust::get<1>(tpl), thrust::get<2>(tpl)); + // NOTE: no Tetens formulas for ice + default: + assert(0); + return 0.; + } + } + }; + template struct common__vterm__visc //: thrust::unary_function// TODO: rename it! (vterm) visc_eta? { @@ -214,12 +258,28 @@ namespace libcloudphxx ); } + // RH 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.begin(), // output detail::RH(opts_init.RH_formula) ); + + 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 + ) + ); + } } // dynamic viscosity diff --git a/src/impl/particles_impl_hskpng_vterm.ipp b/src/impl/particles_impl_hskpng_vterm.ipp index 17c05ac29..02151c597 100644 --- a/src/impl/particles_impl_hskpng_vterm.ipp +++ b/src/impl/particles_impl_hskpng_vterm.ipp @@ -45,7 +45,7 @@ namespace libcloudphxx BOOST_GPU_ENABLED real_t operator()( const real_t &rw2, - const thrust::tuple &tpl + const thrust::tuple &tpl // (T, p, rhod, eta) ) { #if !defined(__NVCC__) using std::sqrt; @@ -106,7 +106,7 @@ namespace libcloudphxx BOOST_GPU_ENABLED real_t operator()( const real_t &rw2, - const thrust::tuple &tpl + const thrust::tuple &tpl // (vt_0, p, rhod, eta) ) { #if !defined(__NVCC__) using std::sqrt; @@ -126,13 +126,64 @@ namespace libcloudphxx return 0.; //sanity checks done in pimpl constructor } } - }; + }; + + template + struct common__vterm__ice + { + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple< + real_t, real_t, real_t, // ice_a, ice_c, rho_i + real_t, real_t, real_t, real_t // T, p, rhod, eta + >& tpl) const + { + const real_t &ice_a = thrust::get<0>(tpl); + const real_t &ice_c = thrust::get<1>(tpl); + const real_t &rho_i = thrust::get<2>(tpl); + const real_t &T = thrust::get<3>(tpl); + const real_t &p = thrust::get<4>(tpl); + const real_t &rhod = thrust::get<5>(tpl); + const real_t &eta = thrust::get<6>(tpl); + + return common::vterm::vt_beard76( + ice_a * si::metres, + T * si::kelvins, + p * si::pascals, + rhod * si::kilograms / si::cubic_metres, + eta * si::pascals * si::seconds + ) / si::metres_per_second + * real_t(common::moist_air::rho_i() / common::moist_air::rho_w()); + } + }; + + template + struct check_vt_invalid_rw2_zero + { + BOOST_GPU_ENABLED + bool operator()(const thrust::tuple& tpl) const + { + return thrust::get<0>(tpl) == real_t(invalid) && + thrust::get<1>(tpl) == real_t(0); + } + }; + + template + struct check_vt_invalid_rw2_nonzero + { + BOOST_GPU_ENABLED + bool operator()(const thrust::tuple& tpl) const + { + return thrust::get<0>(tpl) == real_t(invalid) && + thrust::get<1>(tpl) > real_t(0); + } + }; + }; template void particles_t::impl::hskpng_vterm_invalid() - { + { typedef thrust::permutation_iterator< typename thrust_device::vector::iterator, typename thrust_device::vector::iterator @@ -141,6 +192,7 @@ namespace libcloudphxx namespace arg = thrust::placeholders; + // liquid droplets if(opts_init.terminal_velocity == vt_t::beard77fast) //use cached vt at sea level { auto vt0_bin_g = tmp_device_real_part.get_guard(); // should be thrust_size_t, but we have many real_ available (do we?) @@ -148,10 +200,10 @@ namespace libcloudphxx // get cached bin number thrust::transform_if( rw2.begin(), rw2.end(), - vt.begin(), + thrust::make_zip_iterator(thrust::make_tuple(vt.begin(), rw2.begin())), vt0_bin.begin(), detail::get_vt0_bin(config.vt0_ln_r_min, config.vt0_ln_r_max, config.vt0_n_bin), - arg::_1 == real_t(detail::invalid) + detail::check_vt_invalid_rw2_nonzero() ); // calc the vt thrust::transform_if( @@ -161,11 +213,11 @@ namespace libcloudphxx thrust::make_permutation_iterator(p.begin(), ijk.begin()), thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin(), ijk.begin()) - )), // input - 2nd arg - vt.begin(), // condition argument + )), // input - 2nd arg + thrust::make_zip_iterator(thrust::make_tuple(vt.begin(), rw2.begin())), // stencil vt.begin(), // output detail::common__vterm__vt__cached(opts_init.terminal_velocity), - arg::_1 == real_t(detail::invalid) + detail::check_vt_invalid_rw2_nonzero() ); } // non-cached vt @@ -177,12 +229,36 @@ namespace libcloudphxx thrust::make_permutation_iterator(p.begin(), ijk.begin()), thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin(), ijk.begin()) - )), // input - 2nd arg - vt.begin(), // condition argument + )), // input - 2nd arg + thrust::make_zip_iterator(thrust::make_tuple(vt.begin(), rw2.begin())), // stencil vt.begin(), // output detail::common__vterm__vt(opts_init.terminal_velocity), - arg::_1 == real_t(detail::invalid) + detail::check_vt_invalid_rw2_nonzero() + ); + + // ice terminal velocity + if (opts_init.ice_switch) + { + auto zipped = thrust::make_zip_iterator(thrust::make_tuple( + ice_a.begin(), + ice_c.begin(), + ice_rho.begin(), + thrust::make_permutation_iterator(T.begin(), ijk.begin()), + thrust::make_permutation_iterator(p.begin(), ijk.begin()), + thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), + thrust::make_permutation_iterator(eta.begin(), ijk.begin()) + )); + auto zipped_end = zipped + ice_a.size(); + + thrust::transform_if( + zipped, + zipped_end, + thrust::make_zip_iterator(thrust::make_tuple(vt.begin(), rw2.begin())), + vt.begin(), + detail::common__vterm__ice(), + detail::check_vt_invalid_rw2_zero() ); + } } template @@ -193,7 +269,9 @@ namespace libcloudphxx typename thrust_device::vector::iterator > pi_t; typedef thrust::zip_iterator > zip_it_t; + namespace arg = thrust::placeholders; + // liquid droplets if(opts_init.terminal_velocity == vt_t::beard77fast) //use cached vt at sea level { // auto vt0_bin_g = tmp_device_size_part.get_guard(); @@ -201,26 +279,30 @@ namespace libcloudphxx auto vt0_bin_g = tmp_device_real_part.get_guard(); // should be thrust_size_t, but we have many real_ available (do we?) thrust_device::vector &vt0_bin(vt0_bin_g.get()); // get cached bin number - thrust::transform( + thrust::transform_if( rw2.begin(), rw2.end(), + rw2.begin(), vt0_bin.begin(), - detail::get_vt0_bin(config.vt0_ln_r_min, config.vt0_ln_r_max, config.vt0_n_bin) + detail::get_vt0_bin(config.vt0_ln_r_min, config.vt0_ln_r_max, config.vt0_n_bin), + arg::_1 > real_t(0) ); - // calc the vt - thrust::transform( - rw2.begin(), rw2.end(), // input - 1st arg + // calc the vt for liquid droplets + thrust::transform_if( + rw2.begin(), rw2.end(), // input - 1st arg thrust::make_zip_iterator(thrust::make_tuple( thrust::make_permutation_iterator(vt_0.begin(), vt0_bin.begin()), thrust::make_permutation_iterator(p.begin(), ijk.begin()), thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin(), ijk.begin()) - )), // input - 2nd arg + )), // input - 2nd arg + rw2.begin(), // stencil vt.begin(), // output - detail::common__vterm__vt__cached(opts_init.terminal_velocity) + detail::common__vterm__vt__cached(opts_init.terminal_velocity), + arg::_1 > real_t(0) ); } else - thrust::transform( + thrust::transform_if( rw2.begin(), rw2.end(), // input - 1st arg zip_it_t(thrust::make_tuple( thrust::make_permutation_iterator(T.begin(), ijk.begin()), @@ -228,9 +310,35 @@ namespace libcloudphxx thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), thrust::make_permutation_iterator(eta.begin(), ijk.begin()) )), // input - 2nd arg + rw2.begin(), // stencil vt.begin(), // output - detail::common__vterm__vt(opts_init.terminal_velocity) + detail::common__vterm__vt(opts_init.terminal_velocity), + arg::_1 > real_t(0) ); + + // ice terminal velocity + if (opts_init.ice_switch) + { + auto zipped = thrust::make_zip_iterator(thrust::make_tuple( + ice_a.begin(), + ice_c.begin(), + ice_rho.begin(), + thrust::make_permutation_iterator(T.begin(), ijk.begin()), + thrust::make_permutation_iterator(p.begin(), ijk.begin()), + thrust::make_permutation_iterator(rhod.begin(), ijk.begin()), + thrust::make_permutation_iterator(eta.begin(), ijk.begin()) + )); + auto zipped_end = zipped + ice_a.size(); + + thrust::transform_if( + zipped, + zipped_end, + rw2.begin(), // stencil + vt.begin(), // output + detail::common__vterm__ice(), + arg::_1 == real_t(0) + ); + } } }; }; diff --git a/src/impl/particles_impl_ice_dep.ipp b/src/impl/particles_impl_ice_dep.ipp new file mode 100644 index 000000000..9e987b3b1 --- /dev/null +++ b/src/impl/particles_impl_ice_dep.ipp @@ -0,0 +1,298 @@ +// 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 + +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) + 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); + } + }; + + } + + + 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 + auto drv_ice_g = tmp_device_real_cell.get_guard(); + thrust_device::vector &drv_ice = drv_ice_g.get(); + if(step == 0) + reset_guardp(ice_mass_gp, tmp_device_real_cell); + thrust_device::vector &ice_mass = ice_mass_gp->get(); + + // Compute per-cell 3rd moment of ice before deposition. It is stored in count_mom + if(step == 0) + { + 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) before deposition"); + + // fill with 0s if not all cells have particles + if(count_n!=n_cell) { + thrust::fill(drv_ice.begin(), drv_ice.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(drv_ice.begin(), count_ijk.begin()), + thrust::negate() + ); + } + else // copy ice_mass from previous step + { + // drv = -ice_mass precond + thrust::transform( + ice_mass.begin(), ice_mass.end(), + drv_ice.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, 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 drv_ice + if(step < sstp_cond - 1) + { + 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(), + drv_ice.begin(), + drv_ice.begin(), + thrust::plus() + ); + } + else // last step, calculate change in 3rd moment and update th and rv + { + thrust::transform( + count_mom.begin(), count_mom.begin() + count_n, // input - 1st arg + thrust::make_permutation_iterator(drv_ice.begin(), count_ijk.begin()), // input - 2nd arg + thrust::make_permutation_iterator(drv_ice.begin(), count_ijk.begin()), // output + thrust::plus() + ); + ice_mass_gp.reset(); // destroy guard to tmp array that stored ice_mass + } + + // update th and rv according to change in third specific moment + update_th_rv(drv_ice, impl::phase_change::deposition); + } + }; + } \ No newline at end of file diff --git a/src/impl/particles_impl_ice_nucl_melt.ipp b/src/impl/particles_impl_ice_nucl_melt.ipp new file mode 100644 index 000000000..811ee6ac7 --- /dev/null +++ b/src/impl/particles_impl_ice_nucl_melt.ipp @@ -0,0 +1,230 @@ +// vim:filetype=cpp +/** @file + * @copyright University of Warsaw + * @section LICENSE + * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) + */ + +#include +#include + +namespace libcloudphxx +{ + namespace lgrngn + { + namespace detail + { + // Singular immersion freezing (Shima et al., 2020) + // Functor to update rw2, a, c, rho_i, of frozen droplets + template + class singular_freeze + { + public: + BOOST_GPU_ENABLED + void operator()(thrust::tuple< + real_t&, real_t&, real_t&, real_t&, // to be updated (rw2, a, c, rho_i) + const real_t&, const real_t&, const real_t& // T_freeze, T, RH + > tpl) const + { + auto& rw2 = thrust::get<0>(tpl); + auto& a = thrust::get<1>(tpl); + auto& c = thrust::get<2>(tpl); + auto& rho_i = thrust::get<3>(tpl); + + const real_t T_freeze = thrust::get<4>(tpl); + const real_t T = thrust::get<5>(tpl); + const real_t RH = thrust::get<6>(tpl); + + if (rw2 > real_t(0) && T_freeze >= T && RH >= real_t(1)) // condition for freezing + { + rho_i = common::moist_air::rho_i() * si::cubic_metres / si::kilograms; + a = pow(rw2, real_t(0.5)) * pow(common::moist_air::rho_w() / common::moist_air::rho_i(), real_t(1./3.)); + c = pow(rw2, real_t(0.5)) * pow(common::moist_air::rho_w() / common::moist_air::rho_i(), real_t(1./3.)); + rw2 = real_t(0); + } + } + }; + + // Time-dependent immersion freezing (Arabas et al., 2025) + // Functor to update rw2, a, c, rho_i, of frozen droplets + template + class time_dep_freeze + { + const real_t dt; + const common::ice_nucleation::INP_t INP_type; + + public: + BOOST_GPU_ENABLED + time_dep_freeze(const real_t &dt, const common::ice_nucleation::INP_t INP_type) + : dt(dt), INP_type(INP_type) + {} + + BOOST_GPU_ENABLED + void operator()(thrust::tuple< + real_t&, real_t&, real_t&, real_t&, // to be updated (rw2, a, c, rho_i) + const real_t&, const real_t&, const real_t& // rd2_insol, u01, T + > tpl) const + { + auto& rw2 = thrust::get<0>(tpl); + auto& a = thrust::get<1>(tpl); + auto& c = thrust::get<2>(tpl); + auto& rho_i = thrust::get<3>(tpl); + + const real_t rd2_insol = thrust::get<4>(tpl); + const real_t u01 = thrust::get<5>(tpl); + const real_t T = thrust::get<6>(tpl); + + if (rw2 > real_t(0) && u01 < common::ice_nucleation::p_freeze(INP_type, rd2_insol, rw2, T, dt)) + { + rho_i = common::moist_air::rho_i() * si::cubic_metres / si::kilograms; + a = pow(rw2, real_t(0.5)) * pow(common::moist_air::rho_w() / common::moist_air::rho_i(), real_t(1./3.)); + c = pow(rw2, real_t(0.5)) * pow(common::moist_air::rho_w() / common::moist_air::rho_i(), real_t(1./3.)); + rw2 = real_t(0); + } + } + }; + + // Functor to update rw2, a, c, rho_i, of melted ice + template + class melt + { + public: + BOOST_GPU_ENABLED + void operator()(thrust::tuple< + real_t&, real_t&, real_t&, real_t&, // to be updated (rw2, a, c, rho_i) + const real_t& // ambient T + > tpl) const + { + auto& rw2 = thrust::get<0>(tpl); + auto& a = thrust::get<1>(tpl); + auto& c = thrust::get<2>(tpl); + auto& rho_i = thrust::get<3>(tpl); + + if (a*c > real_t(0) && thrust::get<4>(tpl) > real_t(273.15)) // if T > 0 C + { + rw2 = pow(common::moist_air::rho_i() / common::moist_air::rho_w() * c , real_t(2./3.)) * pow(a , real_t(4./3.)); + rho_i = real_t(0); + a = real_t(0); + c = real_t(0); + } + } + }; + + } + + // Immersion freezing and melting + template + void particles_t::impl::ice_nucl_melt(const real_t &dt) { + + hskpng_sort(); + + // A vector to store liquid 3rd moments + auto drw_g = tmp_device_real_cell.get_guard(); + thrust_device::vector &drw = drw_g.get(); + + // Compute per-cell 3rd moment of liquid droplets (sum of n*r^3) before freezing/melting. It is stored in count_mom + moms_eq0(ice_a.begin()); // choose liquid particles (ice_a=0) + moms_calc(rw2.begin(), real_t(1.5)); + nancheck_range(count_mom.begin(), count_mom.begin() + count_n, "count_mom (3rd wet moment) of droplets before freezing/melting"); + if (count_n != n_cell) + thrust::fill(drw.begin(), drw.end(), real_t(0.)); + + thrust::transform( + count_mom.begin(), count_mom.begin() + count_n, + thrust::make_permutation_iterator(drw.begin(), count_ijk.begin()), + thrust::negate() + ); + + // Change liquid droplets to ice under the freezing condition + + if (opts_init.time_dep_ice_nucl) // time dependent freezing based on Arabas et al., 2025 + { + auto u01g = tmp_device_real_part.get_guard(); + thrust_device::vector &u01 = u01g.get(); + rand_u01(u01, n_part); // random numbers between [0,1] for each particle + thrust::for_each( + thrust::make_zip_iterator(thrust::make_tuple( + rw2.begin(), + ice_a.begin(), + ice_c.begin(), + ice_rho.begin(), + rd2_insol.begin(), + u01.begin(), + thrust::make_permutation_iterator(T.begin(), ijk.begin()) + )), + thrust::make_zip_iterator(thrust::make_tuple( + rw2.begin(), + ice_a.begin(), + ice_c.begin(), + ice_rho.begin(), + rd2_insol.begin(), + u01.begin(), + thrust::make_permutation_iterator(T.begin(), ijk.begin()) + )) + n_part, + detail::time_dep_freeze(dt, opts_init.inp_type) // functor for updating (rw2, a, c, rho_i) if freezing condition satisfied + ); + } + else // singular freezing based on Shima et al., 2020 + { + thrust::for_each( + thrust::make_zip_iterator(thrust::make_tuple( + rw2.begin(), + ice_a.begin(), + ice_c.begin(), + ice_rho.begin(), + T_freeze.begin(), + thrust::make_permutation_iterator(T.begin(), ijk.begin()), + thrust::make_permutation_iterator(RH.begin(), ijk.begin()) + )), + thrust::make_zip_iterator(thrust::make_tuple( + rw2.begin(), + ice_a.begin(), + ice_c.begin(), + ice_rho.begin(), + T_freeze.begin(), + thrust::make_permutation_iterator(T.begin(), ijk.begin()), + thrust::make_permutation_iterator(RH.begin(), ijk.begin()) + )) + n_part, + detail::singular_freeze() // functor for updating (rw2, a, c, rho_i) if freezing condition satisfied + ); + } + + // Change ice to liquid droplets under the melting condition + thrust::for_each( + thrust::make_zip_iterator(thrust::make_tuple( + rw2.begin(), + ice_a.begin(), + ice_c.begin(), + ice_rho.begin(), + thrust::make_permutation_iterator(T.begin(), ijk.begin()) + )), + thrust::make_zip_iterator(thrust::make_tuple( + rw2.begin(), + ice_a.begin(), + ice_c.begin(), + ice_rho.begin(), + thrust::make_permutation_iterator(T.begin(), ijk.begin()) + )) + n_part, + detail::melt() // functor for updating (rw2, a, c, rho_i) if melting condition satisfied + ); + + // Compute per-cell 3rd moment of liquid droplets after freezing/melting. It is stored in count_mom + moms_eq0(ice_a.begin()); // choose liquid particles (ice_a=0) + moms_calc(rw2.begin(), real_t(1.5)); + nancheck_range(count_mom.begin(), count_mom.begin() + count_n, "count_mom (3rd wet moment) of droplets after freezing/melting"); + + // Compute the difference between liquid moments before and after + thrust::transform( + count_mom.begin(), count_mom.begin() + count_n, + thrust::make_permutation_iterator(drw.begin(), count_ijk.begin()), + thrust::make_permutation_iterator(drw.begin(), count_ijk.begin()), + thrust::plus() + ); + + // Update th according to the frozen volume per cell + update_th_freezing(drw); + + } + + } +} diff --git a/src/impl/particles_impl_init_SD_with_distros.ipp b/src/impl/particles_impl_init_SD_with_distros.ipp index 0a20ba6c1..2cea4a1b3 100644 --- a/src/impl/particles_impl_init_SD_with_distros.ipp +++ b/src/impl/particles_impl_init_SD_with_distros.ipp @@ -17,7 +17,7 @@ namespace libcloudphxx // calc sum of ln(rd) ranges of all distributions real_t tot_lnrd_rng = 0.; if(opts_init.sd_conc > 0) - for (typename opts_init_t::dry_distros_t::const_iterator ddi = opts_init.dry_distros.begin(); ddi != opts_init.dry_distros.end(); ++ddi) + for (auto ddi = opts_init.dry_distros.cbegin(); ddi != opts_init.dry_distros.cend(); ++ddi) { dist_analysis_sd_conc( *(ddi->second), @@ -27,7 +27,7 @@ namespace libcloudphxx } // initialize SDs of each kappa-type - for (typename opts_init_t::dry_distros_t::const_iterator ddi = opts_init.dry_distros.begin(); ddi != opts_init.dry_distros.end(); ++ddi) + for (auto ddi = opts_init.dry_distros.cbegin(); ddi != opts_init.dry_distros.cend(); ++ddi) { if(opts_init.sd_conc > 0) { @@ -50,10 +50,20 @@ namespace libcloudphxx // final inits common for tail/sd_conc/const_multi template - void particles_t::impl::init_SD_with_distros_finalize(const real_t &kappa, const bool unravel_ijk_switch) + void particles_t::impl::init_SD_with_distros_finalize(const kappa_rd_insol_t &kpa_rd_insol, const bool unravel_ijk_switch) { // init kappa - init_kappa(kappa); + init_kappa(kpa_rd_insol.kappa); + + if (opts_init.ice_switch) + { + init_insol_dry_sizes(kpa_rd_insol.rd_insol); + init_a_c_rho_ice(); + if (! opts_init.time_dep_ice_nucl) + { + init_T_freeze(); + } + } // initialising wet radii init_wet(); diff --git a/src/impl/particles_impl_init_SD_with_sizes.ipp b/src/impl/particles_impl_init_SD_with_sizes.ipp index d09dee3cd..b462d4bdb 100644 --- a/src/impl/particles_impl_init_SD_with_sizes.ipp +++ b/src/impl/particles_impl_init_SD_with_sizes.ipp @@ -14,20 +14,20 @@ namespace libcloudphxx template void particles_t::impl::init_SD_with_sizes() { - using dry_sizes_t = typename opts_init_t::dry_sizes_t; - using kappa_t = typename dry_sizes_t::key_type; - using size_number_t = typename dry_sizes_t::mapped_type; +// using dry_sizes_t = typename opts_init_t::dry_sizes_t; + // using size_number_t = typename dry_sizes_t::mapped_type; //using conc_multi_t = typename size_number_t::mapped_type; - // loop over kappas - for (typename dry_sizes_t::const_iterator dsi = opts_init.dry_sizes.begin(); dsi != opts_init.dry_sizes.end(); ++dsi) + // loop over (kappa, rd_insol) pairs + for (auto dsi = opts_init.dry_sizes.cbegin(); dsi != opts_init.dry_sizes.cend(); ++dsi) { - const kappa_t &kappa(dsi->first); - const size_number_t &size_number_map(dsi->second); + const real_t &kappa(dsi->first.kappa); + const real_t &rd_insol(dsi->first.rd_insol); + const auto &size_number_map(dsi->second); - // loop over the "size : {concentration, count}" pairs for this kappa - for (typename size_number_t::const_iterator sni = size_number_map.begin(); sni != size_number_map.end(); ++sni) + // loop over the "size : {concentration, count}" pairs for this (kappa, rd_insol) pair + for (auto sni = size_number_map.cbegin(); sni != size_number_map.cend(); ++sni) { // init number of SDs of this kappa in cells init_count_num_dry_sizes(sni->second); @@ -45,8 +45,18 @@ namespace libcloudphxx // initialising dry radii (needs ijk) init_dry_dry_sizes(sni->first); - // init kappa + // init kappa and rd_insol init_kappa(kappa); + + if (opts_init.ice_switch) + { + init_insol_dry_sizes(rd_insol); + init_a_c_rho_ice(); + if (! opts_init.time_dep_ice_nucl) + { + init_T_freeze(); + } + } // init multiplicities init_n_dry_sizes(sni->second.first, sni->second.second); diff --git a/src/impl/particles_impl_init_T_freeze.ipp b/src/impl/particles_impl_init_T_freeze.ipp new file mode 100644 index 000000000..ee0f5da55 --- /dev/null +++ b/src/impl/particles_impl_init_T_freeze.ipp @@ -0,0 +1,34 @@ +// vim:filetype=cpp +/** @file + * @copyright University of Warsaw + * @section LICENSE + * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) + * @brief initialisation routine for super droplets + */ + +#include + +namespace libcloudphxx +{ + namespace lgrngn + { + template + void particles_t::impl::init_T_freeze() { + + using namespace libcloudphxx::common::ice_nucleation; + + // random numbers between [0,1] for sampling + auto u01g = tmp_device_real_part.get_guard(); + thrust_device::vector &u01 = u01g.get(); + 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())), + T_freeze.begin() + n_part_old, + T_freeze_CDF_inv_functor(opts_init.inp_type) + ); + } + }; +}; + diff --git a/src/impl/particles_impl_init_a_c_rho_ice.ipp b/src/impl/particles_impl_init_a_c_rho_ice.ipp new file mode 100644 index 000000000..bf9e4949d --- /dev/null +++ b/src/impl/particles_impl_init_a_c_rho_ice.ipp @@ -0,0 +1,22 @@ +// vim:filetype=cpp +/** @file + * @copyright University of Warsaw + * @section LICENSE + * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) + * @brief initialisation routine for super droplets + */ + +namespace libcloudphxx +{ + namespace lgrngn + { + template + void particles_t::impl::init_a_c_rho_ice() + { + // filling a and c with zeros initially + thrust::fill(ice_a.begin() + n_part_old, ice_a.end(), real_t(0)); + thrust::fill(ice_c.begin() + n_part_old, ice_c.end(), real_t(0)); + thrust::fill(ice_rho.begin() + n_part_old, ice_rho.end(), real_t(0)); + } + }; +}; \ No newline at end of file diff --git a/src/impl/particles_impl_init_hskpng_ncell.ipp b/src/impl/particles_impl_init_hskpng_ncell.ipp index 45cc9ca6a..951ae8054 100644 --- a/src/impl/particles_impl_init_hskpng_ncell.ipp +++ b/src/impl/particles_impl_init_hskpng_ncell.ipp @@ -14,8 +14,10 @@ namespace libcloudphxx { // memory allocation, p already resized in init_sync() T.resize(n_cell); - RH.resize(n_cell); - eta.resize(n_cell); + RH.resize(n_cell); + eta.resize(n_cell); + if (opts_init.ice_switch) + RH_i.resize(n_cell); count_ijk.resize(n_cell); count_num.resize(n_cell); diff --git a/src/impl/particles_impl_init_insol_dry_sizes.ipp b/src/impl/particles_impl_init_insol_dry_sizes.ipp new file mode 100644 index 000000000..3facc1a97 --- /dev/null +++ b/src/impl/particles_impl_init_insol_dry_sizes.ipp @@ -0,0 +1,23 @@ +// vim:filetype=cpp +/** @file + * @copyright University of Warsaw + * @section LICENSE + * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) + * @brief initialisation routine for super droplets + */ + + +namespace libcloudphxx +{ + namespace lgrngn + { + template + void particles_t::impl::init_insol_dry_sizes( + real_t radius + ) + { + real_t rad2 = radius * radius; + thrust::fill(rd2_insol.begin() + n_part_old, rd2_insol.end(), rad2); + } + }; +}; diff --git a/src/impl/particles_impl_init_sanity_check.ipp b/src/impl/particles_impl_init_sanity_check.ipp index 2112b4df3..a27edf1db 100644 --- a/src/impl/particles_impl_init_sanity_check.ipp +++ b/src/impl/particles_impl_init_sanity_check.ipp @@ -57,21 +57,12 @@ namespace libcloudphxx if (opts_init.chem_switch && opts_init.src_type!=src_t::off) throw std::runtime_error("libcloudph++: chemistry and aerosol source are not compatible"); - if (opts_init.src_type!=src_t::off && opts_init.src_dry_distros.empty() && opts_init.src_dry_sizes.empty()) - throw std::runtime_error("libcloudph++: CCN source enabled, but src_dry_distros and src_dry_sizes are empty"); - - if (opts_init.src_type!=src_t::off && opts_init.src_dry_distros.size() > 1) - throw std::runtime_error("libcloudph++: src_dry_distros can only have a single kappa value."); - if (opts_init.src_type==src_t::matching && opts_init.dry_distros.size() > 1) throw std::runtime_error("libcloudph++: For 'matching' CCN source, the initial aerosol distribution can only have one kappa value (na kappa matching done)."); if (opts_init.src_type!=src_t::off && n_dims<2) throw std::runtime_error("libcloudph++: CCN source works in 2D and 3D only."); - if (opts_init.src_type==src_t::matching && !opts_init.src_dry_distros.empty() && - opts_init.src_dry_distros.begin()->first != opts_init.dry_distros.begin()->first) throw std::runtime_error("libcloudph++: For 'matching' CCN source, kappa of the source has to be the same as that of the initial profile (no kappa matching done)"); - if(opts_init.dry_distros.size() > 1 && opts_init.chem_switch) throw std::runtime_error("libcloudph++: chemistry and multiple kappa distributions are not compatible"); @@ -155,6 +146,18 @@ namespace libcloudphxx throw std::runtime_error("libcloudph++: In const_p option, pressure profile must be passed (p in init())"); if(!opts_init.const_p && !p.is_null()) throw std::runtime_error("libcloudph++: pressure profile was passed in init(), but the constant pressure option was not used"); + + if(opts_init.ice_switch) + { + if(opts_init.coal_switch) // because we dont know what to do when ice collides with water + throw std::runtime_error("libcloudph++: coalescence does not work with ice (turn off ice_switch or coal_switch)."); + if(opts_init.rlx_switch) // because we dont account for ice/water when matching and initializing aerosols from relaxation + 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 + throw std::runtime_error("libcloudph++: SGS condensation does not work with ice."); + } } }; }; diff --git a/src/impl/particles_impl_moms.ipp b/src/impl/particles_impl_moms.ipp index b85f4b817..502515447 100644 --- a/src/impl/particles_impl_moms.ipp +++ b/src/impl/particles_impl_moms.ipp @@ -151,7 +151,85 @@ namespace libcloudphxx arg::_1 * (arg::_2 >= 0) // op ); } + selected_before_counting = true; + } + + + // selects particles for which vec[i] > 0 + template + void particles_t::impl::moms_gt0( + const typename thrust_device::vector::iterator &vec_bgn, + const bool cons + ) + { + hskpng_sort(); + + if(!cons) + reset_guardp(n_filtered_gp, tmp_device_real_part); + thrust_device::vector &n_filtered = n_filtered_gp->get(); + + { + namespace arg = thrust::placeholders; + + if(!cons) + { + thrust::transform( + n.begin(), n.end(), // input - 1st arg + vec_bgn, // input - 2nd arg + n_filtered.begin(), // output + arg::_1 * (arg::_2 > 0) // op + ); + } + else + { + thrust::transform( + n_filtered.begin(), n_filtered.end(), // input - 1st arg + vec_bgn, // input - 2nd arg + n_filtered.begin(), // output + arg::_1 * (arg::_2 > 0) // op + ); + } + } + selected_before_counting = true; + } + + // selects particles for which vec[i] = 0 + template + void particles_t::impl::moms_eq0( + const typename thrust_device::vector::iterator &vec_bgn, + const bool cons + ) + { + hskpng_sort(); + + if(!cons) + reset_guardp(n_filtered_gp, tmp_device_real_part); + + thrust_device::vector &n_filtered = n_filtered_gp->get(); + + { + namespace arg = thrust::placeholders; + + if(!cons) + { + thrust::transform( + n.begin(), n.end(), // input - 1st arg + vec_bgn, // input - 2nd arg + n_filtered.begin(), // output + arg::_1 * (arg::_2 == 0) // op + ); + } + else + { + thrust::transform( + n_filtered.begin(), n_filtered.end(), // input - 1st arg + vec_bgn, // input - 2nd arg + n_filtered.begin(), // output + arg::_1 * (arg::_2 == 0) // op + ); + } + } selected_before_counting = true; } @@ -195,8 +273,9 @@ namespace libcloudphxx }; template + template // iterator type void particles_t::impl::moms_calc( - const typename thrust_device::vector::iterator &vec_bgn, + const it_t &vec_bgn, const thrust_size_t npart, const real_t power, const bool specific @@ -206,12 +285,6 @@ namespace libcloudphxx thrust_device::vector &n_filtered = n_filtered_gp->get(); - typedef thrust::permutation_iterator< - typename thrust_device::vector::const_iterator, - typename thrust_device::vector::iterator - > pi_t; - typedef thrust::zip_iterator > zip_it_t; - thrust::pair< thrust_device::vector::iterator, typename thrust_device::vector::iterator @@ -220,9 +293,9 @@ namespace libcloudphxx sorted_ijk.begin(), sorted_ijk.begin()+npart, // input - values thrust::make_transform_iterator( - zip_it_t(thrust::make_tuple( - pi_t(n_filtered.begin(), sorted_id.begin()), - pi_t(vec_bgn, sorted_id.begin()) + thrust::make_zip_iterator(thrust::make_tuple( + thrust::make_permutation_iterator(n_filtered.begin(), sorted_id.begin()), + thrust::make_permutation_iterator(vec_bgn, sorted_id.begin()) )), detail::moment_counter(power) ), @@ -314,8 +387,9 @@ namespace libcloudphxx } template + template // iterator type void particles_t::impl::moms_calc( - const typename thrust_device::vector::iterator &vec_bgn, + const it_t &vec_bgn, const real_t power, const bool specific ) diff --git a/src/impl/particles_impl_reserve_hskpng_npart.ipp b/src/impl/particles_impl_reserve_hskpng_npart.ipp index 2dc024542..a638f08d6 100644 --- a/src/impl/particles_impl_reserve_hskpng_npart.ipp +++ b/src/impl/particles_impl_reserve_hskpng_npart.ipp @@ -33,6 +33,18 @@ namespace libcloudphxx dot_ssp.reserve(opts_init.n_sd_max); } + if(opts_init.ice_switch) + { + rd2_insol.reserve(opts_init.n_sd_max); + ice_a.reserve(opts_init.n_sd_max); + ice_c.reserve(opts_init.n_sd_max); + ice_rho.reserve(opts_init.n_sd_max); + if (! opts_init.time_dep_ice_nucl) + { + T_freeze.reserve(opts_init.n_sd_max); + } + } + vt.reserve(opts_init.n_sd_max); sorted_id.reserve(opts_init.n_sd_max); diff --git a/src/impl/particles_impl_rlx_dry_distros.ipp b/src/impl/particles_impl_rlx_dry_distros.ipp index a8107c1a4..00e270427 100644 --- a/src/impl/particles_impl_rlx_dry_distros.ipp +++ b/src/impl/particles_impl_rlx_dry_distros.ipp @@ -302,7 +302,8 @@ namespace libcloudphxx n_part_to_init = n_part - n_part_old; // hskpng_resize_npart(); - init_SD_with_distros_finalize(kappa, false); // no need to unravel ijk there, becaues i j k are already initialized + init_SD_with_distros_finalize(lgrngn::kappa_rd_insol_t{kappa, real_t(0)}, false); // no need to unravel ijk there, because i j k are already initialized + // TODO: we assume that relaxation produces water (not ice) // TODO: asserts of newly added SD parameters? e.g. how many SD, how big is multiplicity etc. } // end of the distros loop diff --git a/src/impl/particles_impl_src.ipp b/src/impl/particles_impl_src.ipp index 4044b5cab..35512bd8e 100644 --- a/src/impl/particles_impl_src.ipp +++ b/src/impl/particles_impl_src.ipp @@ -10,15 +10,15 @@ namespace libcloudphxx namespace lgrngn { template - void particles_t::impl::src(const real_t &dt) + void particles_t::impl::src(const src_dry_distros_t &sdd, const src_dry_sizes_t &sds) { // ante_adding_SD(); - if(!opts_init.src_dry_distros.empty()) - src_dry_distros(dt); + if(!sdd.empty()) + src_dry_distros(sdd); - if(!opts_init.src_dry_sizes.empty()) - src_dry_sizes(dt); + if(!sds.empty()) + src_dry_sizes(sds); // post_adding_SD(); } diff --git a/src/impl/particles_impl_src_dry_distros.ipp b/src/impl/particles_impl_src_dry_distros.ipp index ffd4ea560..d3c57a40d 100644 --- a/src/impl/particles_impl_src_dry_distros.ipp +++ b/src/impl/particles_impl_src_dry_distros.ipp @@ -14,12 +14,18 @@ namespace libcloudphxx { // create new aerosol particles based on a size distribution template - void particles_t::impl::src_dry_distros(const real_t &dt) + void particles_t::impl::src_dry_distros(const src_dry_distros_t &sdd) { + if (sdd.size() > 1) + throw std::runtime_error("libcloudph++: src_dry_distros can only have a single kappa value."); + + if (opts_init.src_type == src_t::matching && !sdd.empty() && sdd.begin()->first.kappa != opts_init.dry_distros.begin()->first.kappa) + throw std::runtime_error("libcloudph++: For 'matching' CCN source, kappa of the source has to be the same as that of the initial profile (no kappa matching done)"); + if(opts_init.src_type == src_t::matching) - src_dry_distros_matching(dt); + src_dry_distros_matching(sdd); if(opts_init.src_type == src_t::simple) - src_dry_distros_simple(dt); + src_dry_distros_simple(sdd); } }; }; diff --git a/src/impl/particles_impl_src_dry_distros_matching.ipp b/src/impl/particles_impl_src_dry_distros_matching.ipp index 02e93d383..9833865e7 100644 --- a/src/impl/particles_impl_src_dry_distros_matching.ipp +++ b/src/impl/particles_impl_src_dry_distros_matching.ipp @@ -45,11 +45,21 @@ namespace libcloudphxx // create new aerosol particles based on a size distribution // if any SDs with dry radius similar to the one to be added are present, // we increase their multiplicity instead of adding new SDs + // TODO: make it work for sdd.size()>1 template - void particles_t::impl::src_dry_distros_matching(const real_t &dt) + void particles_t::impl::src_dry_distros_matching(const src_dry_distros_t &sdd) { + auto p_sdd = sdd.cbegin(); + + // add the source only once every number of steps + assert(get<2>(p_sdd->second) > 0); + if(src_stp_ctr % get<2>(p_sdd->second) != 0) return; + + const real_t sup_dt = get<2>(p_sdd->second) * opts_init.dt; + // set number of SDs to init; use count_num as storage - init_count_num_src(opts_init.src_sd_conc); + init_count_num_src(get<1>(p_sdd->second)); + // -------- TODO: match not only sizes of old particles, but also kappas and chemical composition... -------- @@ -92,9 +102,9 @@ namespace libcloudphxx // analyze distribution to get rd_min and max needed for bin sizes // TODO: this could be done once at the beginning of the simulation dist_analysis_sd_conc( - *(opts_init.src_dry_distros.begin()->second), - opts_init.src_sd_conc, - dt + *get<0>(p_sdd->second), + get<1>(p_sdd->second), + sup_dt ); // --- see how many already existing SDs match size bins --- @@ -209,16 +219,25 @@ namespace libcloudphxx n_part -= count_bins; hskpng_resize_npart(); - // init other peoperties of SDs that didnt have a match + // init other properties of SDs that didnt have a match init_kappa( - opts_init.src_dry_distros.begin()->first - ); + p_sdd->first.kappa + ); + if (opts_init.ice_switch) + { + init_insol_dry_sizes(p_sdd->first.rd_insol); + init_a_c_rho_ice(); + if (! opts_init.time_dep_ice_nucl) + { + init_T_freeze(); + } + } if(opts_init.diag_incloud_time) init_incloud_time(); init_n_sd_conc( - *(opts_init.src_dry_distros.begin()->second) + *get<0>(p_sdd->second) ); // TODO: document that n_of_lnrd_stp is expected! // init rw @@ -370,7 +389,7 @@ namespace libcloudphxx // init n of the copied SDs, but using the src distribution init_n_sd_conc( - *(opts_init.src_dry_distros.begin()->second) + *get<0>(p_sdd->second) ); // TODO: document that n_of_lnrd_stp is expected! // add the just-initialized multiplicities to the old ones diff --git a/src/impl/particles_impl_src_dry_distros_simple.ipp b/src/impl/particles_impl_src_dry_distros_simple.ipp index 8b7d733b8..3ca1dafbd 100644 --- a/src/impl/particles_impl_src_dry_distros_simple.ipp +++ b/src/impl/particles_impl_src_dry_distros_simple.ipp @@ -12,19 +12,31 @@ namespace libcloudphxx { namespace lgrngn { + using std::get; + // create new aerosol particles based on a size distribution template - void particles_t::impl::src_dry_distros_simple(const real_t &dt) + void particles_t::impl::src_dry_distros_simple(const src_dry_distros_t &sdd) { + // We assume that sdd size is 1 + // TODO: add a loop to allow sdd.size>1 + auto p_sdd = sdd.cbegin(); + + // add the source only once every number of steps + assert(get<2>(p_sdd->second) > 0); + if(src_stp_ctr % get<2>(p_sdd->second) != 0) return; + + const real_t sup_dt = get<2>(p_sdd->second) * opts_init.dt; + // set number of SDs to init; use count_num as storage - init_count_num_src(opts_init.src_sd_conc); + init_count_num_src(get<1>(p_sdd->second)); // analyze distribution to get rd_min and max needed for bin sizes // TODO: this could be done once at the beginning of the simulation dist_analysis_sd_conc( - *(opts_init.src_dry_distros.begin()->second), - opts_init.src_sd_conc, - dt + *(get<0>(p_sdd->second)), + get<1>(p_sdd->second), + sup_dt ); namespace arg = thrust::placeholders; @@ -42,17 +54,26 @@ namespace libcloudphxx init_ijk(); init_dry_sd_conc(); - // init other peoperties of SDs that didnt have a match + // init other properties of SDs that didnt have a match init_kappa( - opts_init.src_dry_distros.begin()->first - ); + p_sdd->first.kappa + ); + if (opts_init.ice_switch) + { + init_insol_dry_sizes(p_sdd->first.rd_insol); + init_a_c_rho_ice(); + if (! opts_init.time_dep_ice_nucl) + { + init_T_freeze(); + } + } if(opts_init.diag_incloud_time) init_incloud_time(); init_n_sd_conc( - *(opts_init.src_dry_distros.begin()->second) - ); // TODO: document that n_of_lnrd_stp is expected! + *get<0>(p_sdd->second) + ); // init rw init_wet(); diff --git a/src/impl/particles_impl_src_dry_sizes.ipp b/src/impl/particles_impl_src_dry_sizes.ipp index 8b44eee62..c3d461cbd 100644 --- a/src/impl/particles_impl_src_dry_sizes.ipp +++ b/src/impl/particles_impl_src_dry_sizes.ipp @@ -12,25 +12,28 @@ namespace libcloudphxx { // initialize SD parameters with dry_radius-concentration pairs (i.e. dry_sizes) template - void particles_t::impl::src_dry_sizes(const real_t &dt) + void particles_t::impl::src_dry_sizes(const src_dry_sizes_t &sds) { - using dry_sizes_t = typename opts_init_t::dry_sizes_t; - using kappa_t = typename dry_sizes_t::key_type; - using size_number_t = typename dry_sizes_t::mapped_type; - //using conc_multi_t = typename size_number_t::mapped_type; - - // loop over kappas - for (typename dry_sizes_t::const_iterator dsi = opts_init.src_dry_sizes.begin(); dsi != opts_init.src_dry_sizes.end(); ++dsi) + // loop over (kappa, rd_insol) pairs + // for (typename dry_sizes_t::const_iterator dsi = opts.src_dry_sizes.begin(); dsi != opts.src_dry_sizes.end(); ++dsi) + for (auto dsi = sds.cbegin(); dsi != sds.cend(); ++dsi) { - const kappa_t &kappa(dsi->first); - const size_number_t &size_number_map(dsi->second); + const real_t &kappa(dsi->first.kappa); + const real_t &rd_insol(dsi->first.rd_insol); + const auto &size_number_map(dsi->second); - // loop over the "size : {concentration per second, multiplicity}" pairs for this kappa - for (typename size_number_t::const_iterator sni = size_number_map.begin(); sni != size_number_map.end(); ++sni) + // loop over the "size : {concentration per second, multiplicity, supstp}" for this (kappa, rd_insol) pair + for (auto sni = size_number_map.cbegin(); sni != size_number_map.cend(); ++sni) { + // add the source only once every number of steps + assert(get<2>(sni->second) > 0); + if(src_stp_ctr % get<2>(sni->second) != 0) continue; + + const real_t sup_dt = get<2>(sni->second) * opts_init.dt; + // init number of SDs of this kappa in cells - init_count_num_src(sni->second.second); + init_count_num_src(get<1>(sni->second)); // update no of particles // TODO: move to a separate function @@ -45,11 +48,21 @@ namespace libcloudphxx // initialising dry radii (needs ijk) init_dry_dry_sizes(sni->first); - // init kappa + // init kappa and rd_insol init_kappa(kappa); + + if (opts_init.ice_switch) + { + init_insol_dry_sizes(rd_insol); + init_a_c_rho_ice(); + if (! opts_init.time_dep_ice_nucl) + { + init_T_freeze(); + } + } // init multiplicities - init_n_dry_sizes(sni->second.first*dt, sni->second.second); + init_n_dry_sizes(get<0>(sni->second)*sup_dt, get<1>(sni->second)); // initialising wet radii init_wet(); diff --git a/src/impl/particles_impl_update_th_rv.ipp b/src/impl/particles_impl_update_th_rv.ipp index 833ebf480..9d49647fb 100644 --- a/src/impl/particles_impl_update_th_rv.ipp +++ b/src/impl/particles_impl_update_th_rv.ipp @@ -13,6 +13,7 @@ namespace libcloudphxx { namespace detail { + // change of th during condensation template struct dth //: thrust::unary_function&, real_t> { @@ -29,28 +30,79 @@ namespace libcloudphxx return drv * common::theta_dry::d_th_d_rv(T, th) / si::kelvins; } }; + // change of th during deposition + template + struct dth_dep //: thrust::unary_function&, real_t> + { + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple &tpl) const + { + const quantity + drv = thrust::get<0>(tpl); + const quantity + T = thrust::get<1>(tpl) * si::kelvins; + const quantity + th = thrust::get<2>(tpl) * si::kelvins; + + return drv * common::theta_dry::d_th_d_rv_dep(T, th) / si::kelvins; + } + }; + + // change of th during freezing + template + struct dth_freezing //: thrust::unary_function&, real_t> + { + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple &tpl) const + { + const quantity + drw = thrust::get<0>(tpl); + const quantity + T = thrust::get<1>(tpl) * si::kelvins; + const quantity + th = thrust::get<2>(tpl) * si::kelvins; + + return drw * common::theta_dry::d_th_d_rw_freeze(T, th) / si::kelvins; + } + }; }; - // update th and rv according to change in 3rd specific wet moments + // update th and rv after condensation / deposition according to change in 3rd specific wet moments // particles have to be sorted template void particles_t::impl::update_th_rv( - thrust_device::vector &drv // change in water vapor mixing ratio + thrust_device::vector &drv, // change in water vapor mixing ratio + phase_change phase // enum for cond/dep, the default is condensation ) { if(!sorted) throw std::runtime_error("libcloudph++: update_th_rv called on an unsorted set"); nancheck(drv, "update_th_rv: input drv"); - // multiplying specific 3rd moms diff by -rho_w*4/3*pi - thrust::transform( - drv.begin(), drv.end(), // input - 1st arg - thrust::make_constant_iterator( // input - 2nd arg - - common::moist_air::rho_w() / si::kilograms * si::cubic_metres - * real_t(4./3) * pi() - ), - drv.begin(), // output - thrust::multiplies() - ); + if (phase == phase_change::condensation) + { + // multiplying specific 3rd moms diff by -rho_w*4/3*pi + thrust::transform( + drv.begin(), drv.end(), // input - 1st arg + thrust::make_constant_iterator( // input - 2nd arg + - common::moist_air::rho_w() / si::kilograms * si::cubic_metres + * real_t(4./3) * pi() + ), + drv.begin(), // output + thrust::multiplies() + ); + } + else if (phase == phase_change::deposition) + { + // drv is already multiplied by 4/3 * pi * rho_i + thrust::transform( + drv.begin(), drv.end(), // input - 1st arg + thrust::make_constant_iterator( // input - 2nd arg + - real_t(1) + ), + drv.begin(), // output + thrust::multiplies() + ); + } // updating rv assert(*thrust::min_element(rv.begin(), rv.end()) >= 0); @@ -72,11 +124,13 @@ namespace libcloudphxx > > zip_it_t; // apply dth + if (phase == phase_change::condensation) + { thrust::transform( th.begin(), th.end(), // input - 1st arg thrust::make_transform_iterator( - zip_it_t(thrust::make_tuple( - drv.begin(), // + zip_it_t(thrust::make_tuple( + drv.begin(), // T.begin(), // dth = drv * d_th_d_rv(T, th) th.begin() // )), @@ -86,9 +140,74 @@ namespace libcloudphxx thrust::plus() ); } + else if (phase == phase_change::deposition) + { + thrust::transform( + th.begin(), th.end(), // input - 1st arg + thrust::make_transform_iterator( + zip_it_t(thrust::make_tuple( + drv.begin(), // + T.begin(), // dth = drv * d_th_d_rv(T, th) + th.begin() // + )), + detail::dth_dep() + ), + th.begin(), // output + thrust::plus() + ); + } + + } nancheck(th, "update_th_rv: th after update"); } + // update th for freezing + // particles have to be sorted + template + void particles_t::impl::update_th_freezing( + thrust_device::vector &drw // change of specific 3rd moment of liquid per cell + ) + { + if(!sorted) throw std::runtime_error("libcloudph++: update_th_freezing called on an unsorted set"); + nancheck(drw, "update_th_freezing: input drw"); + + // Calculating the change of liquid mixing ratio per cell (multiplying specific 3rd mom by rho_w*4/3*pi) + thrust::transform( + drw.begin(), drw.end(), // input - 1st arg + thrust::make_constant_iterator( // input - 2nd arg + common::moist_air::rho_w() / si::kilograms * si::cubic_metres + * real_t(4./3) * pi() + ), + drw.begin(), // output + thrust::multiplies() + ); + + // updating th + { + typedef thrust::zip_iterator::iterator, + typename thrust_device::vector::iterator, + typename thrust_device::vector::iterator + > > zip_it_t; + + // apply dth + thrust::transform( + th.begin(), th.end(), // input - 1st arg + thrust::make_transform_iterator( + zip_it_t(thrust::make_tuple( + drw.begin(), + T.begin(), + th.begin() + )), + detail::dth_freezing() + ), + th.begin(), // output + thrust::plus() + ); + } + nancheck(th, "update_th_freezing: th after update"); + } + // update particle-specific cell state // particles have to be sorted template diff --git a/src/particles.tpp b/src/particles.tpp index a3ad121b1..cc085c444 100644 --- a/src/particles.tpp +++ b/src/particles.tpp @@ -64,6 +64,9 @@ #include "impl/particles_impl_init_dry_const_multi.ipp" #include "impl/particles_impl_init_dry_dry_sizes.ipp" #include "impl/particles_impl_init_kappa.ipp" +#include "impl/particles_impl_init_insol_dry_sizes.ipp" +#include "impl/particles_impl_init_T_freeze.ipp" +#include "impl/particles_impl_init_a_c_rho_ice.ipp" #include "impl/particles_impl_init_incloud_time.ipp" #include "impl/particles_impl_init_n.ipp" #include "impl/particles_impl_init_wet.ipp" @@ -104,6 +107,8 @@ #include "impl/particles_impl_cond_common.ipp" #include "impl/particles_impl_cond.ipp" #include "impl/particles_impl_cond_sstp.ipp" +#include "impl/particles_impl_ice_dep.ipp" +#include "impl/particles_impl_ice_nucl_melt.ipp" #include "impl/particles_impl_sedi.ipp" #include "impl/particles_impl_subs.ipp" #include "impl/particles_impl_coal.ipp" diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp index 29a2df0e9..369592b99 100644 --- a/src/particles_diag.ipp +++ b/src/particles_diag.ipp @@ -74,6 +74,24 @@ namespace libcloudphxx } }; + template + struct precip_rate_ice + { + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple &tpl) // tpl is a tuple (a, c, rho, vt) + { + return real_t(4./3) + #if !defined(__NVCC__) + * pi() + #else + * CUDART_PI + #endif + * thrust::get<0>(tpl) * thrust::get<0>(tpl) * thrust::get<1>(tpl) // a^2 * c + * thrust::get<2>(tpl) // rho + * thrust::get<3>(tpl); // vt + } + }; + template struct RH_minus_Sc { @@ -104,6 +122,25 @@ namespace libcloudphxx return sqrt(rw2); } }; + + template + class ice_mass + { + public: + BOOST_GPU_ENABLED + real_t operator()(const thrust::tuple &tpl) // tpl is a tuple (a, c, rho) + { + return real_t(4./3) + #if !defined(__NVCC__) + * pi() + #else + * CUDART_PI + #endif + * thrust::get<0>(tpl) * thrust::get<0>(tpl) * thrust::get<1>(tpl) // a^2 * c + * thrust::get<2>(tpl); // rho + } + }; + } // records pressure @@ -209,6 +246,20 @@ namespace libcloudphxx pimpl->moms_rng(pow(r_min, 2), pow(r_max, 2), pimpl->rw2.begin(), false); } + // selects particles with (ice_a >= a_min && ice_a < a_max) + template + void particles_t::diag_ice_a_rng(const real_t &a_min, const real_t &a_max) + { + pimpl->moms_rng(a_min, a_max, pimpl->ice_a.begin(), false); + } + + // selects particles with (ice_c >= c_min && ice_c < c_max) + template + void particles_t::diag_ice_c_rng(const real_t &c_min, const real_t &c_max) + { + pimpl->moms_rng(c_min, c_max, pimpl->ice_c.begin(), false); + } + // selects particles with (kpa >= kpa_min && kpa < kpa_max) template void particles_t::diag_kappa_rng(const real_t &kpa_min, const real_t &kpa_max) @@ -216,6 +267,20 @@ namespace libcloudphxx pimpl->moms_rng(kpa_min, kpa_max, pimpl->kpa.begin(), false); } + // selects ice particles + template + void particles_t::diag_ice() + { + pimpl->moms_gt0(pimpl->ice_a.begin()); // ice_a greater than 0 + } + + // selects water particles + template + void particles_t::diag_water() + { + pimpl->moms_eq0(pimpl->ice_a.begin()); // ice_a equal to 0 + } + // selects particles with (r_d >= r_min && r_d < r_max) from particles previously selected template void particles_t::diag_dry_rng_cons(const real_t &r_min, const real_t &r_max) @@ -236,6 +301,20 @@ namespace libcloudphxx pimpl->moms_rng(pow(r_min, 2), pow(r_max, 2), pimpl->rw2.begin(), true); } + // selects particles with (ice_a >= a_min && ice_a < a_max) from particles previously selected + template + void particles_t::diag_ice_a_rng_cons(const real_t &a_min, const real_t &a_max) + { + pimpl->moms_rng(a_min, a_max, pimpl->ice_a.begin(), true); + } + + // selects particles with (ice_c >= c_min && ice_c < c_max) from particles previously selected + template + void particles_t::diag_ice_c_rng_cons(const real_t &c_min, const real_t &c_max) + { + pimpl->moms_rng(c_min, c_max, pimpl->ice_c.begin(), true); + } + // selects particles with (kpa >= kpa_min && kpa < kpa_max) from particles previously selected template void particles_t::diag_kappa_rng_cons(const real_t &kpa_min, const real_t &kpa_max) @@ -243,6 +322,20 @@ namespace libcloudphxx pimpl->moms_rng(kpa_min, kpa_max, pimpl->kpa.begin(), true); } + // selects ice particles from particles previously selected + template + void particles_t::diag_ice_cons() + { + pimpl->moms_gt0(pimpl->ice_a.begin(), true); // ice_a greater than 0 + } + + // selects water particles from particles previously selected + template + void particles_t::diag_water_cons() + { + pimpl->moms_eq0(pimpl->ice_a.begin(), true); // ice_a equal to 0 + } + // selects particles with RH >= Sc (Sc - critical supersaturation) template void particles_t::diag_RH_ge_Sc() @@ -315,6 +408,31 @@ namespace libcloudphxx pimpl->moms_calc(pimpl->rw2.begin(), n/2.); } + // computes n-th moment of the ice equatorial radius spectrum for the selected particles + template + void particles_t::diag_ice_a_mom(const int &n) + { + pimpl->moms_calc(pimpl->ice_a.begin(), n); + } + + // computes n-th moment of the ice polar radius spectrum for the selected particles + template + void particles_t::diag_ice_c_mom(const int &n) + { + pimpl->moms_calc(pimpl->ice_c.begin(), n); + } + + // computes ice mixing ratio + template + void particles_t::diag_ice_mix_ratio() + { + pimpl->moms_calc(thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple(pimpl->ice_a.begin(), pimpl->ice_c.begin(), pimpl->ice_rho.begin())), + detail::ice_mass() + ), + real_t(1)); + } + // compute n-th moment of kappa for selected particles template void particles_t::diag_kappa_mom(const int &n) @@ -445,7 +563,23 @@ namespace libcloudphxx // copy back stored vterm thrust::copy(tmp_vt.begin(), tmp_vt.end(), pimpl->vt.begin()); - } + } + + // compute 1st (non-specific) moment of ice_mass * vt of all SDs + template + void particles_t::diag_precip_rate_ice_mass() + { + // updating terminal velocities + pimpl->hskpng_vterm_all(); + + pimpl->moms_all(); // we need this here, because hskpng_vterm modifies tmp_device_real_part, which is used as n_filtered in moms_calc + pimpl->moms_calc(thrust::make_transform_iterator( + thrust::make_zip_iterator(thrust::make_tuple(pimpl->ice_a.begin(), pimpl->ice_c.begin(), + pimpl->ice_rho.begin(), pimpl->vt.begin())), + detail::precip_rate_ice() + ), + real_t(1), false); + } // get max rw in each cell template diff --git a/src/particles_multi_gpu_diag.ipp b/src/particles_multi_gpu_diag.ipp index 2a096964e..efc05db16 100644 --- a/src/particles_multi_gpu_diag.ipp +++ b/src/particles_multi_gpu_diag.ipp @@ -79,6 +79,18 @@ namespace libcloudphxx pimpl->mcuda_run(&particles_t::diag_kappa_rng, r_mi, r_mx); } + template + void particles_t::diag_ice() + { + pimpl->mcuda_run(&particles_t::diag_ice); + } + + template + void particles_t::diag_water() + { + pimpl->mcuda_run(&particles_t::diag_water); + } + template void particles_t::diag_dry_rng_cons( const real_t &r_mi, const real_t &r_mx @@ -103,6 +115,18 @@ namespace libcloudphxx pimpl->mcuda_run(&particles_t::diag_kappa_rng_cons, r_mi, r_mx); } + template + void particles_t::diag_ice_cons() + { + pimpl->mcuda_run(&particles_t::diag_ice_cons); + } + + template + void particles_t::diag_water_cons() + { + pimpl->mcuda_run(&particles_t::diag_water_cons); + } + template void particles_t::diag_dry_mom(const int &k) { diff --git a/src/particles_step.ipp b/src/particles_step.ipp index c80c80b24..db530eaf6 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -180,6 +180,10 @@ namespace libcloudphxx if (pimpl->opts_init.diag_incloud_time) pimpl->update_incloud_time(pimpl->dt); + // ice nucleation/melting + if (pimpl->opts_init.ice_switch && opts.ice_nucl) + pimpl->ice_nucl_melt(pimpl->dt); + // condensation/evaporation if (opts.cond) { @@ -191,13 +195,16 @@ namespace libcloudphxx if(pimpl->opts_init.exact_sstp_cond && pimpl->sstp_cond > 1) // apply substeps per-particle logic { + if (pimpl->opts_init.ice_switch) + throw std::runtime_error("libcloudph++: deposition works only with per-cell substepping"); + for (int step = 0; step < pimpl->sstp_cond; ++step) { pimpl->sstp_step_exact(step); if(opts.turb_cond) pimpl->sstp_step_ssp(pimpl->dt / pimpl->sstp_cond); - pimpl->cond_sstp(pimpl->dt / pimpl->sstp_cond, opts.RH_max, opts.turb_cond, step); - } + pimpl->cond_sstp(pimpl->dt / pimpl->sstp_cond, opts.RH_max, opts.turb_cond, step); + } // copy sstp_tmp_rv and th to rv and th pimpl->update_state(pimpl->rv, pimpl->sstp_tmp_rv); pimpl->update_state(pimpl->th, pimpl->sstp_tmp_th); @@ -210,8 +217,14 @@ namespace libcloudphxx pimpl->sstp_step(step); if(opts.turb_cond) pimpl->sstp_step_ssp(pimpl->dt / pimpl->sstp_cond); - pimpl->hskpng_Tpr(); + pimpl->hskpng_Tpr(); pimpl->cond(pimpl->dt / pimpl->sstp_cond, opts.RH_max, opts.turb_cond, step); + if (pimpl->opts_init.ice_switch) + { + if (opts.turb_cond) + throw std::runtime_error("libcloudph++:deposition doesnt work with turb_cond"); + pimpl->ice_dep(pimpl->dt / pimpl->sstp_cond, opts.RH_max, step); + } } } @@ -411,11 +424,8 @@ namespace libcloudphxx // sanity check if (pimpl->opts_init.src_type == src_t::off) throw std::runtime_error("libcloudph++: aerosol source was switched off in opts_init"); - // introduce new particles with the given time interval - if(pimpl->src_stp_ctr % pimpl->opts_init.supstp_src == 0) - { - pimpl->src(pimpl->opts_init.supstp_src * pimpl->dt); - } + // introduce new particles + pimpl->src(opts.src_dry_distros, opts.src_dry_sizes); } // aerosol relaxation, in sync since it changes th/rv diff --git a/tests/mpi/mpi_adve_test.cpp b/tests/mpi/mpi_adve_test.cpp index 7f07a2060..b1f2c92f3 100644 --- a/tests/mpi/mpi_adve_test.cpp +++ b/tests/mpi/mpi_adve_test.cpp @@ -121,7 +121,7 @@ void test(backend_t backend, std::string back_name, int ndims, bool dir, int n_d */ opts_init.dry_distros.emplace( - 0.001, // kappa + libcloudphxx::lgrngn::kappa_rd_insol_t{double(0.001), double(0.)}, // kappa, rd_insol std::make_shared>() // distribution ); diff --git a/tests/python/physics/coalescence_golovin.py b/tests/python/physics/coalescence_golovin.py index 9576da3eb..ff35efe92 100644 --- a/tests/python/physics/coalescence_golovin.py +++ b/tests/python/physics/coalescence_golovin.py @@ -68,8 +68,9 @@ def golovin(v,t,n0,v0,b): rv = 0.01 * np.ones((1,)) kappa = 1e-10 +rd_insol = 0. -opts_init.dry_distros = {kappa:expvolumelnr} +opts_init.dry_distros = {(kappa, rd_insol):expvolumelnr} opts_init.kernel = lgrngn.kernel_t.golovin opts_init.terminal_velocity = lgrngn.vt_t.beard77 diff --git a/tests/python/physics/coalescence_hall_davis_no_waals.py b/tests/python/physics/coalescence_hall_davis_no_waals.py index 66f66b9f6..e28db195c 100644 --- a/tests/python/physics/coalescence_hall_davis_no_waals.py +++ b/tests/python/physics/coalescence_hall_davis_no_waals.py @@ -52,8 +52,9 @@ def expvolumelnr(lnr): rv = 0.01 * np.ones((opts_init.nx, opts_init.nz)) kappa = 0 #1e-10 +rd_insol = 0. -opts_init.dry_distros = {kappa:expvolumelnr} +opts_init.dry_distros = {(kappa, rd_insol):expvolumelnr} opts_init.sd_conc = pow(2,14) opts_init.n_sd_max = opts_init.sd_conc diff --git a/tests/python/physics/lgrngn_cond.py b/tests/python/physics/lgrngn_cond.py index 0562d0259..5eef756f7 100644 --- a/tests/python/physics/lgrngn_cond.py +++ b/tests/python/physics/lgrngn_cond.py @@ -29,9 +29,11 @@ def lognormal(lnr): opts_init = lgrngn.opts_init_t() kappa = .61 -opts_init.dry_distros = {kappa:lognormal} +rd_insol = 0. +opts_init.dry_distros = {(kappa, rd_insol):lognormal} opts_init.coal_switch = False opts_init.sedi_switch = False +opts_init.ice_switch = False opts_init.RH_max = 0.999 # to comply with the assert(RH<1) at init opts_init.dt = 1 opts_init.sd_conc = int(1e2) @@ -44,6 +46,7 @@ def lognormal(lnr): opts.cond = True opts.coal = False opts.chem = False +opts.ice_nucl = False #expected theta and rv after condensation: exp_th = { True : 306.9, # constp diff --git a/tests/python/physics/lgrngn_cond_substepping.py b/tests/python/physics/lgrngn_cond_substepping.py index 55b29eb3c..0b453f9b8 100644 --- a/tests/python/physics/lgrngn_cond_substepping.py +++ b/tests/python/physics/lgrngn_cond_substepping.py @@ -41,9 +41,10 @@ def lognormal2(lnr): opts = lgrngn.opts_t() opts_init = lgrngn.opts_init_t() -opts_init.dry_distros = {.61:lognormal, 1.28:lognormal2} # normal mode + GCCNs +opts_init.dry_distros = {(.61,0.):lognormal, (1.28,0.):lognormal2} # normal mode + GCCNs opts_init.coal_switch = False opts_init.sedi_switch = False +opts_init.ice_switch = False opts_init.RH_max = 1.0001 opts_init.dt = 1 opts_init.sd_conc = int(1e3) @@ -56,6 +57,7 @@ def lognormal2(lnr): opts.cond = True opts.coal = False opts.chem = False +opts.ice_nucl = False #opts.RH_max = 1.005 #expected theta and rv after condensation: diff --git a/tests/python/physics/puddle.py b/tests/python/physics/puddle.py index 1ef20f2dd..e80b93525 100644 --- a/tests/python/physics/puddle.py +++ b/tests/python/physics/puddle.py @@ -17,9 +17,11 @@ def lognormal(lnr): Opts_init = lgrngn.opts_init_t() kappa = .61 -Opts_init.dry_distros = {kappa:lognormal} +rd_insol = 0. +Opts_init.dry_distros = {(kappa, rd_insol):lognormal} Opts_init.coal_switch = False Opts_init.sedi_switch = True +Opts_init.ice_switch = False Opts_init.terminal_velocity = lgrngn.vt_t.beard76 Opts_init.dt = 1 @@ -44,6 +46,7 @@ def lognormal(lnr): Opts.coal = False Opts.chem = False Opts.rcyc = False +Opts.ice_nucl = False Rhod = 1. * np.ones((Opts_init.nx, Opts_init.nz)) Th = 300. * np.ones((Opts_init.nx, Opts_init.nz)) @@ -68,7 +71,10 @@ def lognormal(lnr): assert(tab_out[0][0] == 0.) -puddle_expected_per_cell = {0: 0.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0, 5: 0.0, 6: 0.0, 7: 0.0, 8: 7.087802417148837e-05, 9: 5.630090090571395e-06, 10: 815411.5, 11: 0.0} +puddle_expected_per_cell = {'HNO3': 0.0, 'NH3': 0.0, 'CO2': 0.0, 'SO2': 0.0, 'H2O2': 0.0, 'O3': 0.0, 'S_VI': 0.0, 'H': 0.0, + 'liquid_volume': 7.087802417148837e-05, 'dry_volume': 5.630090090571395e-06, + 'particle_number': 815411.5, 'liquid_number': 815411.5, + 'ice_mass': 0.0, 'ice_number': 0.0} for a in puddle: print(puddle[a], Opts_init.nx * puddle_expected_per_cell[a]) diff --git a/tests/python/physics/test_coal.py b/tests/python/physics/test_coal.py index f841a083a..8c2273a5a 100644 --- a/tests/python/physics/test_coal.py +++ b/tests/python/physics/test_coal.py @@ -29,8 +29,9 @@ def expvolumelnr(lnr): kappa1 = 0.1 kappa2 = 0.9 +rd_insol = 0. -opts_init.dry_distros = {kappa1:expvolumelnr, kappa2:expvolumelnr} +opts_init.dry_distros = {(kappa1, rd_insol):expvolumelnr, (kappa2, rd_insol):expvolumelnr} opts_init.sd_conc = pow(2,14) opts_init.n_sd_max = pow(2,14) diff --git a/tests/python/unit/CMakeLists.txt b/tests/python/unit/CMakeLists.txt index 14bf67684..badc92dbc 100644 --- a/tests/python/unit/CMakeLists.txt +++ b/tests/python/unit/CMakeLists.txt @@ -1,5 +1,5 @@ # non-pytest tests -foreach(test api_blk_1m api_blk_2m api_lgrngn api_common segfault_20150216 col_kernels terminal_velocities uniform_init source sstp_cond multiple_kappas adve_scheme lgrngn_subsidence sat_adj_blk_1m diag_incloud_time relax blk_1m_ice) +foreach(test api_blk_1m api_blk_2m api_lgrngn api_common segfault_20150216 col_kernels terminal_velocities uniform_init source sstp_cond multiple_kappas adve_scheme lgrngn_subsidence sat_adj_blk_1m diag_incloud_time relax blk_1m_ice ice_SD) #TODO: indicate that tests depend on the lib add_test( diff --git a/tests/python/unit/SD_removal.py b/tests/python/unit/SD_removal.py index b09b5f61c..5abba2b91 100755 --- a/tests/python/unit/SD_removal.py +++ b/tests/python/unit/SD_removal.py @@ -36,9 +36,10 @@ def expvolumelnr(lnr): r=np.exp(lnr) return n_zero * 3.*np.power(r,3)/np.power(r_zero,3)*np.exp(- np.power((r/r_zero),3)); -kappa = .01 +kappa = .01 +rd_insol = 0. -opts_init.dry_distros = {kappa:expvolumelnr} +opts_init.dry_distros = {(kappa, rd_insol):expvolumelnr} opts_init.sd_conc = 64 opts_init.n_sd_max = 64 diff --git a/tests/python/unit/adve_scheme.py b/tests/python/unit/adve_scheme.py index 9448d7fad..8c01a560a 100644 --- a/tests/python/unit/adve_scheme.py +++ b/tests/python/unit/adve_scheme.py @@ -21,8 +21,9 @@ def lognormal(lnr): ) / log(stdev) / sqrt(2*pi); kappa = .61 +rd_insol = 0. -opts_init.dry_distros = {kappa:lognormal} +opts_init.dry_distros = {(kappa, rd_insol):lognormal} opts_init.sd_conc = 50 opts_init.n_sd_max = 50 diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index 5b466a50e..e3e78fa3a 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -97,7 +97,7 @@ def test_rhs_cell(opts, name): assert flux == 0 assert dot_rr == dot_rr_old # no rain water -> no precip -#test ice, TODO: separate file/test? +#test RHS cellwise ice def test_rhs_cell_ice(opts, name): print("Testing RHS cellwise ice with " + name) th = arr_t([230.]) #testing ice physics diff --git a/tests/python/unit/api_lgrngn.py b/tests/python/unit/api_lgrngn.py index 40a7709e3..4da852d4c 100644 --- a/tests/python/unit/api_lgrngn.py +++ b/tests/python/unit/api_lgrngn.py @@ -24,8 +24,9 @@ def lognormal(lnr): kappa1 = .61 kappa2 = 1.28 kappa3 = 0.8 +rd_insol = 0.5e-6 rho_stp = 1.2248 -opts_init.dry_distros = {kappa1:lognormal, kappa2:lognormal} +opts_init.dry_distros = {(kappa1, rd_insol): lognormal} opts_init.kernel = lgrngn.kernel_t.geometric opts_init.terminal_velocity = lgrngn.vt_t.beard76 opts_init.adve_scheme = lgrngn.as_t.euler @@ -35,9 +36,7 @@ def lognormal(lnr): opts_init.rng_seed = 396 opts_init.rng_seed_init = 456 opts_init.rng_seed_init_switch = True -opts_init.src_dry_distros = {kappa1:lognormal} opts_init.rlx_dry_distros = {kappa1:[lognormal, [0,2], [0,opts_init.dz]]} -opts_init.src_sd_conc = 64 opts_init.src_z1 = opts_init.dz opts_init.diag_incloud_time = True @@ -69,6 +68,7 @@ def lognormal(lnr): print("coal_switch = ", opts_init.coal_switch) print("sedi_switch = ", opts_init.sedi_switch) print("subs_switch = ", opts_init.subs_switch) +print("ice_switch = ", opts_init.ice_switch) print("src_type = ", opts_init.src_type) print("exact_sstp_cond = ", opts_init.exact_sstp_cond) @@ -106,6 +106,8 @@ def lognormal(lnr): print("chem_gas[SO2] = ", opts.chem_gas[lgrngn.chem_species_t.SO2]) print("chem_gas = ", opts.chem_gas) +opts.src_dry_distros = {(kappa1, rd_insol):(lognormal, 64, 1)} + # --------- test runs ----------- # ---------- @@ -241,8 +243,8 @@ def lognormal(lnr): # 0D dry_sizes init with two kappas print("0D dry sizes") opts_init.dry_distros = dict() -opts_init.dry_sizes = {kappa1 : {1.e-6 : [30. * rho_stp, 15], 15.e-6 : [10. * rho_stp, 10]}, - kappa2 : {1.2e-6 : [20. * rho_stp, 10], 12.e-6 : [15. * rho_stp, 15]}} +opts_init.dry_sizes = {(kappa1, rd_insol) : {1.e-6 : [30. * rho_stp, 15], 15.e-6 : [10. * rho_stp, 10]}, + (kappa2, rd_insol) : {1.2e-6 : [20. * rho_stp, 10], 12.e-6 : [15. * rho_stp, 15]}} sd_conc_old = opts_init.sd_conc opts_init.sd_conc = 0 @@ -301,14 +303,14 @@ def lognormal(lnr): # go back to distros init opts_init.sd_conc = sd_conc_old opts_init.dry_sizes = dict() -opts_init.dry_distros = {kappa1:lognormal, kappa2:lognormal} +opts_init.dry_distros = {(kappa1, rd_insol):lognormal, (kappa2, rd_insol):lognormal} # ---------- # 0D dry_sizes + sd_conc init print("0D dry_sizes + sd_conc") -opts_init.dry_sizes = {kappa3 : {1.e-6 : [30. * rho_stp, 15], 15.e-6 : [10. * rho_stp, 5]}} +opts_init.dry_sizes = {(kappa3, rd_insol) : {1.e-6 : [30. * rho_stp, 15], 15.e-6 : [10. * rho_stp, 5]}} prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) @@ -326,7 +328,7 @@ def lognormal(lnr): # ---------- # 0D dry_sizes + sd_conc + tail print("0D dry_sizes + sd_conc + tail") -opts_init.dry_sizes = {kappa3 : {1.e-6 : [30. * rho_stp, 15], 15.e-6 : [10. * rho_stp, 5]}} +opts_init.dry_sizes = {(kappa3, rd_insol): {1.e-6 : [30. * rho_stp, 15], 15.e-6 : [10. * rho_stp, 5]}} opts_init.sd_conc_large_tail = 1 prtcls = lgrngn.factory(backend, opts_init) @@ -346,7 +348,7 @@ def lognormal(lnr): # ---------- # 0D dry_sizes + const_multi init print("0D dry_sizes + const_multi") -opts_init.dry_sizes = {kappa3 : {1.e-6 : [30. * rho_stp, 15], 15.e-6 : [10. * rho_stp, 5]}} +opts_init.dry_sizes = {(kappa3, rd_insol) : {1.e-6 : [30. * rho_stp, 15], 15.e-6 : [10. * rho_stp, 5]}} opts_init.sd_conc = 0 prtcls_per_cell = 2 * n_tot / rho_stp #rhod=1; 2* because of two distributions opts_init.sd_const_multi = int(prtcls_per_cell / 64) @@ -386,6 +388,29 @@ def lognormal(lnr): opts_init.turb_coal_switch=False opts.turb_coal=False +# ---------- +# 0D ice +print("0D ice") +th = arr_t([263.]) +opts_init.ice_switch = True +opts_init.coal_switch = False +opts.ice_nucl = True +prtcls = lgrngn.factory(backend, opts_init) +prtcls.init(th, rv, rhod) +prtcls.step_sync(opts, th, rv) +prtcls.step_async(opts) + +prtcls.diag_all() +prtcls.diag_sd_conc() +assert len(frombuffer(prtcls.outbuf())) == 1 +print(frombuffer(prtcls.outbuf())) +assert (frombuffer(prtcls.outbuf()) > 0).all() +assert sum(frombuffer(prtcls.outbuf())) == 64 + +opts_init.ice_switch=False +opts.ice_nucl=False +opts_init.coal_switch = True + # ---------- @@ -617,7 +642,7 @@ def lognormal(lnr): # 3D dry_sizes init print("3D dry sizes") opts_init.dry_distros = dict() -opts_init.dry_sizes = {kappa1 : {1.e-6 : [30./ cell_vol * rho_stp, 30], 15.e-6 : [10. / cell_vol * rho_stp, 10]}} +opts_init.dry_sizes = {(kappa1, rd_insol) : {1.e-6 : [30./ cell_vol * rho_stp, 30], 15.e-6 : [10. / cell_vol * rho_stp, 10]}} prtcls = lgrngn.factory(backend, opts_init) prtcls.init(th, rv, rhod) @@ -649,8 +674,8 @@ def lognormal(lnr): # ---------- # 3D dry_sizes + sd_conc init print("3D dry_sizes + sd_conc") -opts_init.dry_distros = {kappa1:lognormal, kappa2:lognormal} -opts_init.dry_sizes = {kappa1 : {1.e-6 : [30./ cell_vol * rho_stp, 15], 15.e-6 : [10. / cell_vol * rho_stp, 5]}} +opts_init.dry_distros = {(kappa1, rd_insol):lognormal, (kappa2, rd_insol):lognormal} +opts_init.dry_sizes = {(kappa1, rd_insol) : {1.e-6 : [30./ cell_vol * rho_stp, 15], 15.e-6 : [10. / cell_vol * rho_stp, 5]}} opts_init.sd_conc = sd_conc_old opts_init.sd_const_multi = 0 @@ -664,8 +689,12 @@ def lognormal(lnr): # test if get_attr work and if kappas are set correctly kappa = asarray(prtcls.get_attr("kappa")) -assert (kappa[:(32*opts_init.nx*opts_init.ny*opts_init.nz)] == kappa2).all() -assert (kappa[(32*opts_init.nx*opts_init.ny*opts_init.nz):] == kappa1).all() +# assert (kappa[:(32*opts_init.nx*opts_init.ny*opts_init.nz)] == kappa2).all() +# assert (kappa[(32*opts_init.nx*opts_init.ny*opts_init.nz):] == kappa1).all() +n = 32 * opts_init.nx * opts_init.ny * opts_init.nz +assert (kappa[:n] == kappa1).all() +assert (kappa[n:2*n] == kappa2).all() +assert (kappa[2*n:] == kappa1).all() # ---------- @@ -692,7 +721,7 @@ def lognormal(lnr): # ---------- # 3D dry_sizes + const_multi init print("3D dry_sizes + const_multi") -opts_init.dry_sizes = {kappa1 : {1.e-6 : [30./ cell_vol * rho_stp, 15], 15.e-6 : [10. / cell_vol * rho_stp, 5]}} +opts_init.dry_sizes = {(kappa1, rd_insol) : {1.e-6 : [30./ cell_vol * rho_stp, 15], 15.e-6 : [10. / cell_vol * rho_stp, 5]}} opts_init.sd_conc = 0 prtcls_per_cell = 2 * n_tot * cell_vol / rho_stp #rhod=1; 2* because of two distributions opts_init.sd_const_multi = int(prtcls_per_cell / 64) @@ -709,3 +738,28 @@ def lognormal(lnr): opts_init.sd_conc = sd_conc_old opts_init.sd_const_multi = 0 opts_init.dry_sizes = dict() + + +# ---------- +# 3D ice +print("3D ice") +th_arr = arr_t([[263., 263. ], [263., 263. ]]) +th = arr_t([th_arr, th_arr ]) +opts_init.ice_switch = True +opts_init.coal_switch = False +opts.ice_nucl = True +prtcls = lgrngn.factory(backend, opts_init) +prtcls.init(th, rv, rhod) +prtcls.step_sync(opts, th, rv) +prtcls.step_async(opts) + +prtcls.diag_all() +prtcls.diag_sd_conc() +assert len(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny +print(frombuffer(prtcls.outbuf())) +assert (frombuffer(prtcls.outbuf()) > 0).all() +assert sum(frombuffer(prtcls.outbuf())) == opts_init.nz * opts_init.nx * opts_init.ny * opts_init.sd_conc + +opts_init.ice_switch = False +opts_init.coal_switch = True +opts.ice_nucl = False diff --git a/tests/python/unit/chem_coal.py b/tests/python/unit/chem_coal.py index 84bd1d9d0..bb0be47b8 100755 --- a/tests/python/unit/chem_coal.py +++ b/tests/python/unit/chem_coal.py @@ -34,8 +34,9 @@ def expvolumelnr(lnr): return n_zero * 3.*np.power(r,3)/np.power(r_zero,3)*np.exp(- np.power((r/r_zero),3)); kappa = 0.1 +rd_insol = 0. -opts_init.dry_distros = {kappa:expvolumelnr} +opts_init.dry_distros = {(kappa, rd_insol):expvolumelnr} opts_init.sd_conc = 64 opts_init.n_sd_max = 64 diff --git a/tests/python/unit/col_kernels.py b/tests/python/unit/col_kernels.py index 9f3c1ba1b..7a433012b 100644 --- a/tests/python/unit/col_kernels.py +++ b/tests/python/unit/col_kernels.py @@ -24,13 +24,14 @@ def lognormal(lnr): ) / log(stdev) / sqrt(2*pi); kappa = .61 +rd_insol = 0. count = 0 for kernel in [lgrngn.kernel_t.geometric, lgrngn.kernel_t.geometric, lgrngn.kernel_t.long, lgrngn.kernel_t.hall, lgrngn.kernel_t.hall_davis_no_waals, lgrngn.kernel_t.golovin, lgrngn.kernel_t.onishi_hall, lgrngn.kernel_t.onishi_hall_davis_no_waals, lgrngn.kernel_t.vohl_davis_no_waals, lgrngn.kernel_t.hall_pinsky_cumulonimbus, lgrngn.kernel_t.hall_pinsky_stratocumulus]: print(kernel) opts_init = lgrngn.opts_init_t() opts_init.dt = 1 - opts_init.dry_distros = {kappa:lognormal} + opts_init.dry_distros = {(kappa, rd_insol):lognormal} opts_init.sd_conc = 50 opts_init.n_sd_max = 50 opts_init.terminal_velocity=lgrngn.vt_t.beard76 diff --git a/tests/python/unit/diag_incloud_time.py b/tests/python/unit/diag_incloud_time.py index 64a2a4f6d..a848f38df 100644 --- a/tests/python/unit/diag_incloud_time.py +++ b/tests/python/unit/diag_incloud_time.py @@ -20,7 +20,7 @@ def lognormal(lnr): opts = lgrngn.opts_t() opts_init = lgrngn.opts_init_t() -opts_init.dry_distros = {.61:lognormal, 1.28:lognormal} +opts_init.dry_distros = {(.61, 0.):lognormal, (1.28, 0.):lognormal} opts_init.coal_switch = False opts_init.sedi_switch = False opts_init.RH_max = 0.999 # to comply with the assert(RH<1) at init diff --git a/tests/python/unit/ice_SD.py b/tests/python/unit/ice_SD.py new file mode 100644 index 000000000..c22c9e19e --- /dev/null +++ b/tests/python/unit/ice_SD.py @@ -0,0 +1,63 @@ +import sys + +sys.path.insert(0, "../../bindings/python/") +sys.path.insert(0, "../../../build/bindings/python/") + +from libcloudphxx import lgrngn, common +import numpy as np +from math import exp, log, sqrt, pi + +n_tot = 60e6 +def lognormal(lnr): + mean_r = .04e-6 / 2 + stdev = 1.4 + return n_tot * exp( + -pow((lnr - log(mean_r)), 2) / 2 / pow(log(stdev),2) + ) / log(stdev) / sqrt(2*pi); + +opts_init = lgrngn.opts_init_t() +backend = lgrngn.backend_t.serial +opts = lgrngn.opts_t() + +kappa = .61 +rd_insol = 0.5e-6 +opts_init.dry_distros = {(kappa, rd_insol): lognormal} +opts_init.dt = 0.1 +opts_init.sd_conc = 100 +opts_init.n_sd_max = 1000 +opts_init.RH_max = 0.95 + +opts_init.ice_switch = True +opts_init.coal_switch = False +opts_init.sedi_switch = False +opts.ice_nucl = True +opts.cond = True +opts.coal = False +opts.adve = False +opts.sedi = False + + +for time_dep_switch in [True, False]: + print("time dependent ice nucleation = ", time_dep_switch) + p = 80000. + T = 243. + RH = 1. + rv = np.array([RH * common.r_vs(T, p)]) + th = np.array([T / common.exner(p)]) + rhod = np.array([common.rhod(p, th[0], rv[0])]) + + opts_init.time_dep_ice_nucl = time_dep_switch + prtcls = lgrngn.factory(backend, opts_init) + prtcls.init(th, rv, rhod) + for _ in range(500): + prtcls.step_sync(opts, th, rv, rhod) + prtcls.step_async(opts) + prtcls.diag_all() + prtcls.diag_ice_mix_ratio() + ri = np.frombuffer(prtcls.outbuf()) + print("ice mixing ratio ", ri[0] * 1e3, "g/kg") + print("water vapor mixing ratio ", rv[0] * 1e3, "g/kg") + assert np.isnan(ri[0]) == False + assert np.isnan(rv[0]) == False + assert rv[0] >= 0 + assert ri[0] >= 0 \ No newline at end of file diff --git a/tests/python/unit/lgrngn_adve.py b/tests/python/unit/lgrngn_adve.py index 509cfbf90..b3c518742 100644 --- a/tests/python/unit/lgrngn_adve.py +++ b/tests/python/unit/lgrngn_adve.py @@ -24,7 +24,8 @@ def lognormal(lnr): Opts_init = lgrngn.opts_init_t() kappa = .61 -Opts_init.dry_distros = {kappa:lognormal} +rd_insol = 0. +Opts_init.dry_distros = {(kappa, rd_insol):lognormal} Opts_init.coal_switch = False Opts_init.sedi_switch = False diff --git a/tests/python/unit/lgrngn_subsidence.py b/tests/python/unit/lgrngn_subsidence.py index 59e0cfccb..c55da7e32 100644 --- a/tests/python/unit/lgrngn_subsidence.py +++ b/tests/python/unit/lgrngn_subsidence.py @@ -17,7 +17,8 @@ def lognormal(lnr): Opts_init = lgrngn.opts_init_t() kappa = .61 -Opts_init.dry_distros = {kappa:lognormal} +rd_insol = 0. +Opts_init.dry_distros = {(kappa, rd_insol):lognormal} Opts_init.coal_switch = False Opts_init.sedi_switch = False Opts_init.subs_switch = True diff --git a/tests/python/unit/multiple_kappas.py b/tests/python/unit/multiple_kappas.py index 44d41f137..5f332fbac 100644 --- a/tests/python/unit/multiple_kappas.py +++ b/tests/python/unit/multiple_kappas.py @@ -33,15 +33,15 @@ def check_kappa_conc(prtcls, eps): opts_init = lgrngn.opts_init_t() kappa1 = .61 kappa2 = 1.28 +rd_insol = 0. rho_stp = 1.2248 -opts_init.dry_distros = {kappa1:lognormal, kappa2:lognormal} +opts_init.dry_distros = {(kappa1, rd_insol):lognormal, (kappa2, rd_insol):lognormal} opts_init.kernel = lgrngn.kernel_t.geometric opts_init.terminal_velocity = lgrngn.vt_t.beard76 opts_init.dt = 1 opts_init.sd_conc = 64 opts_init.n_sd_max = 512 opts_init.rng_seed = 396 -opts_init.src_sd_conc = 64 opts_init.src_z1 = opts_init.dz opts_init.sedi_switch = False diff --git a/tests/python/unit/parcel/test.py b/tests/python/unit/parcel/test.py index b2fc56f19..291633601 100644 --- a/tests/python/unit/parcel/test.py +++ b/tests/python/unit/parcel/test.py @@ -35,7 +35,8 @@ def lognormal(lnr): return n_tot * exp( -pow((lnr - log(mean_r)), 2) / 2 / pow(log(stdev),2) ) / log(stdev) / sqrt(2*pi); -kappa = .61 #TODO - check +kappa = .61 +rd_insol = 0. # rho = 1.8 # TODO chem_gas = { @@ -47,7 +48,7 @@ def lognormal(lnr): # running all three for rhs in [ rhs_blk_2m(dt), - rhs_lgrngn(dt, sd_conc, {kappa:lognormal}, chem_gas) + rhs_lgrngn(dt, sd_conc, {(kappa, rd_insol):lognormal}, chem_gas) ]: parcel(p_d, th_d, r_v, w, dt, nt, rhs) print(p_d, arr_t([th_dry2std(th_d[0], r_v[0])]), r_v) diff --git a/tests/python/unit/relax.py b/tests/python/unit/relax.py index ee0289794..3f27547a4 100644 --- a/tests/python/unit/relax.py +++ b/tests/python/unit/relax.py @@ -103,8 +103,9 @@ def test(opts_init): # test source with dry_distros kappa = .61 +rd_insol = 0. opts_init = lgrngn.opts_init_t() -opts_init.dry_distros = {kappa:lognormal} +opts_init.dry_distros = {(kappa, rd_insol):lognormal} opts_init.rlx_dry_distros = {kappa: [lognormal_rlx, [0,2],[0,opts_init.dz]]} opts_init.sd_conc = 1024 opts_init.rlx_bins = 1024 diff --git a/tests/python/unit/segfault_20150216.py b/tests/python/unit/segfault_20150216.py index 1d60ee6ad..7e0b3b23b 100644 --- a/tests/python/unit/segfault_20150216.py +++ b/tests/python/unit/segfault_20150216.py @@ -24,8 +24,9 @@ def lognormal(lnr): -pow((lnr - log(mean_r)), 2) / 2 / pow(log(stdev),2) ) / log(stdev) / sqrt(2*pi); kappa = .61 +rd_insol = 0. -opts_init.dry_distros = {kappa:lognormal} +opts_init.dry_distros = {(kappa, rd_insol):lognormal} opts_init.sd_conc = 50 opts_init.n_sd_max = 50 diff --git a/tests/python/unit/source.py b/tests/python/unit/source.py index 2b12e5ff5..04f306609 100644 --- a/tests/python/unit/source.py +++ b/tests/python/unit/source.py @@ -29,8 +29,10 @@ def lognormal_src(lnr): -pow((lnr - log(mean_r)), 2) / 2 / pow(log(stdev),2) ) / log(stdev) / sqrt(2*pi); -def test(opts_init): - opts_init.supstp_src = 50 +kappa = .61 +rd_insol = 0.5e-6 + +def test(opts_init, opts): opts_init.rng_seed = int(time()) opts_init.dt = 1 opts_init.nx = 2; @@ -51,14 +53,14 @@ def test(opts_init): opts_init.adve_switch = 0; opts_init.cond_switch = 0; opts_init.sedi_switch = 0; - - opts = lgrngn.opts_t() + opts_init.ice_switch = 0; opts.adve = 0; opts.chem = 0; opts.sedi = 0; opts.coal = 0; opts.cond = 0; + opts.ice_nucl = 0; rhod = arr_t([[ 1., 1. ],[ 1., 1. ]]) th = arr_t([[300., 300. ],[ 300., 300. ]]) @@ -91,19 +93,19 @@ def test(opts_init): return sd_conc, wet_mom0, wet_mom1 -kappa = .61 - # ----------- test source with dry_distros simple ------------------ print(' --- dry_distros simple src ---') opts_init = lgrngn.opts_init_t() -opts_init.dry_distros = {kappa:lognormal} -opts_init.src_dry_distros = {kappa:lognormal_src} +opts = lgrngn.opts_t() +opts_init.dry_distros = {(kappa, rd_insol):lognormal} opts_init.sd_conc = 1024 -opts_init.src_sd_conc = 512 -opts_init.n_sd_max = int((opts_init.sd_conc * 2 + opts_init.src_sd_conc * 2) * 2) # assuming nx=nz=2 -opts_init.src_type = lgrngn.src_t.simple; +src_sd_conc = 512 +supstp_src = 50 +opts.src_dry_distros = {(kappa, rd_insol):(lognormal_src, src_sd_conc, supstp_src)} +opts_init.n_sd_max = int((opts_init.sd_conc * 2 + src_sd_conc * 2) * 2) # assuming nx=nz=2 +opts_init.src_type = lgrngn.src_t.simple -sd_conc, wet_mom0, wet_mom1 = test(opts_init) +sd_conc, wet_mom0, wet_mom1 = test(opts_init, opts) print('diag_sd_conc', sd_conc) if not(sd_conc[0] == 2048 and sd_conc[2] == 2048): @@ -122,14 +124,16 @@ def test(opts_init): # --------------- test source with dry_distros matching ------------------ print(' --- dry_distros matching src ---') opts_init = lgrngn.opts_init_t() -opts_init.dry_distros = {kappa:lognormal} -opts_init.src_dry_distros = {kappa:lognormal_src} +opts = lgrngn.opts_t() +opts_init.dry_distros = {(kappa, rd_insol):lognormal} opts_init.sd_conc = 1024 -opts_init.src_sd_conc = 512 -opts_init.n_sd_max = int((opts_init.sd_conc * 2 + opts_init.src_sd_conc * 2) * 2) # assuming nx=nz=2 -opts_init.src_type = lgrngn.src_t.matching; +src_sd_conc = 512 +supstp_src = 50 +opts.src_dry_distros = {(kappa, rd_insol):(lognormal_src, src_sd_conc, supstp_src)} +opts_init.n_sd_max = int((opts_init.sd_conc * 2 + src_sd_conc * 2) * 2) # assuming nx=nz=2 +opts_init.src_type = lgrngn.src_t.matching -sd_conc, wet_mom0, wet_mom1 = test(opts_init) +sd_conc, wet_mom0, wet_mom1 = test(opts_init, opts) print('diag_sd_conc', sd_conc) if not((sd_conc[0] == 1164 or sd_conc[0] == 1165) and (sd_conc[2] == 1164 or sd_conc[2] == 1165)): @@ -148,13 +152,15 @@ def test(opts_init): # --------- test source with dry_sizes ------------ print(' --- dry_sizes src ---') opts_init = lgrngn.opts_init_t() -opts_init.dry_sizes = {kappa : {1.e-6 : [30., 20], 15.e-6 : [10., 10]}} -opts_init.src_dry_sizes = {kappa : {1.e-6 : [0.3, 10], 15.e-6 : [0.1, 5]}} +opts = lgrngn.opts_t() +supstp_src = 50 +opts_init.dry_sizes = {(kappa, rd_insol) : {1.e-6 : [30., 20], 15.e-6 : [10., 10]}} +opts.src_dry_sizes = {(kappa, rd_insol) : {1.e-6 : [0.3, 10, supstp_src], 15.e-6 : [0.1, 5, supstp_src]}} opts_init.n_sd_max=240 -opts_init.src_type = lgrngn.src_t.simple; # dry sizes works the same for simple and matching (no matching done) +opts_init.src_type = lgrngn.src_t.simple # dry sizes works the same for simple and matching (no matching done) -sd_conc, wet_mom0, wet_mom1 = test(opts_init) +sd_conc, wet_mom0, wet_mom1 = test(opts_init, opts) print('diag_sd_conc', sd_conc) if not((sd_conc[0] == 60) and (sd_conc[2] == 60)): diff --git a/tests/python/unit/sstp_cond.py b/tests/python/unit/sstp_cond.py index 629f79127..583a2ce52 100644 --- a/tests/python/unit/sstp_cond.py +++ b/tests/python/unit/sstp_cond.py @@ -20,9 +20,11 @@ def test(turb_cond): print('turb_cond = ', turb_cond) opts_init = lgrngn.opts_init_t() kappa = .61 - opts_init.dry_distros = {kappa:lognormal} + rd_insol = 0. + opts_init.dry_distros = {(kappa, rd_insol):lognormal} opts_init.coal_switch=0 opts_init.sedi_switch=0 + opts_init.ice_switch=0 opts_init.dt = 1 opts_init.sd_conc = 64 opts_init.n_sd_max = 512 @@ -87,6 +89,8 @@ def test(turb_cond): prtcls.diag_wet_mom(3); wet_post_spin = copy(frombuffer(prtcls.outbuf()).reshape(opts_init.nx, opts_init.nz)) water_post_spin = 1000. * 4./3. * pi * wet_post_spin + rv + print("water post spin: ", water_post_spin) + print("water post init: ", water_post_init) assert allclose(water_post_spin, water_post_init, atol=0, rtol=1e-10) #some discrepancy due to water density #advect SDs @@ -119,7 +123,7 @@ def test(turb_cond): prtcls.diag_all() prtcls.diag_wet_mom(3); wet_post_adve_cond = copy(frombuffer(prtcls.outbuf()).reshape(opts_init.nx, opts_init.nz)) - print(wet_post_adve, wet_post_adve_cond) + #print(wet_post_adve, wet_post_adve_cond) assert allclose(wet_post_adve, wet_post_adve_cond, atol=0, rtol=5e-2) test(False) diff --git a/tests/python/unit/terminal_velocities.py b/tests/python/unit/terminal_velocities.py index 470a64f1b..4b300c917 100644 --- a/tests/python/unit/terminal_velocities.py +++ b/tests/python/unit/terminal_velocities.py @@ -26,8 +26,9 @@ def lognormal(lnr): ) / log(stdev) / sqrt(2*pi); kappa = .61 +rd_insol = 0. -opts_init.dry_distros = {kappa:lognormal} +opts_init.dry_distros = {(kappa, rd_insol):lognormal} opts_init.sd_conc = 50 opts_init.n_sd_max = 50 diff --git a/tests/python/unit/uniform_init.py b/tests/python/unit/uniform_init.py index c053ecd13..6f7c9f0dc 100644 --- a/tests/python/unit/uniform_init.py +++ b/tests/python/unit/uniform_init.py @@ -42,8 +42,9 @@ def expvolumelnr(lnr): rhod = 1. * np.ones((opts_init.nx, opts_init.ny, opts_init.nz)) + .1 * np.mgrid[1:1+opts_init.nx, 1:1+opts_init.ny, 1:1+opts_init.nz][1] # different densities, hence different water content kappa = 1e-6 +rd_insol = 0. -opts_init.dry_distros = {kappa:expvolumelnr} +opts_init.dry_distros = {(kappa, rd_insol):expvolumelnr} opts_init.sd_conc = 64 opts_init.n_sd_max = opts_init.sd_conc * opts_init.nx * opts_init.ny * opts_init.nz