diff --git a/bindings/python/blk_1m.hpp b/bindings/python/blk_1m.hpp index 163fba10..c4dc24d9 100644 --- a/bindings/python/blk_1m.hpp +++ b/bindings/python/blk_1m.hpp @@ -18,6 +18,7 @@ namespace libcloudphxx { namespace b1m = libcloudphxx::blk_1m; namespace blk_1m { +/* template void adj_cellwise( const b1m::opts_t &opts, @@ -71,26 +72,32 @@ namespace libcloudphxx { dt ); } + */ template - void adj_cellwise_nwtrph( + void adj_cellwise( const b1m::opts_t &opts, + const bp_array &rhod, const bp_array &p, bp_array &th, bp_array &rv, bp_array &rc, + bp_array &rr, const typename arr_t::T_numtype &dt ) { arr_t np2bz_th(np2bz(th)), np2bz_rv(np2bz(rv)), - np2bz_rc(np2bz(rc)); - b1m::adj_cellwise_nwtrph( + np2bz_rc(np2bz(rc)), + np2bz_rr(np2bz(rr)); + b1m::adj_cellwise( opts, + np2bz(rhod), np2bz(p), np2bz_th, np2bz_rv, np2bz_rc, + np2bz_rr, dt ); } @@ -116,7 +123,7 @@ namespace libcloudphxx { } template - void rhs_cellwise_nwtrph( + void rhs_cellwise_revap( const b1m::opts_t &opts, bp_array &dot_th, bp_array &dot_rv, @@ -135,7 +142,7 @@ namespace libcloudphxx { np2bz_dot_rr(np2bz(dot_rr)), np2bz_dot_rv(np2bz(dot_rv)), np2bz_dot_th(np2bz(dot_th)); - b1m::rhs_cellwise_nwtrph( + b1m::rhs_cellwise_revap( opts, np2bz_dot_th, np2bz_dot_rv, @@ -153,7 +160,7 @@ namespace libcloudphxx { template - void rhs_cellwise_nwtrph_ice( + void rhs_cellwise_ice( const b1m::opts_t &opts, bp_array &dot_th, bp_array &dot_rv, @@ -178,7 +185,7 @@ namespace libcloudphxx { np2bz_dot_rr(np2bz(dot_rr)), np2bz_dot_rv(np2bz(dot_rv)), np2bz_dot_th(np2bz(dot_th)); - b1m::rhs_cellwise_nwtrph_ice( + b1m::rhs_cellwise_ice( opts, np2bz_dot_th, np2bz_dot_rv, diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index e2e6b1dd..edb76312 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -168,13 +168,14 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("r_c0", &b1m::opts_t::r_c0) .def_readwrite("r_eps", &b1m::opts_t::r_eps) .def_readwrite("nwtrph_iters", &b1m::opts_t::nwtrph_iters) + .def_readwrite("const_p", &b1m::opts_t::const_p) + .def_readwrite("th_dry", &b1m::opts_t::th_dry) + .def_readwrite("adj_nwtrph", &b1m::opts_t::adj_nwtrph) ; bp::def("adj_cellwise", blk_1m::adj_cellwise); - bp::def("adj_cellwise_constp", blk_1m::adj_cellwise_constp); - bp::def("adj_cellwise_nwtrph", blk_1m::adj_cellwise_nwtrph); bp::def("rhs_cellwise", blk_1m::rhs_cellwise); - bp::def("rhs_cellwise_nwtrph", blk_1m::rhs_cellwise_nwtrph); - bp::def("rhs_cellwise_nwtrph_ice", blk_1m::rhs_cellwise_nwtrph_ice); + bp::def("rhs_cellwise_revap", blk_1m::rhs_cellwise_revap); + bp::def("rhs_cellwise_ice", blk_1m::rhs_cellwise_ice); bp::def("rhs_columnwise", blk_1m::rhs_columnwise); // TODO: handle the returned flux bp::def("rhs_columnwise_ice", blk_1m::rhs_columnwise_ice); @@ -197,6 +198,8 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("sedi", &b2m::opts_t::sedi) .def_readwrite("RH_max", &b2m::opts_t::RH_max) .add_property("dry_distros", &blk_2m::get_dd, &blk_2m::set_dd) + .def_readwrite("const_p", &b2m::opts_t::const_p) + .def_readwrite("th_dry", &b2m::opts_t::th_dry) ; bp::def("rhs_cellwise", blk_2m::rhs_cellwise); bp::def("rhs_columnwise", blk_2m::rhs_columnwise); // TODO: handle the returned flux @@ -344,6 +347,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("const_p", &lgr::opts_init_t::const_p) + .def_readwrite("th_dry", &lgr::opts_init_t::th_dry) ; bp::class_/*, boost::noncopyable*/>("particles_proto_t") .add_property("opts_init", &lgrngn::get_oi) diff --git a/include/libcloudph++/blk_1m/adj_cellwise.hpp b/include/libcloudph++/blk_1m/adj_cellwise.hpp index 5b8f22e6..9eeebd58 100644 --- a/include/libcloudph++/blk_1m/adj_cellwise.hpp +++ b/include/libcloudph++/blk_1m/adj_cellwise.hpp @@ -20,71 +20,77 @@ namespace libcloudphxx template class rhs { - private: + private: quantity rhod; - bool const_p; // true if a constant pressure profile is used (e.g. anelastic model), otherwise pressure is deduced from rhod, rv and T + bool const_p, // true if a constant pressure profile is used (e.g. anelastic model) + th_dry; // true if th is dry potential temperature; "standard" potential temperature otherwise public: - quantity p; // total pressure + quantity p; // total pressure void init( - const quantity &_rhod, - const quantity &_p, - const quantity &th, - const quantity &rv, - const bool _const_p - ) - { + const quantity &_rhod, + const quantity &_p, + const quantity &th, + const quantity &rv, + const bool _const_p, + const bool _th_dry + ) + { const_p = _const_p; - rhod = _rhod; + th_dry = _th_dry; + rhod = _rhod; p = _p; - update(th, rv); - } + update(th, rv); + } - quantity r, rs; - quantity T; + quantity r, rs; + quantity T; - private: + private: void update( - const quantity &th, - const quantity &rv - ) - { + const quantity &th, + const quantity &rv + ) + { r = rv; - if(!const_p) + + if(!const_p && th_dry) { - T = common::theta_dry::T(th, rhod); - p = common::theta_dry::p(rhod, rv, T); + T = common::theta_dry::T(th, rhod); + p = common::theta_dry::p(rhod, rv, T); } - else + else if(const_p && !th_dry) { - T = th * common::theta_std::exner(p); + T = th * common::theta_std::exner(p); } + else throw std::runtime_error("adj_cellwise: one (and only one) of opts.const_p and opts.th_dry must be true"); - rs = common::const_cp::r_vs(T, p); - } + rs = common::const_cp::r_vs(T, p); + } - public: + public: - // F = d th / d rv + // F = d th / d rv void operator()( - const quantity &th, - quantity &F, - const quantity &rv - ) - { - update(th, rv); - F = common::theta_dry::d_th_d_rv(T, th); - } + const quantity &th, + quantity &F, + const quantity &rv + ) + { + update(th, rv); + F = common::theta_dry::d_th_d_rv(T, th); + } }; } template void adj_cellwise_nwtrph( const opts_t &opts, + const cont_t &rhod_cont, const cont_t &p_cont, cont_t &th_cont, cont_t &rv_cont, @@ -95,51 +101,70 @@ namespace libcloudphxx using namespace common; if (!opts.cond) return; // ignoring values of opts.cevp - for (auto tup : zip(p_cont, th_cont, rv_cont, rc_cont)) + for (auto tup : zip(rhod_cont, p_cont, th_cont, rv_cont, rc_cont)) { - const quantity - p = std::get<0>(tup) * si::pascals; - real_t - &th = std::get<1>(tup), - &rv = std::get<2>(tup), - &rc = std::get<3>(tup); + &th = std::get<2>(tup), + &rv = std::get<3>(tup), + &rc = std::get<4>(tup); // double-checking.... - //assert(th >= 273.15); // TODO: that's theta, not T! - assert(rc >= 0); - assert(rv >= 0); + //assert(th >= 273.15); // TODO: that's theta, not T! + assert(rc >= 0); + assert(rv >= 0); real_t drc = 0; real_t rv_tmp = rv; quantity th_tmp = th * si::kelvins; - auto exner_p = theta_std::exner(p); + quantity T, T_tmp; + quantity rhod; + quantity p; + quantity exner; - quantity T = th_tmp * exner_p; + if(!opts.const_p && opts.th_dry) + { + rhod = std::get<0>(tup) * si::kilograms / si::cubic_metres; + T = common::theta_dry::T(th_tmp, rhod); + p = common::theta_dry::p(rhod, rv, T); + } + else if(opts.const_p && !opts.th_dry) + { + p = std::get<1>(tup) * si::pascals; + exner = common::theta_std::exner(p); + T = th_tmp * exner; + } + else throw std::runtime_error("adj_cellwise: one (and only one) of opts.const_p and opts.th_dry must be true"); // constant l_v used in theta update auto L0 = const_cp::l_v(T); + T_tmp = T; for (int iter = 0; iter < opts.nwtrph_iters; ++iter) { // TODO: use the approximate Tetens formulas for p_vs and r_vs from tetens.hpp? - quantity p_vs = const_cp::p_vs(T); + quantity p_vs = const_cp::p_vs(T_tmp); // tricky, constant L0 comes from theta = theta + L0 / (c_pd * exner_p) * drc // while variable L comes from dp_vs/dT - auto L = const_cp::l_v(T); - real_t coeff = L * L0 / (moist_air::c_pd() * moist_air::R_v()) / (T * T) / (1 - p_vs / p); + auto L = const_cp::l_v(T_tmp); + real_t coeff = L * L0 / (moist_air::c_pd() * moist_air::R_v()) / (T_tmp * T_tmp) / (1 - p_vs / p); - real_t r_vs = const_cp::r_vs(T, p); + real_t r_vs = const_cp::r_vs(T_tmp, p); drc += (rv_tmp - r_vs) / (1 + coeff * r_vs); rv_tmp = rv - drc; - th_tmp = th * si::kelvins + L0 / (moist_air::c_pd() * exner_p) * drc; + th_tmp = th * si::kelvins + th_tmp / T_tmp * L0 / (moist_air::c_pd()) * drc; - T = th_tmp * exner_p; + if(opts.th_dry) + T_tmp = common::theta_dry::T(th_tmp, rhod); + else + T_tmp = th_tmp * exner; + + if(!opts.const_p) + p = common::theta_dry::p(rhod, rv_tmp, T_tmp); } // limiting @@ -147,27 +172,28 @@ namespace libcloudphxx rv -= drc; rc += drc; - th += L0 / (moist_air::c_pd() * exner_p) * drc / si::kelvins; + th += th / T * L0 / (moist_air::c_pd()) * drc; - // triple-checking.... - //assert(th >= 273.15); // that is theta, not T ! TODO - assert(rc >= 0); - assert(rv >= 0); + // triple-checking.... + //assert(th >= 273.15); // that is theta, not T ! TODO + assert(rc >= 0); + assert(rv >= 0); } } +// Below are saturation adjustment functions that use RK4 integration (from boost.odeint) + // template - void adj_cellwise_hlpr( + void adj_cellwise_rk4( const opts_t &opts, const cont_t &rhod_cont, - const cont_t &p_cont, // value not used if const_p = false - cont_t &th_cont, // if const_p == false, its th_dry, otherwise its th_std + const cont_t &p_cont, // value not used if opts.const_p = false + cont_t &th_cont, // if opts._thdry, its th_dry, otherwise its th_std cont_t &rv_cont, cont_t &rc_cont, cont_t &rr_cont, - const real_t &dt, - const bool const_p + const real_t &dt ) // { @@ -177,152 +203,140 @@ namespace libcloudphxx // odeint::euler< // TODO: opcja? odeint::runge_kutta4< - quantity, // state_type - real_t, // value_type - quantity, // deriv_type - quantity, // time_type - odeint::vector_space_algebra, - odeint::default_operations, - odeint::never_resizer + quantity, // state_type + real_t, // value_type + quantity, // deriv_type + quantity, // time_type + odeint::vector_space_algebra, + odeint::default_operations, + odeint::never_resizer > S; // TODO: would be better to instantiate in the ctor (but what about thread safety! :() typename detail::rhs F; for (auto tup : zip(rhod_cont, p_cont, th_cont, rv_cont, rc_cont, rr_cont)) { const real_t - &rhod = std::get<0>(tup), - &p = std::get<1>(tup); + &rhod = std::get<0>(tup); +// p = opts.const_p ? std::get<1>(tup) : 0; real_t &th = std::get<2>(tup), &rv = std::get<3>(tup), &rc = std::get<4>(tup), - &rr = std::get<5>(tup); + &rr = std::get<5>(tup), + p = 0; + + if(opts.const_p) p = std::get<1>(tup); - // double-checking.... - //assert(th >= 273.15); // TODO: that's theta, not T! - assert(rc >= 0); - assert(rv >= 0); - assert(rr >= 0); + // double-checking.... + //assert(th >= 273.15); // TODO: that's theta, not T! + assert(rc >= 0); + assert(rv >= 0); + assert(rr >= 0); - F.init( - rhod * si::kilograms / si::cubic_metres, + F.init( + rhod * si::kilograms / si::cubic_metres, p * si::pascals, - th * si::kelvins, - rv * si::dimensionless(), - const_p - ); - - real_t vapour_excess; - real_t drr_max = 0; - if (F.rs > F.r && rr > 0 && opts.revp) + th * si::kelvins, + rv * si::dimensionless(), + opts.const_p, + opts.th_dry + ); + + real_t vapour_excess; + real_t drr_max = 0; + if (F.rs > F.r && rr > 0 && opts.revp) { drr_max = (dt * si::seconds) * formulae::evaporation_rate( F.r, F.rs, rr * si::dimensionless(), rhod * si::kilograms / si::cubic_metres, F.p - ); + ); } - bool incloud; + bool incloud; - // TODO: rethink and document r_eps!!! - while ( - // condensation of cloud water if supersaturated more than a threshold - (vapour_excess = rv - F.rs) > opts.r_eps - || + // TODO: rethink and document r_eps!!! + while ( + // condensation of cloud water if supersaturated more than a threshold + (vapour_excess = rv - F.rs) > opts.r_eps + || ( opts.cevp && vapour_excess < -opts.r_eps && ( // or if subsaturated and - (incloud = (rc > 0)) // in cloud (then cloud evaporation first) - || // or + (incloud = (rc > 0)) // in cloud (then cloud evaporation first) + || // or (opts.revp && rr > 0 && drr_max > 0) // in rain shaft (rain evaporation out-of-cloud) - ) + ) ) - ) - { + ) + { // an arbitrary initial guess for drv - real_t drv = - copysign(std::min(.5 * opts.r_eps, .5 * vapour_excess), vapour_excess); + real_t drv = - copysign(std::min(.5 * opts.r_eps, .5 * vapour_excess), vapour_excess); // preventing negative mixing ratios if evaporating - if (vapour_excess < 0) drv = + if (vapour_excess < 0) drv = incloud ? std::min(rc, drv) // limiting by rc - : std::min(drr_max, std::min(rr, drv)); // limiting by rr and drr_max - assert(drv != 0); // otherwise it should not pass the while condition! + : std::min(drr_max, std::min(rr, drv)); // limiting by rr and drr_max + assert(drv != 0); // otherwise it should not pass the while condition! - // theta is modified by do_step, and hence we cannot pass an expression and we need a temp. var. - quantity tmp = th * si::kelvins; + // theta is modified by do_step, and hence we cannot pass an expression and we need a temp. var. + quantity tmp = th * si::kelvins; - // integrating the First Law for moist air - S.do_step( - boost::ref(F), - tmp, - rv * si::dimensionless(), - drv * si::dimensionless() - ); + // integrating the First Law for moist air + S.do_step( + boost::ref(F), + tmp, + rv * si::dimensionless(), + drv * si::dimensionless() + ); - // latent heat source/sink due to evaporation/condensation - th = tmp / si::kelvins; + // latent heat source/sink due to evaporation/condensation + th = tmp / si::kelvins; - // updating rv - rv += drv; - assert(rv >= 0); + // updating rv + rv += drv; + assert(rv >= 0); - if (vapour_excess > 0 || incloud) - { + if (vapour_excess > 0 || incloud) + { // condensation or evaporation of cloud water - rc -= drv; - assert(rc >= 0); - } - else - { + rc -= drv; + assert(rc >= 0); + } + else + { // evaporation of rain water - assert(opts.revp); // should be guaranteed by the while() condition above - rr -= drv; - assert(rr >= 0); - if ((drr_max -= drv) == 0) break; // but not more than Kessler allows - } - } - - // hopefully true for RK4 - assert(F.r == rv); - // triple-checking.... - //assert(th >= 273.15); // that is theta, not T ! TODO - assert(rc >= 0); - assert(rv >= 0); - assert(rr >= 0); + assert(opts.revp); // should be guaranteed by the while() condition above + rr -= drv; + assert(rr >= 0); + if ((drr_max -= drv) == 0) break; // but not more than Kessler allows + } + } + + // hopefully true for RK4 + assert(F.r == rv); + // triple-checking.... + //assert(th >= 273.15); // that is theta, not T ! TODO + assert(rc >= 0); + assert(rv >= 0); + assert(rr >= 0); } } -// saturation adjustment with variable pressure calculated using rhod, rv and T -//- template - void adj_cellwise( - const opts_t &opts, - const cont_t &rhod_cont, - cont_t &th_cont, // th_dry - cont_t &rv_cont, - cont_t &rc_cont, - cont_t &rr_cont, - const real_t &dt - ) -// - { - // rhod_cont passed twice on purpose - it's a dummy var for pressures to make the tuple compile - adj_cellwise_hlpr(opts, rhod_cont, rhod_cont, th_cont, rv_cont, rc_cont, rr_cont, dt, false); - } -// saturation adjustment with a constant pressure profile (e.g. anleastic model) -// needs a different name, because boost python got confused (TODO: fix it) // template - void adj_cellwise_constp( + void adj_cellwise( const opts_t &opts, - const cont_t &rhod_cont, - const cont_t &p_cont, - cont_t &th_cont, // th_std + const cont_t &rhod_cont, // used only if opts.th_dry == true + const cont_t &p_cont, // used only if opts.const_p == true + cont_t &th_cont, cont_t &rv_cont, cont_t &rc_cont, - cont_t &rr_cont, + cont_t &rr_cont, // used only if opts_adj_nwtrph == false const real_t &dt ) // { - adj_cellwise_hlpr(opts, rhod_cont, p_cont, th_cont, rv_cont, rc_cont, rr_cont, dt, true); + if(opts.adj_nwtrph) + adj_cellwise_nwtrph(opts, rhod_cont, p_cont, th_cont, rv_cont, rc_cont, dt); + else + adj_cellwise_rk4(opts, rhod_cont, p_cont, th_cont, rv_cont, rc_cont, rr_cont, dt); } } }; diff --git a/include/libcloudph++/blk_1m/options.hpp b/include/libcloudph++/blk_1m/options.hpp index 55e7f9e1..308b0341 100644 --- a/include/libcloudph++/blk_1m/options.hpp +++ b/include/libcloudph++/blk_1m/options.hpp @@ -35,7 +35,14 @@ namespace libcloudphxx r_c0 = 5e-4, // autoconv. threshold k_acnv = 0.001, // Kessler autoconversion (eq. 5a in Grabowski & Smolarkiewicz 1996) r_eps = 2e-5; // absolute tolerance + + bool adj_nwtrph = true; // if true, use simpler Newton-Raphson iteration in saturation adjustment; otherwise use RK4 from boost.odeint int nwtrph_iters = 3; // number of iterations in Newton-Raphson saturation adjustment + + // NOTE: only tested combinations are: th_dry == true && const_p == false; th_dry == false && const_p == true + // NOTE: th_dry == true && const_p == false doesn't work very well with Newton-Raphson, e.g. sat_adj_blk_1m test (TODO: probably Newton-Raphson needs to be fixed) + bool th_dry = true, // if true, input and output theta are dry-air potential temperature; if false, they are 'standard' potential temperature + const_p = false; // if true, pressure is equal to a supplied profile except for solving velocity (e.g. anelastic model); if false, pressure comes from the gas equation }; // }; diff --git a/include/libcloudph++/blk_1m/rhs_cellwise.hpp b/include/libcloudph++/blk_1m/rhs_cellwise.hpp index 0b1f6812..a92d173f 100644 --- a/include/libcloudph++/blk_1m/rhs_cellwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_cellwise.hpp @@ -14,16 +14,14 @@ namespace libcloudphxx { namespace blk_1m { - // template - void rhs_cellwise( + void rhs_cellwise_hlpr( const opts_t &opts, cont_t &dot_rc_cont, cont_t &dot_rr_cont, const cont_t &rc_cont, const cont_t &rr_cont ) -// { for (auto tup : zip(dot_rc_cont, dot_rr_cont, rc_cont, rr_cont)) { @@ -64,7 +62,20 @@ namespace libcloudphxx } template - void rhs_cellwise_nwtrph( + void rhs_cellwise( + const opts_t &opts, + cont_t &dot_rc_cont, + cont_t &dot_rr_cont, + const cont_t &rc_cont, + const cont_t &rr_cont + ) + { + assert(!opts.adj_nwtrph && "libcloudph++: rhs_cellwise without rain evap requires RK4 in adj_cellwise"); + rhs_cellwise_hlpr(opts, dot_rc_cont, dot_rr_cont, rc_cont, rr_cont); + } + + template + void rhs_cellwise_revap( const opts_t &opts, cont_t &dot_th_cont, cont_t &dot_rv_cont, @@ -79,7 +90,8 @@ namespace libcloudphxx const real_t &dt ) { - rhs_cellwise(opts, dot_rc_cont, dot_rr_cont, rc_cont, rr_cont); + assert(opts.adj_nwtrph && "libcloudph++: rhs_cellwise with rain evap requires Newton-Raphson in adj_cellwise"); + rhs_cellwise_hlpr(opts, dot_rc_cont, dot_rr_cont, rc_cont, rr_cont); // rain evaporation treated as a force in Newthon-Raphson saturation adjustment for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rr_cont, rhod_cont, p_cont, th_cont, rv_cont, rr_cont)) @@ -95,8 +107,8 @@ namespace libcloudphxx const quantity rhod = std::get<3>(tup) * si::kilograms / si::cubic_metres; - const quantity - p = std::get<4>(tup) * si::pascals; + // const quantity + // p = std::get<4>(tup) * si::pascals; const quantity th = std::get<5>(tup) * si::kelvins; @@ -105,7 +117,22 @@ namespace libcloudphxx &rv = std::get<6>(tup), &rr = std::get<7>(tup); - quantity T = th * theta_std::exner(p); + // quantity T = th * theta_std::exner(p); + quantity T; + quantity p; + + if(!opts.const_p && opts.th_dry) + { + T = common::theta_dry::T(th, rhod); + p = common::theta_dry::p(rhod, rv, T); + } + else if(opts.const_p && !opts.th_dry) + { + p = std::get<4>(tup) * si::pascals; + T = th * common::theta_std::exner(p); + } + else throw std::runtime_error("rhs_cellwise: one (and only one) of opts.const_p and opts.th_dry must be true"); + real_t r_vs = const_cp::r_vs(T, p); rr_to_rv += ( @@ -123,12 +150,13 @@ namespace libcloudphxx dot_rv += rr_to_rv; dot_rr -= rr_to_rv; - dot_th -= const_cp::l_v(T) / (moist_air::c_pd() * theta_std::exner(p)) * rr_to_rv / si::kelvins; + //dot_th -= const_cp::l_v(T) / (moist_air::c_pd() * theta_std::exner(p)) * rr_to_rv / si::kelvins; + dot_th += common::theta_dry::d_th_d_rv(T, th) * rr_to_rv / si::kelvins; } } template - void rhs_cellwise_nwtrph_ice( + void rhs_cellwise_ice( const opts_t &opts, cont_t &dot_th_cont, cont_t &dot_rv_cont, @@ -148,8 +176,11 @@ namespace libcloudphxx ) { // autoconversion, collection and rain evaporation: - rhs_cellwise_nwtrph(opts, dot_th_cont, dot_rv_cont, dot_rc_cont, dot_rr_cont, rhod_cont, p_cont, - th_cont, rv_cont, rc_cont, rr_cont, dt); + if(opts.adj_nwtrph) + rhs_cellwise_revap(opts, dot_th_cont, dot_rv_cont, dot_rc_cont, dot_rr_cont, rhod_cont, p_cont, + th_cont, rv_cont, rc_cont, rr_cont, dt); + else + rhs_cellwise(opts, dot_rc_cont, dot_rr_cont, rc_cont, rr_cont); for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rc_cont, dot_rr_cont, dot_ria_cont, dot_rib_cont, rhod_cont, p_cont, th_cont, rv_cont, rc_cont, rr_cont, ria_cont, rib_cont)) @@ -175,8 +206,8 @@ namespace libcloudphxx const quantity rhod = std::get<6>(tup) * si::kilograms / si::cubic_metres; - const quantity - p = std::get<7>(tup) * si::pascals; + // const quantity + // p = std::get<7>(tup) * si::pascals; const quantity th = std::get<8>(tup) * si::kelvins; @@ -188,7 +219,22 @@ namespace libcloudphxx &ria = std::get<12>(tup), &rib = std::get<13>(tup); - quantity T = th * theta_std::exner(p); + // quantity T = th * theta_std::exner(p); + quantity T; + quantity p; + + if(!opts.const_p && opts.th_dry) + { + T = common::theta_dry::T(th, rhod); + p = common::theta_dry::p(rhod, rv, T); + } + else if(opts.const_p && !opts.th_dry) + { + p = std::get<7>(tup) * si::pascals; + T = th * common::theta_std::exner(p); + } + else throw std::runtime_error("rhs_cellwise: one (and only one) of opts.const_p and opts.th_dry must be true"); + real_t rvs = const_cp::r_vs(T, p); real_t rvsi = const_cp::r_vsi(T, p); @@ -369,9 +415,9 @@ namespace libcloudphxx dot_rr += ria_to_rr - rr_to_rib + rib_to_rr; dot_ria += rc_to_ria + rv_to_ria - ria_to_rib - ria_to_rr; dot_rib += rr_to_rib + ria_to_rib + rv_to_rib + rc_to_rib - rib_to_rr; - dot_th += const_cp::l_s(T) / (moist_air::c_pd() * theta_std::exner(p)) * (rv_to_ria + rv_to_rib) / + dot_th += th / T * const_cp::l_s(T) / (moist_air::c_pd()) * (rv_to_ria + rv_to_rib) / si::kelvins; //heat of sublimation - dot_th += const_cp::l_f(T) / (moist_air::c_pd() * theta_std::exner(p)) * (rc_to_ria + rc_to_rib + + dot_th += th / T * const_cp::l_f(T) / (moist_air::c_pd()) * (rc_to_ria + rc_to_rib + rr_to_rib - rib_to_rr - ria_to_rr) / si::kelvins; //heat of freezing } diff --git a/include/libcloudph++/blk_2m/options.hpp b/include/libcloudph++/blk_2m/options.hpp index 0b3f4e64..d8ddaa8f 100644 --- a/include/libcloudph++/blk_2m/options.hpp +++ b/include/libcloudph++/blk_2m/options.hpp @@ -45,6 +45,10 @@ namespace libcloudphxx chem_b; // [1] }; std::vector dry_distros; + + // NOTE: only working combinations are: th_dry == true && const_p == false; th_dry == false && const_p == true + bool th_dry = true, // if true, input and output theta are dry-air potential temperature; if false, they are 'standard' potential temperature + const_p = false; // if true, pressure is equal to a supplied profile except for solving velocity (e.g. anelastic model); if false, pressure comes from the gas equation }; // } diff --git a/include/libcloudph++/blk_2m/rhs_cellwise.hpp b/include/libcloudph++/blk_2m/rhs_cellwise.hpp index 81877fd3..8017897d 100644 --- a/include/libcloudph++/blk_2m/rhs_cellwise.hpp +++ b/include/libcloudph++/blk_2m/rhs_cellwise.hpp @@ -34,7 +34,6 @@ namespace libcloudphxx const cont_t &rr_cont, const cont_t &nr_cont, const real_t &dt, - const bool const_p = false, // is pressure constant (e.g. anelastic approximation on UWLCM)? const cont_t &p_cont = cont_t() // pressure, required if const_p == true ) // @@ -46,8 +45,9 @@ namespace libcloudphxx assert(min(rr_cont) >= 0); assert(min(nc_cont) >= 0); assert(min(nr_cont) >= 0); - assert(!const_p || p_cont.size() == th_cont.size()); - assert(!const_p || min(p_cont) > 0); + assert(!opts.const_p || p_cont.size() == th_cont.size()); + assert(!opts.const_p || min(p_cont) > 0); + assert(opts.const_p != opts.th_dry); // either const_p or th_dry using namespace formulae; using namespace common::moist_air; @@ -96,16 +96,17 @@ namespace libcloudphxx quantity T; quantity p; - if(!const_p) + if(!opts.const_p && opts.th_dry) { T = common::theta_dry::T(th, rhod); p = common::theta_dry::p(rhod, rv, T); } - else + else if(opts.const_p && !opts.th_dry) { p = std::get<13>(tup) * si::pascals; T = th * common::theta_std::exner(p); } + else throw std::runtime_error("rhs_cellwise: one (and only one) of opts.const_p and opts.th_dry must be true"); // rhs only due to rhs_cellwise microphysics functions (needed for limiting) real_t local_dot_rc = 0, @@ -140,7 +141,7 @@ namespace libcloudphxx quantity::type, real_t> tmp = activation_rate(n_ccn, nc / si::kilograms, dt * si::seconds); - local_dot_nc += tmp * si::kilograms * si::seconds; + local_dot_nc += tmp * si::kilograms * si::seconds; local_dot_rc += tmp * ccnmass() * si::seconds; } @@ -198,10 +199,7 @@ namespace libcloudphxx } dot_rv -= (local_dot_rc + local_dot_rr); - if(!const_p) - dot_th -= (local_dot_rc + local_dot_rr) * d_th_d_rv(T, th) / si::kelvins; - else - dot_th += common::const_cp::l_v(T) / (common::moist_air::c_pd() * common::theta_std::exner(p)) * (local_dot_rc + local_dot_rr) / si::kelvins; + dot_th -= (local_dot_rc + local_dot_rr) * d_th_d_rv(T, th) / si::kelvins; dot_rc += local_dot_rc; dot_rr += local_dot_rr; dot_nc += local_dot_nc; diff --git a/include/libcloudph++/common/theta_dry.hpp b/include/libcloudph++/common/theta_dry.hpp index db5da170..cf1d49c8 100644 --- a/include/libcloudph++/common/theta_dry.hpp +++ b/include/libcloudph++/common/theta_dry.hpp @@ -58,7 +58,7 @@ namespace libcloudphxx BOOST_GPU_ENABLED quantity d_th_d_rv( const quantity &T, - const quantity &th // theta dry!!! + const quantity &th // theta dry OR theta std (depends how T was calculated from th) ) { return - th / T * const_cp::l_v(T) / c_pd(); } diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp index 8d4de231..6a8d2232 100644 --- a/include/libcloudph++/lgrngn/opts_init.hpp +++ b/include/libcloudph++/lgrngn/opts_init.hpp @@ -188,6 +188,10 @@ namespace libcloudphxx // relaxation time scale [s] real_t rlx_timescale; + // NOTE: only tested combinations are: th_dry == true && const_p == false; th_dry == false && const_p == true + bool th_dry = true, // if true, input and output theta are dry-air potential temperature; if false, they are 'standard' potential temperature + const_p = false; // if true, pressure is equal to a supplied profile except for solving velocity (e.g. anelastic model); if false, pressure comes from the gas equation + // -- ctors --- // ctor with defaults (C++03 compliant) ... diff --git a/models/kinematic_2D/src/kin_cloud_2d_blk_1m.hpp b/models/kinematic_2D/src/kin_cloud_2d_blk_1m.hpp index 824fe19b..f2005440 100644 --- a/models/kinematic_2D/src/kin_cloud_2d_blk_1m.hpp +++ b/models/kinematic_2D/src/kin_cloud_2d_blk_1m.hpp @@ -29,7 +29,7 @@ class kin_cloud_2d_blk_1m : public kin_cloud_2d_common rhod = (*this->mem->G)(this->ijk); libcloudphxx::blk_1m::adj_cellwise( - opts, rhod, th, rv, rc, rr, this->dt + opts, rhod, typename parent_t::arr_t(), th, rv, rc, rr, this->dt ); this->mem->barrier(); } diff --git a/models/kinematic_2D/src/opts_blk_1m.hpp b/models/kinematic_2D/src/opts_blk_1m.hpp index c489dc42..a4759416 100644 --- a/models/kinematic_2D/src/opts_blk_1m.hpp +++ b/models/kinematic_2D/src/opts_blk_1m.hpp @@ -57,4 +57,8 @@ void setopts_micro( {solver_t::ix::rc, {"rc", "[kg kg-1]"}}, {solver_t::ix::rr, {"rr", "[kg kg-1]"}} }; + + rt_params.cloudph_opts.adj_nwtrph = false; + rt_params.cloudph_opts.th_dry = true; + rt_params.cloudph_opts.const_p = false; } diff --git a/models/kinematic_2D/src/opts_blk_2m.hpp b/models/kinematic_2D/src/opts_blk_2m.hpp index 90b1248a..903e5ae7 100644 --- a/models/kinematic_2D/src/opts_blk_2m.hpp +++ b/models/kinematic_2D/src/opts_blk_2m.hpp @@ -59,6 +59,9 @@ void setopts_micro( .chem_b = setup.chem_b }); + rt_params.cloudph_opts.th_dry = true; + rt_params.cloudph_opts.const_p = false; + // output variables rt_params.outvars = { // : make it common among all three micro? diff --git a/models/kinematic_2D/src/opts_lgrngn.hpp b/models/kinematic_2D/src/opts_lgrngn.hpp index 63e06329..ddae0d1a 100644 --- a/models/kinematic_2D/src/opts_lgrngn.hpp +++ b/models/kinematic_2D/src/opts_lgrngn.hpp @@ -272,6 +272,9 @@ void setopts_micro( else if (backend_str == "serial") rt_params.backend = libcloudphxx::lgrngn::serial; else if (backend_str == "multi_CUDA") rt_params.backend = libcloudphxx::lgrngn::multi_CUDA; + rt_params.cloudph_opts_init.th_dry = true; + rt_params.cloudph_opts_init.const_p = false; + rt_params.async = vm["async"].as(); rt_params.cloudph_opts_init.sd_conc = vm["sd_conc"].as(); diff --git a/models/kinematic_2D/tests/paper_GMD_2015/fig_a/CMakeLists.txt b/models/kinematic_2D/tests/paper_GMD_2015/fig_a/CMakeLists.txt index a8cfa132..ef746c64 100644 --- a/models/kinematic_2D/tests/paper_GMD_2015/fig_a/CMakeLists.txt +++ b/models/kinematic_2D/tests/paper_GMD_2015/fig_a/CMakeLists.txt @@ -81,8 +81,8 @@ add_test(travis_2D_kin_cloud_diff_blk_2m bash -c " h5diff --relative=1e-9 -v2 $dir/timestep0000000000.h5 ${CMAKE_CURRENT_SOURCE_DIR}/refdata/$dir/timestep0000000000.h5 /th && echo 'comparing timestep0000009000.h5' && h5diff --relative=0.02 -v2 $dir/timestep0000009000.h5 ${CMAKE_CURRENT_SOURCE_DIR}/refdata/$dir/timestep0000009000.h5 /rv && - h5diff --delta=12e-6 -v2 $dir/timestep0000009000.h5 ${CMAKE_CURRENT_SOURCE_DIR}/refdata/$dir/timestep0000009000.h5 /rr && - h5diff --delta=4e-6 -v2 $dir/timestep0000009000.h5 ${CMAKE_CURRENT_SOURCE_DIR}/refdata/$dir/timestep0000009000.h5 /rc && + h5diff --delta=12e-6 -v2 $dir/timestep0000009000.h5 ${CMAKE_CURRENT_SOURCE_DIR}/refdata/$dir/timestep0000009000.h5 /rr && + h5diff --delta=4.5e-6 -v2 $dir/timestep0000009000.h5 ${CMAKE_CURRENT_SOURCE_DIR}/refdata/$dir/timestep0000009000.h5 /rc && h5diff --delta=0.4 -v2 $dir/timestep0000009000.h5 ${CMAKE_CURRENT_SOURCE_DIR}/refdata/$dir/timestep0000009000.h5 /th || exit 1; done ") diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index 204ae5ca..deef715d 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -159,9 +159,6 @@ namespace libcloudphxx // is it a pure const_multi run, i.e. no sd_conc bool pure_const_multi; - // is a constant, external pressure profile used? (e.g. anelastic model) - bool const_p; - // is it allowed to do substepping, if not, some memory can be saved bool allow_sstp_cond, allow_sstp_chem; @@ -438,7 +435,7 @@ namespace libcloudphxx { resize_real_vctrs.insert(&tmp_device_real_part3); resize_real_vctrs.insert(&tmp_device_real_part4); - if(const_p) + if(opts_init.const_p) resize_real_vctrs.insert(&tmp_device_real_part5); } diff --git a/src/impl/particles_impl_cond_sstp.ipp b/src/impl/particles_impl_cond_sstp.ipp index 8431de4c..a4f08d89 100644 --- a/src/impl/particles_impl_cond_sstp.ipp +++ b/src/impl/particles_impl_cond_sstp.ipp @@ -98,7 +98,7 @@ namespace libcloudphxx thrust_device::vector &Tp(tmp_device_real_part3); // calc Tp - if(!const_p) // variable pressure + if(opts_init.th_dry) { thrust::transform( sstp_tmp_th.begin(), sstp_tmp_th.end(), // input - first arg @@ -107,7 +107,7 @@ namespace libcloudphxx detail::common__theta_dry__T_rhod() ); } - else // external pressure profile + else // th_std { // T = dry2std(th_d, rv) * exner(p_tot) thrust::transform( @@ -117,13 +117,13 @@ namespace libcloudphxx sstp_tmp_p.begin() // input - third arg )), Tp.begin(), // output - detail::common__theta_dry__T_p() + detail::common__theta_std__T_p() ); } // calculating drop growth in a timestep using backward Euler - if(!const_p) + if(!opts_init.const_p) { // particle-specific pressure iterator, used twice // TODO: store the value somewhere? @@ -166,7 +166,7 @@ namespace libcloudphxx ) ); } - else // const_p + else // opts_init.const_p { if(turb_cond) cond_sstp_hlpr(dt, RH_max, Tp, diff --git a/src/impl/particles_impl_hskpng_Tpr.ipp b/src/impl/particles_impl_hskpng_Tpr.ipp index 6d90e8fb..0ce9e70d 100644 --- a/src/impl/particles_impl_hskpng_Tpr.ipp +++ b/src/impl/particles_impl_hskpng_Tpr.ipp @@ -37,7 +37,7 @@ namespace libcloudphxx }; template - struct common__theta_dry__T_p + struct common__theta_std__T_p { template BOOST_GPU_ENABLED @@ -173,9 +173,9 @@ namespace libcloudphxx template void particles_t::impl::hskpng_Tpr() { - if(!const_p) // variable pressure + if(opts_init.th_dry) // th is th_dry { - // T = common::theta_dry::T(th, rhod); + // T = common::theta_dry::T(th, rhod), so th is assumed to be the dry-air potential temperature thrust::transform( th.begin(), th.end(), // input - first arg rhod.begin(), // input - second arg @@ -183,14 +183,14 @@ namespace libcloudphxx detail::common__theta_dry__T_rhod() ); } - else // external pressure profile + else // th is th_std { - // T = th * exner(p_tot) + // T = th * exner(p_tot), so th is considered to be the "standard" potential temperature thrust::transform( th.begin(), th.end(), // input - first arg thrust::make_zip_iterator(thrust::make_tuple(rv.begin(), p.begin())), // input - second and third args T.begin(), // output - detail::common__theta_dry__T_p() + detail::common__theta_std__T_p() ); } @@ -203,9 +203,9 @@ namespace libcloudphxx > > zip_it_t; - if(!const_p) + if(!opts_init.const_p) { - // p = common::theta_dry::p(rhod, r, T); + // p = common::theta_dry::p(rhod, r, T); works for both th = th_dry and th = th_std thrust::transform( zip_it_t(thrust::make_tuple(rhod.begin(), rv.begin(), T.begin())), // input - begin zip_it_t(thrust::make_tuple(rhod.end(), rv.end(), T.end() )), // input - end diff --git a/src/impl/particles_impl_init_sanity_check.ipp b/src/impl/particles_impl_init_sanity_check.ipp index 91865576..2112b4df 100644 --- a/src/impl/particles_impl_init_sanity_check.ipp +++ b/src/impl/particles_impl_init_sanity_check.ipp @@ -150,6 +150,11 @@ namespace libcloudphxx throw std::runtime_error("libcloudph++: rlx_timescale <= 0"); if(opts_init.rlx_switch && opts_init.chem_switch) throw std::runtime_error("libcloudph++: CCN relaxation does not work with chemistry"); + + if(opts_init.const_p && p.is_null()) + 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"); } }; }; diff --git a/src/impl/particles_impl_reserve_hskpng_npart.ipp b/src/impl/particles_impl_reserve_hskpng_npart.ipp index 27ec4369..6c6530b0 100644 --- a/src/impl/particles_impl_reserve_hskpng_npart.ipp +++ b/src/impl/particles_impl_reserve_hskpng_npart.ipp @@ -67,7 +67,7 @@ namespace libcloudphxx sstp_tmp_rv.reserve(opts_init.n_sd_max); sstp_tmp_th.reserve(opts_init.n_sd_max); sstp_tmp_rh.reserve(opts_init.n_sd_max); - if(const_p) // in const_p pressure is not diagnostic (it's constant) - in per-particle sub-stepping it has to be substepped and we need two vectors to do that + if(opts_init.const_p) // in const_p pressure is not diagnostic (it's constant) - in per-particle sub-stepping it has to be substepped and we need two vectors to do that { sstp_tmp_p.reserve(opts_init.n_sd_max); tmp_device_real_part5.reserve(opts_init.n_sd_max); diff --git a/src/impl/particles_impl_sstp.ipp b/src/impl/particles_impl_sstp.ipp index 28ce6191..e2dee1ea 100644 --- a/src/impl/particles_impl_sstp.ipp +++ b/src/impl/particles_impl_sstp.ipp @@ -20,7 +20,7 @@ namespace libcloudphxx thrust_device::vector *fr[n] = { &rv, &th, &rhod, &p }, *to[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh, &sstp_tmp_p }; - for (int ix = 0; ix < ( (const_p && opts_init.exact_sstp_cond) ? n : n-1); ++ix) // TODO: var_rho + for (int ix = 0; ix < ( (opts_init.const_p && opts_init.exact_sstp_cond) ? n : n-1); ++ix) // TODO: var_rho { if(opts_init.exact_sstp_cond) // per-particle version thrust::copy( @@ -47,7 +47,7 @@ namespace libcloudphxx thrust_device::vector *fr[n] = { &rv, &th, &rhod, &p }, *to[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh, &sstp_tmp_p }; - for (int ix = 0; ix < ( (const_p) ? n : n-1); ++ix) // TODO: var_rho + for (int ix = 0; ix < ( (opts_init.const_p) ? n : n-1); ++ix) // TODO: var_rho { thrust::copy( thrust::make_permutation_iterator(fr[ix]->begin(), ijk.begin()+n_part_old), @@ -119,7 +119,7 @@ namespace libcloudphxx *tmp[n] = { &sstp_tmp_rv, &sstp_tmp_th, &sstp_tmp_rh, &sstp_tmp_p }, *dlt[n] = { &tmp_device_real_part, &tmp_device_real_part1, &tmp_device_real_part2, &tmp_device_real_part5 }; - for (int ix = 0; ix < (const_p ? n : n-1); ++ix) + for (int ix = 0; ix < (opts_init.const_p ? n : n-1); ++ix) { const real_t sstp = sstp_cond; if (step == 0) diff --git a/src/particles_init.ipp b/src/particles_init.ipp index bbbf09e1..5460e87d 100644 --- a/src/particles_init.ipp +++ b/src/particles_init.ipp @@ -14,10 +14,10 @@ namespace libcloudphxx // init template void particles_t::init( - const arrinfo_t th, + const arrinfo_t th, // dry-air (if opts_init.th_dry==true) or (else) 'standard' potential temperature const arrinfo_t rv, const arrinfo_t rhod, - const arrinfo_t p, // pressure profile [in Pascals], needed if pressure perturbations are neglected in condensation (e.g. anelastic model) + const arrinfo_t p, // pressure profile [in Pascals], needed if opts_init.const_p == true // defaults to NULL-NULL pair (variable pressure) const arrinfo_t courant_x, // might be NULL const arrinfo_t courant_y, // might be NULL @@ -31,10 +31,8 @@ namespace libcloudphxx if(pimpl->opts_init.rng_seed_init_switch) pimpl->rng.reseed(pimpl->opts_init.rng_seed_init); - // is a constant pressure profile used? - pimpl->const_p = !p.is_null(); // if pressure comes from a profile, sstp_tmp_p also needs to be copied between distributed memories - if(pimpl->const_p && pimpl->allow_sstp_cond && pimpl->opts_init.exact_sstp_cond) + if(pimpl->opts_init.const_p && pimpl->allow_sstp_cond && pimpl->opts_init.exact_sstp_cond) pimpl->distmem_real_vctrs.insert({&pimpl->sstp_tmp_p, detail::no_initial_value}); // initialising Eulerian-Lagrangian coupling @@ -42,7 +40,7 @@ namespace libcloudphxx pimpl->init_e2l(th, &pimpl->th); pimpl->init_e2l(rv, &pimpl->rv); pimpl->init_e2l(rhod, &pimpl->rhod); - if(pimpl->const_p) + if(pimpl->opts_init.const_p) pimpl->init_e2l(p, &pimpl->p); #if !defined(__NVCC__) diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index 8943571f..5b466a50 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -25,7 +25,18 @@ print("melB =", opts.melB) print("r_c0 =", opts.r_c0) print("r_eps =", opts.r_eps) +print("adj_nwtrph =", opts.adj_nwtrph) +print("th_dry =", opts.th_dry) +print("const_p =", opts.const_p) +test_cases = { \ + "RK4, variable pressure and dry theta" : [0, 1, 0], # adj_nwtrph, th_dry, const_p + "RK4, constant pressure and 'standard' theta" : [0, 0, 1], + "Newton-Raphson, variable pressure and dry theta" : [1, 1, 0], + "Newton-Raphson, constant pressure and 'standard' theta" : [1, 0, 1], +} + +# test saturation adjustment rhod = arr_t([1. ]) p = arr_t([1.e5]) th = arr_t([300.]) @@ -35,63 +46,82 @@ dt = 1 dz = 1 -# sat adjustment with variable pressure -th_old = th.copy() -rv_old = rv.copy() -rc_old = rc.copy() -rr_old = rr.copy() -blk_1m.adj_cellwise(opts, rhod, th, rv, rc, rr, dt) -assert th != th_old # some water should have evaporated -assert rv != rv_old -assert rc != rc_old -assert rr == rr_old +def test_sat_adj(opts, name): + print("Testing saturation adjustment with " + name) + th_new = th.copy() + rv_new = rv.copy() + rc_new = rc.copy() + rr_new = rr.copy() + blk_1m.adj_cellwise(opts, rhod, p, th_new, rv_new, rc_new, rr_new, dt) + assert th != th_new # some water should have evaporated + assert rv != rv_new + assert rc != rc_new + assert rr == rr_new -# sat adjustment with constant pressure -th = arr_t([300.]) -rv = arr_t([0. ]) -rc = arr_t([0.01]) -rr = arr_t([0. ]) -blk_1m.adj_cellwise_constp(opts, rhod, p, th, rv, rc, rr, dt) -assert th != th_old # some water should have evaporated -assert rv != rv_old -assert rc != rc_old -assert rr == rr_old +for name, opt in test_cases.items(): + opts.adj_nwtrph = opt[0] + opts.th_dry = opt[1] + opts.const_p = opt[2] + test_sat_adj(opts, name) -# sat adjustment using Newton-Raphson under constant pressure -th = arr_t([300.]) -rv = arr_t([0. ]) -rc = arr_t([0.01]) -blk_1m.adj_cellwise_nwtrph(opts, p, th, rv, rc, dt) -assert th != th_old # some water should have evaporated -assert rv != rv_old -assert rc != rc_old -assert rr == rr_old +# test RHS cellwise +def test_rhs_cell(opts, name): + print("Testing RHS cellwise with " + name) + dot_th = arr_t([0.]) + dot_rv = arr_t([0.]) + dot_rc = arr_t([0.]) + dot_rr = arr_t([0.]) + if opts.adj_nwtrph: + rr = arr_t([0.01]) + blk_1m.rhs_cellwise_revap(opts, dot_th, dot_rv, dot_rc, dot_rr, rhod, p, th, rv, rc, rr, dt) + assert dot_th != 0 # some rain should have evaporated + assert dot_rv != 0 + else: + rr = arr_t([0.00]) + blk_1m.rhs_cellwise(opts, dot_rc, dot_rr, rc, rr) + assert dot_rc != 0 # some water should have coalesced + assert dot_rr != 0 -dot_rc = arr_t([0.]) -dot_rr = arr_t([0.]) -blk_1m.rhs_cellwise(opts, dot_rc, dot_rr, rc, rr) -assert dot_rc != 0 # some water should have coalesced -assert dot_rr != 0 - -# forcing using Newton-Raphson saturation adjustment include rain evaporation -rr = arr_t([0.01]) -dot_th = arr_t([0.]) -dot_rv = arr_t([0.]) -dot_rc = arr_t([0.]) -dot_rr = arr_t([0.]) - -blk_1m.rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, rhod, p, th, rv, rc, rr, dt) -assert dot_rc != 0 # some water should have coalesced -assert dot_rr != 0 -assert dot_th != 0 # some rain should have evaporated -assert dot_rv != 0 +for name, opt in test_cases.items(): + opts.adj_nwtrph = opt[0] + opts.th_dry = opt[1] + opts.const_p = opt[2] + test_rhs_cell(opts, name) +# test RHS columnwise +print("Testing RHS columnwise") rr = arr_t([0.]) +dot_rr = arr_t([0.]) dot_rr_old = dot_rr.copy() flux = blk_1m.rhs_columnwise(opts, dot_rr, rhod, rr, dz) assert flux == 0 assert dot_rr == dot_rr_old # no rain water -> no precip +#test ice, TODO: separate file/test? +def test_rhs_cell_ice(opts, name): + print("Testing RHS cellwise ice with " + name) + th = arr_t([230.]) #testing ice physics + ria = arr_t([0.1]) + rib = arr_t([0.1]) + dot_th = arr_t([0.]) + dot_rv = arr_t([0.]) + dot_rc = arr_t([0.]) + dot_rr = arr_t([0.]) + dot_rv = arr_t([0.]) + dot_ria = arr_t([0.]) + dot_rib = arr_t([0.]) + blk_1m.rhs_cellwise_ice(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, dot_rib, rhod, p, th, rv, rc, rr, ria, rib, dt) + assert dot_ria != 0 + assert dot_rib != 0 + +for name, opt in test_cases.items(): + opts.adj_nwtrph = opt[0] + opts.th_dry = opt[1] + opts.const_p = opt[2] + test_rhs_cell_ice(opts, name) + +#testing sedimentation of ice +print("Testing RHS columnwise ice") th = arr_t([230.]) #testing ice physics ria = arr_t([0.1]) rib = arr_t([0.1]) @@ -100,12 +130,7 @@ dot_rv = arr_t([0.]) dot_ria = arr_t([0.]) dot_rib = arr_t([0.]) -blk_1m.rhs_cellwise_nwtrph_ice(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, dot_rib, rhod, p, th, rv, rc, rr, ria, rib, dt) -assert dot_ria != 0 -assert dot_rib != 0 - -#testing sedimentation of ice flux_iceA = blk_1m.rhs_columnwise_ice(opts, dot_ria, rhod, ria, dz, ice_t.iceA) flux_iceB = blk_1m.rhs_columnwise_ice(opts, dot_rib, rhod, rib, dz, ice_t.iceB) assert flux_iceA != 0 -assert flux_iceB != 0 \ No newline at end of file +assert flux_iceB != 0 diff --git a/tests/python/unit/blk_1m_ice.py b/tests/python/unit/blk_1m_ice.py index 6cc4e1df..89c5fcfd 100644 --- a/tests/python/unit/blk_1m_ice.py +++ b/tests/python/unit/blk_1m_ice.py @@ -7,6 +7,9 @@ from libcloudphxx import blk_1m opts = blk_1m.opts_t() +opts.adj_nwtrph = True +opts.const_p = True +opts.th_dry = False ice_t = blk_1m.ice_t() rhod = arr_t([1. ]) @@ -27,7 +30,7 @@ dot_rv = arr_t([0.]) dot_ria = arr_t([0.]) dot_rib = arr_t([0.]) - blk_1m.rhs_cellwise_nwtrph_ice(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, dot_rib, rhod, p, th, rv, rc, rr, ria, rib, dt) + blk_1m.rhs_cellwise_ice(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, dot_rib, rhod, p, th, rv, rc, rr, ria, rib, dt) flux_rain = blk_1m.rhs_columnwise(opts, dot_rr, rhod, rr, dz) flux_iceA = blk_1m.rhs_columnwise_ice(opts, dot_ria, rhod, ria, dz, ice_t.iceA) flux_iceB = blk_1m.rhs_columnwise_ice(opts, dot_rib, rhod, rib, dz, ice_t.iceB) @@ -66,4 +69,4 @@ assert rib >= 0 assert np.isclose(ria, 2.7564e-05, rtol=1e-5, atol=1e-8) -assert np.isclose(rib, 3.2808e-06, rtol=1e-5, atol=1e-8) \ No newline at end of file +assert np.isclose(rib, 3.2808e-06, rtol=1e-5, atol=1e-8) diff --git a/tests/python/unit/sat_adj_blk_1m.py b/tests/python/unit/sat_adj_blk_1m.py index 9dbacd38..4cd93243 100644 --- a/tests/python/unit/sat_adj_blk_1m.py +++ b/tests/python/unit/sat_adj_blk_1m.py @@ -1,5 +1,7 @@ import sys sys.path.insert(0, "../../bindings/python/") +from timeit import default_timer as timer + from numpy import array as arr_t # ndarray dtype default to float64, while array's is int64! @@ -17,7 +19,7 @@ def supersaturation(T, p, rv): def initial_state(init_sup_sat): rhod = arr_t([1. ]) - th = arr_t([300.]) + th_dry = arr_t([300.]) if init_sup_sat: rv = arr_t([0.02]) else: @@ -26,75 +28,104 @@ def initial_state(init_sup_sat): rr = arr_t([0. ]) dt = 1 - T = common.T(th[0], rhod[0]) - p = common.p(rhod[0], rv[0], T) - ss = supersaturation(T, p, rv[0]) - print("initial supersaturation", ss) - - return rhod, th, rv, rc, rr, dt - -def test_adj_cellwise(init_sup_sat, r_eps = r_eps_def): - opts.r_eps = r_eps - print("[standard adj_cellwise]") - rhod, th, rv, rc, rr, dt = initial_state(init_sup_sat) - blk_1m.adj_cellwise(opts, rhod, th, rv, rc, rr, dt) - - T = common.T(th[0], rhod[0]) - p = common.p(rhod[0], rv[0], T) - ss = supersaturation(T, p, rv[0]) - print("final supersaturation", ss, th[0], rv[0]) - return ss - -def test_adj_cellwise_constp(init_sup_sat, r_eps = r_eps_def): - opts.r_eps = r_eps - print("[constp adj_cellwise]") - rhod, th, rv, rc, rr, dt = initial_state(init_sup_sat) - - # define pressure consistent with adj_cellwise to compare results - p = arr_t([common.p(rhod[0], rv[0], common.T(th[0], rhod[0]))]) + T = common.T(th_dry[0], rhod[0]) + p = arr_t([common.p(rhod[0], rv[0], T)]) + ss = supersaturation(T, p[0], rv[0]) + print("initial supersaturation", ss, "th_dry", th_dry[0], "th_std", common.th_dry2std(th_dry[0], rv[0]), "rv", rv[0], "T", T, "p", p[0], "rc", rc[0]) + + return rhod, th_dry, rv, rc, rr, p, dt + +def test_adj_cellwise(init_sup_sat, _opts, name, r_eps = r_eps_def, nwtrph_iters = nwtrph_iters_def): + _opts.r_eps = r_eps + _opts.nwtrph_iters = nwtrph_iters + rhod, th_dry, rv, rc, rr, p, dt = initial_state(init_sup_sat) + th_std = arr_t([common.th_dry2std(th_dry[0], rv[0])]) + + # only two combinations allowed + if _opts.th_dry and not _opts.const_p: + start = timer() + blk_1m.adj_cellwise(_opts, rhod, p, th_dry, rv, rc, rr, dt) + end = timer() + T = common.T(th_dry[0], rhod[0]) + p = arr_t([common.p(rhod[0], rv[0], T)]) + th_std = arr_t([common.th_dry2std(th_dry[0], rv[0])]) + elif not _opts.th_dry and _opts.const_p: + start = timer() + blk_1m.adj_cellwise(_opts, rhod, p, th_std, rv, rc, rr, dt) + end = timer() + T = common.exner(p[0]) * th_std[0] + th_dry = arr_t([common.th_std2dry(th_std[0], rv[0])]) - #constp requires th_std input - th_std = arr_t([common.th_dry2std(th[0], rv[0])]) - blk_1m.adj_cellwise_constp(opts, rhod, p, th_std, rv, rc, rr, dt) - - T = common.exner(p[0]) * th_std[0] ss = supersaturation(T, p[0], rv[0]) - print("final supersaturation", ss, th_std[0], rv[0]) + print("final supersaturation", ss, "th_dry", th_dry[0], "th_std", th_std[0], "r_v", rv[0], "T", T, "p", p[0], "rc", rc[0]) + print(f"Execution time: {(end - start) * 1e6:.0e} microseconds") return ss -def test_adj_cellwise_nwtrph(init_sup_sat, nwtrph_iters = nwtrph_iters_def): - opts.nwtrph_iters = nwtrph_iters - print("[nwtrph adj_cellwise]") - rhod, th, rv, rc, rr, dt = initial_state(init_sup_sat) - # define pressure consistent with adj_cellwise to compare results - p = arr_t([common.p(rhod[0], rv[0], common.T(th[0], rhod[0]))]) - - #nwtrph requires th_std input - th_std = arr_t([common.th_dry2std(th[0], rv[0])]) - blk_1m.adj_cellwise_nwtrph(opts, p, th_std, rv, rc, dt) - - T = common.exner(p[0]) * th_std[0] - ss = supersaturation(T, p[0], rv[0]) - print("final supersaturation", ss, th_std[0], rv[0]) - return ss +test_cases = { + "RK4, variable pressure and dry theta" : [0, 1, 0], # adj_nwtrph, th_dry, const_p + "RK4, constant pressure and 'standard' theta" : [0, 0, 1], + "Newton-Raphson, variable pressure and dry theta" : [1, 1, 0], + "Newton-Raphson, constant pressure and 'standard' theta" : [1, 0, 1], +} +# precision eps = { # supersaturation - True : { 'org' : 3e-2, 'constp' : 0.1 , 'nwtrph' : 1e-2, - 'org_prec' : 6e-4, 'constp_prec' : 6e-4, 'nwtrph_prec' : 9e-4 }, - # subsaturation - False : { 'org' : 0.5 , 'constp' : 0.5 , 'nwtrph' : 5e-3, - 'org_prec' : 2e-3, 'constp_prec' : 2e-3, 'nwtrph_prec' : 2e-6 } + True : { + "RK4, variable pressure and dry theta" : 3e-2, + "RK4, constant pressure and 'standard' theta" : 3e-2, + "Newton-Raphson, variable pressure and dry theta" : 3, # !! + "Newton-Raphson, constant pressure and 'standard' theta" : 1, + }, + False : { + "RK4, variable pressure and dry theta" : 0.5, + "RK4, constant pressure and 'standard' theta" : 0.5, + "Newton-Raphson, variable pressure and dry theta" : 0.8, # !! + "Newton-Raphson, constant pressure and 'standard' theta" : 5e-3, + } + } +# increased precision +eps_crank = { + # supersaturation + True : { + "RK4, variable pressure and dry theta" : 6e-4, + "RK4, constant pressure and 'standard' theta" : 6e-4, + "Newton-Raphson, variable pressure and dry theta" : 0.8, # !! + "Newton-Raphson, constant pressure and 'standard' theta" : 9e-4, + }, + False : { + "RK4, variable pressure and dry theta" : 2e-3, + "RK4, constant pressure and 'standard' theta" : 2e-3, + "Newton-Raphson, variable pressure and dry theta" : 0.8, # !! + "Newton-Raphson, constant pressure and 'standard' theta" : 2e-6, + } } -for init_sup_sat in [True, False]: +for name, opt in test_cases.items(): + opts.adj_nwtrph = opt[0] + opts.th_dry = opt[1] + opts.const_p = opt[2] + + for init_sup_sat in [True, False]: # default precision - assert(abs(test_adj_cellwise(init_sup_sat)) < eps[init_sup_sat]['org']) - assert(abs(test_adj_cellwise_constp(init_sup_sat)) < eps[init_sup_sat]['constp']) - assert(abs(test_adj_cellwise_nwtrph(init_sup_sat)) < eps[init_sup_sat]['nwtrph']) - - # cranked up precision - assert(abs(test_adj_cellwise(init_sup_sat, 1e-7)) < eps[init_sup_sat]['org_prec']) - assert(abs(test_adj_cellwise_constp(init_sup_sat, 1e-7)) < eps[init_sup_sat]['constp_prec']) - assert(abs(test_adj_cellwise_nwtrph(init_sup_sat, 5)) < eps[init_sup_sat]['nwtrph_prec']) + if init_sup_sat: + print("[ " + name + " default precision ] initial supersaturation") + else: + print("[ " + name + " default precision ] initial subsaturation") + print("maximum final supersaturation allowed (absolute value): ", eps[init_sup_sat][name]) + assert(abs(test_adj_cellwise(init_sup_sat, opts, name)) < eps[init_sup_sat][name]) + + # cranked up precision + if init_sup_sat: + print("[ " + name + " increased precision ] initial supersaturation") + else: + print("[ " + name + " increased precision ] initial subsaturation") + print("maximum final supersaturation allowed (absolute value): ", eps_crank[init_sup_sat][name]) + if opts.adj_nwtrph: + ss = test_adj_cellwise(init_sup_sat, opts, name, nwtrph_iters = 5) + else: + ss = test_adj_cellwise(init_sup_sat, opts, name, r_eps = 1e-7) + assert(abs(ss) < eps_crank[init_sup_sat][name]) + + print("")