From 3a14e33cdfa12bc0ca64bca29ba14d6247cf7b66 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 23 Apr 2018 17:43:33 +0200 Subject: [PATCH 01/74] use T_e and theta_e when calculating change od theta --- src/bicycles.cpp | 11 ++++++----- src/calc_forces.hpp | 6 +++--- src/cases/CasesCommon.hpp | 2 +- src/cases/DYCOMS98.hpp | 10 ++++------ src/cases/DryThermalGMD2015.hpp | 3 ++- src/cases/LasherTrapp2001.hpp | 6 +++++- src/cases/MoistThermalGrabowskiClark99.hpp | 4 ++-- src/slvr_common.hpp | 2 +- 8 files changed, 24 insertions(+), 20 deletions(-) diff --git a/src/bicycles.cpp b/src/bicycles.cpp index 195cd7614f..f46114eca0 100644 --- a/src/bicycles.cpp +++ b/src/bicycles.cpp @@ -23,11 +23,12 @@ // copy external profiles into rt_parameters // TODO: more elegant way template -void copy_profiles(setup::arr_1D_t &th_e, setup::arr_1D_t &rv_e, setup::arr_1D_t &th_ref, setup::arr_1D_t &pre_ref, setup::arr_1D_t &rhod, setup::arr_1D_t &w_LS, setup::arr_1D_t &hgt_v, setup::arr_1D_t &hgt_s, params_t &p) +void copy_profiles(setup::arr_1D_t &th_e, setup::arr_1D_t &T_e, setup::arr_1D_t &rv_e, setup::arr_1D_t &th_ref, setup::arr_1D_t &pre_ref, setup::arr_1D_t &rhod, setup::arr_1D_t &w_LS, setup::arr_1D_t &hgt_v, setup::arr_1D_t &hgt_s, params_t &p) { p.hgt_fctr_sclr = new setup::arr_1D_t(hgt_s.dataFirst(), hgt_s.shape(), blitz::neverDeleteData); p.hgt_fctr_vctr = new setup::arr_1D_t(hgt_v.dataFirst(), hgt_v.shape(), blitz::neverDeleteData); p.th_e = new setup::arr_1D_t(th_e.dataFirst(), th_e.shape(), blitz::neverDeleteData); + p.T_e = new setup::arr_1D_t(T_e.dataFirst(), T_e.shape(), blitz::neverDeleteData); p.rv_e = new setup::arr_1D_t(rv_e.dataFirst(), rv_e.shape(), blitz::neverDeleteData); p.th_ref = new setup::arr_1D_t(th_ref.dataFirst(), th_ref.shape(), blitz::neverDeleteData); p.pre_ref = new setup::arr_1D_t(pre_ref.dataFirst(), pre_ref.shape(), blitz::neverDeleteData); @@ -102,13 +103,13 @@ void run(int nx, int nz, const user_params_t &user_params) setopts_micro(p, user_params, case_ptr); // reference profiles shared among threads - setup::arr_1D_t th_e(nz), rv_e(nz), th_ref(nz), pre_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr_vctr(nz), hgt_fctr_sclr(nz); + setup::arr_1D_t th_e(nz), T_e(nz), rv_e(nz), th_ref(nz), pre_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr_vctr(nz), hgt_fctr_sclr(nz); // rhod needs to be bigger, cause it divides vertical courant number, TODO: should have a halo both up and down, not only up like now; then it should be interpolated in courant calculation // assign their values - case_ptr->env_prof(th_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, nz, user_params); + case_ptr->env_prof(th_e, T_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, nz, user_params); // pass them to rt_params - copy_profiles(th_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, p); + copy_profiles(th_e, T_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, p); // set outvars p.outvars.insert({solver_t::ix::rv, {"rv", "[kg kg-1]"}}); @@ -217,7 +218,7 @@ void run(int nx, int ny, int nz, const user_params_t &user_params) // rhod needs to be bigger, cause it divides vertical courant number, TODO: should have a halo both up and down, not only up like now; then it should be interpolated in courant calculation // assign their values - case_ptr->env_prof(th_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, nz, user_params); + case_ptr->env_prof(th_e, T_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, nz, user_params); // pass them to rt_params copy_profiles(th_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, p); diff --git a/src/calc_forces.hpp b/src/calc_forces.hpp index c68d8e8b97..e325f695b8 100644 --- a/src/calc_forces.hpp +++ b/src/calc_forces.hpp @@ -77,10 +77,10 @@ void slvr_common::th_src(typename parent_t::arr_t &rv) this->vert_grad_fwd(F, alpha, params.dz); nancheck(alpha(ijk), "sum of th flux"); - // change of theta[K/s] = heating[W/m^3] * theta[K] / T[K] / c_p[J/K/kg] / this->rhod[kg/m^3] - alpha(ijk).reindex(this->zero) *= this->state(ix::th)(ijk).reindex(this->zero) / + // change of theta[K/s] = heating[W/m^3] * theta_e[K] / T_e[K] / c_p[J/K/kg] / this->rhod[kg/m^3] + alpha(ijk).reindex(this->zero) *= (*params.th_e)(this->vert_idx) / calc_c_p()(rv(ijk).reindex(this->zero)) / - calc_T()(this->state(ix::th)(ijk).reindex(this->zero), (*params.rhod)(this->vert_idx)) / + (*params.T_e)(this->vert_idx) / (*params.rhod)(this->vert_idx); nancheck2(alpha(ijk), this->state(ix::th)(ijk), "change of theta"); diff --git a/src/cases/CasesCommon.hpp b/src/cases/CasesCommon.hpp index 0eb057f1b1..c1a195ecb1 100644 --- a/src/cases/CasesCommon.hpp +++ b/src/cases/CasesCommon.hpp @@ -64,7 +64,7 @@ namespace setup virtual void setopts(typename concurr_t::solver_t::rt_params_t ¶ms, int nx, int nz, const user_params_t &user_params) {assert(false);}; virtual void setopts(typename concurr_t::solver_t::rt_params_t ¶ms, int nx, int ny, int nz, const user_params_t &user_params) {assert(false);}; virtual void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, int rng_seed) =0; - virtual void env_prof(arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) =0; + virtual void env_prof(arr_1D_t &th_e, , arr_1D_t &T_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) =0; virtual void update_surf_flux_sens(typename concurr_t::solver_t::arr_sub_t &surf_flux_sens, int timestep, real_t dt) {if(timestep==0) surf_flux_sens = 0.;}; virtual void update_surf_flux_lat(typename concurr_t::solver_t::arr_sub_t &surf_flux_lat, int timestep, real_t dt) {if(timestep==0) surf_flux_lat = 0.;}; diff --git a/src/cases/DYCOMS98.hpp b/src/cases/DYCOMS98.hpp index 214d9f05ec..63e82c36cc 100644 --- a/src/cases/DYCOMS98.hpp +++ b/src/cases/DYCOMS98.hpp @@ -169,7 +169,7 @@ namespace setup // calculate the initial environmental theta and rv profiles // alse set w_LS and hgt_fctrs // like in Wojtek's BabyEulag - void env_prof(arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) + void env_prof(arr_1D_t &th_e, arr_1D_t &T_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) { using libcloudphxx::common::moist_air::R_d_over_c_pd; using libcloudphxx::common::moist_air::c_pd; @@ -179,13 +179,11 @@ namespace setup // pressure profile arr_1D_t pre(nz); - // temperature profile - arr_1D_t T(nz); real_t dz = (Z / si::metres) / (nz-1); r_t rt; // input sounding at z=0, for moist air, no liquid water - T(0) = th_l(0.) / si::kelvins * pow(p_0 / p_1000(), R_d_over_c_pd()); + T_e(0) = th_l(0.) / si::kelvins * pow(p_0 / p_1000(), R_d_over_c_pd()); pre(0) = p_0 / si::pascals; th_e(0) = th_l(0.) / si::kelvins; rv_e(0) = rt(0.); @@ -201,7 +199,7 @@ namespace setup for(int k=1; k() / si::joules * si::kelvins * si::kilograms * T(k-1) * (1 + 0.61 * rv_e(k-1)); // (p / rho) of moist air at k-1 + real_t bottom = R_d() / si::joules * si::kelvins * si::kilograms * T_e(k-1) * (1 + 0.61 * rv_e(k-1)); // (p / rho) of moist air at k-1 real_t rho1 = pre(k-1) / bottom; // rho at k-1 pre(k) = pre(k-1) - rho1 * 9.81 * dz; // estimate of pre at k (dp = -g * rho * dz) real_t thetme = pow(p_1000() / si::pascals / pre(k), f); // 1/Exner @@ -216,7 +214,7 @@ namespace setup if(delta < 0.) delta = 0.; rv_e(k) = rt(k*dz) - delta; th_e(k) = th_l(k*dz) / si::kelvins + c * thetme * delta; - T(k) = th_e(k) * pow(pre(k) / (p_1000() / si::pascals), f); + T_e(k) = th_e(k) * pow(pre(k) / (p_1000() / si::pascals), f); } // compute reference state theta and rhod diff --git a/src/cases/DryThermalGMD2015.hpp b/src/cases/DryThermalGMD2015.hpp index 9a314b3823..4a1aa7a4e7 100644 --- a/src/cases/DryThermalGMD2015.hpp +++ b/src/cases/DryThermalGMD2015.hpp @@ -76,10 +76,11 @@ namespace setup // calculate the initial environmental theta and rv profiles // alse set w_LS and hgt_fctrs // like in Wojtek's BabyEulag - void env_prof(arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) + void env_prof(arr_1D_t &th_e, arr_1D_t &T_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) { rhod = 1; th_e = 300; + T_e = 300; // TODO: not really? th_ref = 300; } }; diff --git a/src/cases/LasherTrapp2001.hpp b/src/cases/LasherTrapp2001.hpp index c45fb9c1d3..56670112ad 100644 --- a/src/cases/LasherTrapp2001.hpp +++ b/src/cases/LasherTrapp2001.hpp @@ -34,6 +34,7 @@ namespace setup // env profiles of th and rv from the sounding arr_1D_t th_dry_env; arr_1D_t th_std_env; + arr_1D_t T_env; arr_1D_t rv_env; @@ -116,7 +117,7 @@ namespace setup // calculate the initial environmental theta and rv profiles // alse set w_LS and hgt_fctrs - void env_prof(arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) + void env_prof(arr_1D_t &th_e, arr_1D_t &T_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) { using libcloudphxx::common::moist_air::R_d_over_c_pd; using libcloudphxx::common::moist_air::c_pd; @@ -179,14 +180,17 @@ namespace setup // create 1D blitz arrays to wrap the derived profiles, store the for use in intcond_hlpr th_dry_env.resize(nz); th_std_env.resize(nz); + T_env.resize(nz); rv_env.resize(nz); th_dry_env = arr_1D_t(th_dry.data(), blitz::shape(nz), blitz::neverDeleteData).copy(); th_std_env = arr_1D_t(th_std.data(), blitz::shape(nz), blitz::neverDeleteData).copy(); + T_env = arr_1D_t(temp_si.data(), blitz::shape(nz), blitz::neverDeleteData).copy(); rv_env = arr_1D_t(rv.data(), blitz::shape(nz), blitz::neverDeleteData).copy(); // TODO: calc hydrostatic env profiles like in dycoms? w kodzie od S. L-T tego jednak nie ma... rv_e = rv_env; + T_e = T_env; th_e = th_std_env; // temp to calc rhod // compute reference state theta and rhod diff --git a/src/cases/MoistThermalGrabowskiClark99.hpp b/src/cases/MoistThermalGrabowskiClark99.hpp index a8517c00bb..f0bc33fd40 100644 --- a/src/cases/MoistThermalGrabowskiClark99.hpp +++ b/src/cases/MoistThermalGrabowskiClark99.hpp @@ -206,7 +206,7 @@ namespace setup public: // calculate the initial environmental theta and rv profiles as Wojtek does it // i.e. for stable virtual standard potential temperature - void env_prof(arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) + void env_prof(arr_1D_t &th_e, arr_1D_t &T_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) // pre_ref - total pressure // th_e - dry potential temp // th_ref - dry potential temp refrence profile @@ -292,7 +292,7 @@ namespace setup T(k)=T(k)/(1.+a*rv_e(k)); } rv_e(k) = RH_T_p_to_rv(env_RH, T(k) * si::kelvins, pre_ref(k) * si::pascals); // cheating! - + T_e(k) = T(k); } //th_ref = th_std_fctr(th_std_0 / si::kelvins)(k * dz); diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index e9dc4d9f4a..3b4d5c551f 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -311,7 +311,7 @@ class slvr_common : public slvr_dim int spinup = 0, // number of timesteps during which autoconversion is to be turned off nt; // total number of timesteps bool rv_src, th_src, uv_src, w_src, subsidence, friction, buoyancy_wet, radiation; - setup::arr_1D_t *th_e, *rv_e, *th_ref, *pre_ref, *rhod, *w_LS, *hgt_fctr_sclr, *hgt_fctr_vctr; + setup::arr_1D_t *th_e, *T_e, *rv_e, *th_ref, *pre_ref, *rhod, *w_LS, *hgt_fctr_sclr, *hgt_fctr_vctr; typename ct_params_t::real_t dz; // vertical grid size setup::ForceParameters_t ForceParameters; user_params_t user_params; // copy od user_params needed only for output to const.h5, since the output has to be done at the end of hook_ante_loop From 1b554eaa7eb7b5f83db398b43a755041436856fe Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 24 Apr 2018 11:41:45 +0200 Subject: [PATCH 02/74] use p_e instead of T_e to make condensation easier --- src/bicycles.cpp | 12 +++++------ src/calc_forces.hpp | 12 ++++++++--- src/cases/CasesCommon.hpp | 2 +- src/cases/DYCOMS98.hpp | 24 +++++++++++----------- src/cases/DryThermalGMD2015.hpp | 9 ++++++-- src/cases/LasherTrapp2001.hpp | 11 +++++----- src/cases/MoistThermalGrabowskiClark99.hpp | 4 ++-- src/slvr_blk_1m.hpp | 2 +- src/slvr_common.hpp | 2 +- 9 files changed, 44 insertions(+), 34 deletions(-) diff --git a/src/bicycles.cpp b/src/bicycles.cpp index f46114eca0..c4cb38376c 100644 --- a/src/bicycles.cpp +++ b/src/bicycles.cpp @@ -23,12 +23,12 @@ // copy external profiles into rt_parameters // TODO: more elegant way template -void copy_profiles(setup::arr_1D_t &th_e, setup::arr_1D_t &T_e, setup::arr_1D_t &rv_e, setup::arr_1D_t &th_ref, setup::arr_1D_t &pre_ref, setup::arr_1D_t &rhod, setup::arr_1D_t &w_LS, setup::arr_1D_t &hgt_v, setup::arr_1D_t &hgt_s, params_t &p) +void copy_profiles(setup::arr_1D_t &th_e, setup::arr_1D_t &p_e, setup::arr_1D_t &rv_e, setup::arr_1D_t &th_ref, setup::arr_1D_t &pre_ref, setup::arr_1D_t &rhod, setup::arr_1D_t &w_LS, setup::arr_1D_t &hgt_v, setup::arr_1D_t &hgt_s, params_t &p) { p.hgt_fctr_sclr = new setup::arr_1D_t(hgt_s.dataFirst(), hgt_s.shape(), blitz::neverDeleteData); p.hgt_fctr_vctr = new setup::arr_1D_t(hgt_v.dataFirst(), hgt_v.shape(), blitz::neverDeleteData); p.th_e = new setup::arr_1D_t(th_e.dataFirst(), th_e.shape(), blitz::neverDeleteData); - p.T_e = new setup::arr_1D_t(T_e.dataFirst(), T_e.shape(), blitz::neverDeleteData); + p.p_e = new setup::arr_1D_t(p_e.dataFirst(), p_e.shape(), blitz::neverDeleteData); p.rv_e = new setup::arr_1D_t(rv_e.dataFirst(), rv_e.shape(), blitz::neverDeleteData); p.th_ref = new setup::arr_1D_t(th_ref.dataFirst(), th_ref.shape(), blitz::neverDeleteData); p.pre_ref = new setup::arr_1D_t(pre_ref.dataFirst(), pre_ref.shape(), blitz::neverDeleteData); @@ -103,13 +103,13 @@ void run(int nx, int nz, const user_params_t &user_params) setopts_micro(p, user_params, case_ptr); // reference profiles shared among threads - setup::arr_1D_t th_e(nz), T_e(nz), rv_e(nz), th_ref(nz), pre_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr_vctr(nz), hgt_fctr_sclr(nz); + setup::arr_1D_t th_e(nz), p_e(nz), rv_e(nz), th_ref(nz), pre_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr_vctr(nz), hgt_fctr_sclr(nz); // rhod needs to be bigger, cause it divides vertical courant number, TODO: should have a halo both up and down, not only up like now; then it should be interpolated in courant calculation // assign their values - case_ptr->env_prof(th_e, T_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, nz, user_params); + case_ptr->env_prof(th_e, p_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, nz, user_params); // pass them to rt_params - copy_profiles(th_e, T_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, p); + copy_profiles(th_e, p_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, p); // set outvars p.outvars.insert({solver_t::ix::rv, {"rv", "[kg kg-1]"}}); @@ -218,7 +218,7 @@ void run(int nx, int ny, int nz, const user_params_t &user_params) // rhod needs to be bigger, cause it divides vertical courant number, TODO: should have a halo both up and down, not only up like now; then it should be interpolated in courant calculation // assign their values - case_ptr->env_prof(th_e, T_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, nz, user_params); + case_ptr->env_prof(th_e, p_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, nz, user_params); // pass them to rt_params copy_profiles(th_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, p); diff --git a/src/calc_forces.hpp b/src/calc_forces.hpp index e325f695b8..38140af5eb 100644 --- a/src/calc_forces.hpp +++ b/src/calc_forces.hpp @@ -17,6 +17,13 @@ struct calc_T BZ_DECLARE_FUNCTOR2(calc_T) }; +struct calc_exner +{ + setup::real_t operator()(setup::real_t p) const + {return libcloudphxx::common::theta_dry::exner(p) / si::pascals;} + BZ_DECLARE_FUNCTOR(calc_exner) +}; + // forcing functions // TODO: make functions return blitz arrays to avoid unnecessary copies template @@ -77,10 +84,9 @@ void slvr_common::th_src(typename parent_t::arr_t &rv) this->vert_grad_fwd(F, alpha, params.dz); nancheck(alpha(ijk), "sum of th flux"); - // change of theta[K/s] = heating[W/m^3] * theta_e[K] / T_e[K] / c_p[J/K/kg] / this->rhod[kg/m^3] - alpha(ijk).reindex(this->zero) *= (*params.th_e)(this->vert_idx) / + // change of theta[K/s] = heating[W/m^3] * exner / c_p[J/K/kg] / this->rhod[kg/m^3] + alpha(ijk).reindex(this->zero) *= calc_exner()((*params.p_e)(this->vert_idx)) / calc_c_p()(rv(ijk).reindex(this->zero)) / - (*params.T_e)(this->vert_idx) / (*params.rhod)(this->vert_idx); nancheck2(alpha(ijk), this->state(ix::th)(ijk), "change of theta"); diff --git a/src/cases/CasesCommon.hpp b/src/cases/CasesCommon.hpp index c1a195ecb1..6534ecf5b9 100644 --- a/src/cases/CasesCommon.hpp +++ b/src/cases/CasesCommon.hpp @@ -64,7 +64,7 @@ namespace setup virtual void setopts(typename concurr_t::solver_t::rt_params_t ¶ms, int nx, int nz, const user_params_t &user_params) {assert(false);}; virtual void setopts(typename concurr_t::solver_t::rt_params_t ¶ms, int nx, int ny, int nz, const user_params_t &user_params) {assert(false);}; virtual void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, int rng_seed) =0; - virtual void env_prof(arr_1D_t &th_e, , arr_1D_t &T_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) =0; + virtual void env_prof(arr_1D_t &th_e, arr_1D_t &p_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) =0; virtual void update_surf_flux_sens(typename concurr_t::solver_t::arr_sub_t &surf_flux_sens, int timestep, real_t dt) {if(timestep==0) surf_flux_sens = 0.;}; virtual void update_surf_flux_lat(typename concurr_t::solver_t::arr_sub_t &surf_flux_lat, int timestep, real_t dt) {if(timestep==0) surf_flux_lat = 0.;}; diff --git a/src/cases/DYCOMS98.hpp b/src/cases/DYCOMS98.hpp index 63e82c36cc..6ddc2b3924 100644 --- a/src/cases/DYCOMS98.hpp +++ b/src/cases/DYCOMS98.hpp @@ -169,7 +169,7 @@ namespace setup // calculate the initial environmental theta and rv profiles // alse set w_LS and hgt_fctrs // like in Wojtek's BabyEulag - void env_prof(arr_1D_t &th_e, arr_1D_t &T_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) + void env_prof(arr_1D_t &th_e, arr_1D_t &p_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) { using libcloudphxx::common::moist_air::R_d_over_c_pd; using libcloudphxx::common::moist_air::c_pd; @@ -177,14 +177,14 @@ namespace setup using libcloudphxx::common::const_cp::l_tri; using libcloudphxx::common::theta_std::p_1000; - // pressure profile - arr_1D_t pre(nz); + // temp profile + arr_1D_t T(nz); real_t dz = (Z / si::metres) / (nz-1); r_t rt; // input sounding at z=0, for moist air, no liquid water - T_e(0) = th_l(0.) / si::kelvins * pow(p_0 / p_1000(), R_d_over_c_pd()); - pre(0) = p_0 / si::pascals; + T(0) = th_l(0.) / si::kelvins * pow(p_0 / p_1000(), R_d_over_c_pd()); + p_e(0) = p_0 / si::pascals; th_e(0) = th_l(0.) / si::kelvins; rv_e(0) = rt(0.); @@ -199,22 +199,22 @@ namespace setup for(int k=1; k() / si::joules * si::kelvins * si::kilograms * T_e(k-1) * (1 + 0.61 * rv_e(k-1)); // (p / rho) of moist air at k-1 - real_t rho1 = pre(k-1) / bottom; // rho at k-1 - pre(k) = pre(k-1) - rho1 * 9.81 * dz; // estimate of pre at k (dp = -g * rho * dz) - real_t thetme = pow(p_1000() / si::pascals / pre(k), f); // 1/Exner + real_t bottom = R_d() / si::joules * si::kelvins * si::kilograms * T(k-1) * (1 + 0.61 * rv_e(k-1)); // (p / rho) of moist air at k-1 + real_t rho1 = p_e(k-1) / bottom; // rho at k-1 + p_e(k) = p_e(k-1) - rho1 * 9.81 * dz; // estimate of pre at k (dp = -g * rho * dz) + real_t thetme = pow(p_1000() / si::pascals / p_e(k), f); // 1/Exner real_t thi = 1. / (th_l(k * dz) / si::kelvins); // 1/theta_std real_t y = b * thetme * tt0 * thi; real_t ees = ee0 * exp(b-y); // saturation vapor pressure (Tetens equation or what?) - real_t qvs = a * ees / (pre(k) - ees); // saturation vapor mixing ratio = R_d / R_v * ees / p_d + real_t qvs = a * ees / (p_e(k) - ees); // saturation vapor mixing ratio = R_d / R_v * ees / p_d // calculate linearized condensation rate real_t cf1 = thetme*thetme*thi*thi; // T^{-2} - cf1 *= c * d * pre(k) / (pre(k) - ees); // = l_tri^2 / (C_pd * R_v * T^2) * p/p_d + cf1 *= c * d * p_e(k) / (p_e(k) - ees); // = l_tri^2 / (C_pd * R_v * T^2) * p/p_d real_t delta = (rt(k*dz) - qvs) / (1 + qvs * cf1); // how much supersaturated is the air (divided by sth) if(delta < 0.) delta = 0.; rv_e(k) = rt(k*dz) - delta; th_e(k) = th_l(k*dz) / si::kelvins + c * thetme * delta; - T_e(k) = th_e(k) * pow(pre(k) / (p_1000() / si::pascals), f); + T(k) = th_e(k) * pow(p_e(k) / (p_1000() / si::pascals), f); } // compute reference state theta and rhod diff --git a/src/cases/DryThermalGMD2015.hpp b/src/cases/DryThermalGMD2015.hpp index 4a1aa7a4e7..509f76c2cf 100644 --- a/src/cases/DryThermalGMD2015.hpp +++ b/src/cases/DryThermalGMD2015.hpp @@ -76,11 +76,16 @@ namespace setup // calculate the initial environmental theta and rv profiles // alse set w_LS and hgt_fctrs // like in Wojtek's BabyEulag - void env_prof(arr_1D_t &th_e, arr_1D_t &T_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) + void env_prof(arr_1D_t &th_e, arr_1D_t &p_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) { + using libcloudphxx::common::theta_std::p_1000 + using libcloudphxx::common::moist_air::R_d; + using libcloudphxx::common::moist_air::c_pd; + using libcloudphxx::common::moist_air::R_d_over_c_pd; + rhod = 1; th_e = 300; - T_e = 300; // TODO: not really? + p_e = pow((rhod * si::kilograms / si::cubic_metres) * R_d() * (th_e * si::kelvins) * pow(p_1000(), R_d_over_c_pd()), 1 / (R_d_over_c_pd() + 1)) / si::pascals; th_ref = 300; } }; diff --git a/src/cases/LasherTrapp2001.hpp b/src/cases/LasherTrapp2001.hpp index 56670112ad..c9e956f775 100644 --- a/src/cases/LasherTrapp2001.hpp +++ b/src/cases/LasherTrapp2001.hpp @@ -34,7 +34,7 @@ namespace setup // env profiles of th and rv from the sounding arr_1D_t th_dry_env; arr_1D_t th_std_env; - arr_1D_t T_env; + arr_1D_t p_env; arr_1D_t rv_env; @@ -117,7 +117,7 @@ namespace setup // calculate the initial environmental theta and rv profiles // alse set w_LS and hgt_fctrs - void env_prof(arr_1D_t &th_e, arr_1D_t &T_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) + void env_prof(arr_1D_t &th_e, arr_1D_t &p_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) { using libcloudphxx::common::moist_air::R_d_over_c_pd; using libcloudphxx::common::moist_air::c_pd; @@ -180,17 +180,16 @@ namespace setup // create 1D blitz arrays to wrap the derived profiles, store the for use in intcond_hlpr th_dry_env.resize(nz); th_std_env.resize(nz); - T_env.resize(nz); + p_env.resize(nz); rv_env.resize(nz); th_dry_env = arr_1D_t(th_dry.data(), blitz::shape(nz), blitz::neverDeleteData).copy(); th_std_env = arr_1D_t(th_std.data(), blitz::shape(nz), blitz::neverDeleteData).copy(); - T_env = arr_1D_t(temp_si.data(), blitz::shape(nz), blitz::neverDeleteData).copy(); + p_env = arr_1D_t(pres_si.data(), blitz::shape(nz), blitz::neverDeleteData).copy(); rv_env = arr_1D_t(rv.data(), blitz::shape(nz), blitz::neverDeleteData).copy(); // TODO: calc hydrostatic env profiles like in dycoms? w kodzie od S. L-T tego jednak nie ma... - + p_e = p_env; rv_e = rv_env; - T_e = T_env; th_e = th_std_env; // temp to calc rhod // compute reference state theta and rhod diff --git a/src/cases/MoistThermalGrabowskiClark99.hpp b/src/cases/MoistThermalGrabowskiClark99.hpp index f0bc33fd40..7162d6d95a 100644 --- a/src/cases/MoistThermalGrabowskiClark99.hpp +++ b/src/cases/MoistThermalGrabowskiClark99.hpp @@ -206,7 +206,7 @@ namespace setup public: // calculate the initial environmental theta and rv profiles as Wojtek does it // i.e. for stable virtual standard potential temperature - void env_prof(arr_1D_t &th_e, arr_1D_t &T_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) + void env_prof(arr_1D_t &th_e, arr_1D_t &p_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) // pre_ref - total pressure // th_e - dry potential temp // th_ref - dry potential temp refrence profile @@ -292,7 +292,7 @@ namespace setup T(k)=T(k)/(1.+a*rv_e(k)); } rv_e(k) = RH_T_p_to_rv(env_RH, T(k) * si::kelvins, pre_ref(k) * si::pascals); // cheating! - T_e(k) = T(k); + p_e(k) = pre_ref(k); } //th_ref = th_std_fctr(th_std_0 / si::kelvins)(k * dz); diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index c2b272c8d0..8f54c6c31b 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -28,7 +28,7 @@ class slvr_blk_1m_common : public slvr_common rhod = (*this->mem->G)(this->ijk); libcloudphxx::blk_1m::adj_cellwise( - opts, rhod, th, rv, rc, rr, this->dt + opts, rhod, p_e, th, rv, rc, rr, this->dt ); this->mem->barrier(); } diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index 3b4d5c551f..12e74195a0 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -311,7 +311,7 @@ class slvr_common : public slvr_dim int spinup = 0, // number of timesteps during which autoconversion is to be turned off nt; // total number of timesteps bool rv_src, th_src, uv_src, w_src, subsidence, friction, buoyancy_wet, radiation; - setup::arr_1D_t *th_e, *T_e, *rv_e, *th_ref, *pre_ref, *rhod, *w_LS, *hgt_fctr_sclr, *hgt_fctr_vctr; + setup::arr_1D_t *th_e, *p_e, *rv_e, *th_ref, *pre_ref, *rhod, *w_LS, *hgt_fctr_sclr, *hgt_fctr_vctr; typename ct_params_t::real_t dz; // vertical grid size setup::ForceParameters_t ForceParameters; user_params_t user_params; // copy od user_params needed only for output to const.h5, since the output has to be done at the end of hook_ante_loop From c9711f9c9b5763fc34031c89d5198041869f77e4 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 24 Apr 2018 12:42:55 +0200 Subject: [PATCH 03/74] a 2D/3D arr for p_e in blk_1m --- src/bicycles.cpp | 4 ++-- src/calc_forces.hpp | 2 +- src/cases/DryThermalGMD2015.hpp | 8 ++++++-- src/slvr_blk_1m.hpp | 11 ++++++++++- src/slvr_lgrngn.hpp | 1 + 5 files changed, 20 insertions(+), 6 deletions(-) diff --git a/src/bicycles.cpp b/src/bicycles.cpp index c4cb38376c..0daebe3144 100644 --- a/src/bicycles.cpp +++ b/src/bicycles.cpp @@ -214,13 +214,13 @@ void run(int nx, int ny, int nz, const user_params_t &user_params) setopts_micro(p, user_params, case_ptr); // reference profiles shared among threads - setup::arr_1D_t th_e(nz), rv_e(nz), th_ref(nz), pre_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr_vctr(nz), hgt_fctr_sclr(nz); + setup::arr_1D_t th_e(nz), p_e(nz), rv_e(nz), th_ref(nz), pre_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr_vctr(nz), hgt_fctr_sclr(nz); // rhod needs to be bigger, cause it divides vertical courant number, TODO: should have a halo both up and down, not only up like now; then it should be interpolated in courant calculation // assign their values case_ptr->env_prof(th_e, p_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, nz, user_params); // pass them to rt_params - copy_profiles(th_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, p); + copy_profiles(th_e, p_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, p); // set outvars p.outvars.insert({solver_t::ix::rv, {"rv", "[kg kg-1]"}}); diff --git a/src/calc_forces.hpp b/src/calc_forces.hpp index 38140af5eb..78d9222a17 100644 --- a/src/calc_forces.hpp +++ b/src/calc_forces.hpp @@ -20,7 +20,7 @@ struct calc_T struct calc_exner { setup::real_t operator()(setup::real_t p) const - {return libcloudphxx::common::theta_dry::exner(p) / si::pascals;} + {return libcloudphxx::common::theta_dry::exner(p * si::pascals);} BZ_DECLARE_FUNCTOR(calc_exner) }; diff --git a/src/cases/DryThermalGMD2015.hpp b/src/cases/DryThermalGMD2015.hpp index 509f76c2cf..f28cf33ddd 100644 --- a/src/cases/DryThermalGMD2015.hpp +++ b/src/cases/DryThermalGMD2015.hpp @@ -78,14 +78,18 @@ namespace setup // like in Wojtek's BabyEulag void env_prof(arr_1D_t &th_e, arr_1D_t &p_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) { - using libcloudphxx::common::theta_std::p_1000 + using libcloudphxx::common::theta_std::p_1000; using libcloudphxx::common::moist_air::R_d; using libcloudphxx::common::moist_air::c_pd; using libcloudphxx::common::moist_air::R_d_over_c_pd; rhod = 1; th_e = 300; - p_e = pow((rhod * si::kilograms / si::cubic_metres) * R_d() * (th_e * si::kelvins) * pow(p_1000(), R_d_over_c_pd()), 1 / (R_d_over_c_pd() + 1)) / si::pascals; + //p_e = pow((rhod * si::kilograms / si::cubic_metres) * R_d() * (th_e * si::kelvins) * pow(p_1000(), R_d_over_c_pd()), 1. / (R_d_over_c_pd() + 1)) / si::pascals; +// const quantity p_theta = (quantity(rhod * si::kilograms / si::cubic_metres) * R_d() * quantity(th_e * si::kelvins))/si::pascals; + const quantity p_theta = (quantity(1. * si::kilograms / si::cubic_metres) * R_d() * quantity(300. * si::kelvins)); + /*const quantity*/ real_t p = pow( real_t(p_theta / si::pascals) * pow(real_t(p_1000() / si::pascals), R_d_over_c_pd()), 1. / (R_d_over_c_pd() + 1)) ;// * pow(p_1000(), R_d_over_c_pd()));//, 1. / (R_d_over_c_pd() + 1)) / si::pascals; + p_e = p; th_ref = 300; } }; diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index 8f54c6c31b..29071dceed 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -17,6 +17,10 @@ class slvr_blk_1m_common : public slvr_common using arr_sub_t = typename parent_t::arr_sub_t; private: + // a 2D/3D array with a copy of the environmental pressure profile, + // it's needed by adj_cellwise_constp + blitz::Array p_e; + void condevap() { auto @@ -27,7 +31,7 @@ class slvr_blk_1m_common : public slvr_common auto const rhod = (*this->mem->G)(this->ijk); - libcloudphxx::blk_1m::adj_cellwise( + libcloudphxx::blk_1m::adj_cellwise_constp( opts, rhod, p_e, th, rv, rc, rr, this->dt ); this->mem->barrier(); @@ -56,6 +60,11 @@ class slvr_blk_1m_common : public slvr_common parent_t::hook_ante_loop(nt); // forcings after adjustments + // init the p_e array + p_e.resize(this->shape(this->ijk)); + p_e = (*this->params.p_e)(this->vert_idx); + p_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? + // recording parameters if(this->rank==0) { diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index baa1bd1a66..afa078567d 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -503,6 +503,7 @@ class slvr_lgrngn : public slvr_common private: // per-thread copy of params + // TODO: but slvr_common also has a copy of it's params.... rt_params_t params; public: From 3d527118924bd02ddae704bd6cd83a3597ef13b9 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Tue, 24 Apr 2018 14:07:25 +0200 Subject: [PATCH 04/74] sllkvr blk_1m: init p_e b4 first call to condevap --- src/slvr_blk_1m.hpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index 29071dceed..cbd8879022 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -30,7 +30,7 @@ class slvr_blk_1m_common : public slvr_common rr = this->state(ix::rr)(this->ijk); // rain water mixing ratio auto const rhod = (*this->mem->G)(this->ijk); - + libcloudphxx::blk_1m::adj_cellwise_constp( opts, rhod, p_e, th, rv, rc, rr, this->dt ); @@ -55,16 +55,16 @@ class slvr_blk_1m_common : public slvr_common zero_if_uninitialised(ix::rc); zero_if_uninitialised(ix::rr); - // deal with initial supersaturation - condevap(); - - parent_t::hook_ante_loop(nt); // forcings after adjustments - // init the p_e array p_e.resize(this->shape(this->ijk)); p_e = (*this->params.p_e)(this->vert_idx); p_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? + // deal with initial supersaturation + condevap(); + + parent_t::hook_ante_loop(nt); // forcings after adjustments + // recording parameters if(this->rank==0) { From 52151c0f5b0b94d44736d7fc21a177c6f334bd61 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Tue, 24 Apr 2018 15:54:34 +0200 Subject: [PATCH 05/74] blk_1m: also pass p_d_env to adj_cellwise --- src/cases/DryThermalGMD2015.hpp | 2 +- src/slvr_blk_1m.hpp | 18 ++++++++++++++++-- 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/src/cases/DryThermalGMD2015.hpp b/src/cases/DryThermalGMD2015.hpp index f28cf33ddd..6073a51743 100644 --- a/src/cases/DryThermalGMD2015.hpp +++ b/src/cases/DryThermalGMD2015.hpp @@ -89,7 +89,7 @@ namespace setup // const quantity p_theta = (quantity(rhod * si::kilograms / si::cubic_metres) * R_d() * quantity(th_e * si::kelvins))/si::pascals; const quantity p_theta = (quantity(1. * si::kilograms / si::cubic_metres) * R_d() * quantity(300. * si::kelvins)); /*const quantity*/ real_t p = pow( real_t(p_theta / si::pascals) * pow(real_t(p_1000() / si::pascals), R_d_over_c_pd()), 1. / (R_d_over_c_pd() + 1)) ;// * pow(p_1000(), R_d_over_c_pd()));//, 1. / (R_d_over_c_pd() + 1)) / si::pascals; - p_e = p; + p_e = p; // total env pressure th_ref = 300; } }; diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index cbd8879022..8f8238272d 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -6,6 +6,13 @@ #include #include +struct calc_p_v +{ + setup::real_t operator()(setup::real_t p, setup::real_t rv) const + {return libcloudphxx::common::moist_air::p_v(p * si::pascals, rv) / si::pascals;} + BZ_DECLARE_FUNCTOR2(calc_p_v) +}; + template class slvr_blk_1m_common : public slvr_common { @@ -17,9 +24,10 @@ class slvr_blk_1m_common : public slvr_common using arr_sub_t = typename parent_t::arr_sub_t; private: - // a 2D/3D array with a copy of the environmental pressure profile, + // a 2D/3D arrays with copies of the environmental total pressure/partial pressure of dry air, // it's needed by adj_cellwise_constp blitz::Array p_e; + blitz::Array p_d_e; void condevap() { @@ -32,7 +40,7 @@ class slvr_blk_1m_common : public slvr_common rhod = (*this->mem->G)(this->ijk); libcloudphxx::blk_1m::adj_cellwise_constp( - opts, rhod, p_e, th, rv, rc, rr, this->dt + opts, rhod, p_e, p_d_e, th, rv, rc, rr, this->dt ); this->mem->barrier(); } @@ -60,6 +68,12 @@ class slvr_blk_1m_common : public slvr_common p_e = (*this->params.p_e)(this->vert_idx); p_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? + // init the p_d_e array + p_d_e.resize(this->shape(this->ijk)); + // p_d_e = p_e - p_v_e + p_d_e = (*this->params.p_e)(this->vert_idx) - calc_p_v()((*this->params.p_e)(this->vert_idx), (*this->params.rv_e)(this->vert_idx)); + p_d_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? + // deal with initial supersaturation condevap(); From 8f0898396cfb1321de87c9ab8d693515b8bc8162 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 25 Apr 2018 12:09:23 +0200 Subject: [PATCH 06/74] blk_1m: condensation pre-advection, post-half-rhs; r_l stored in update_rhs, i.e. used post-advection, post-half-rhs --- src/slvr_blk_1m.hpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index 8f8238272d..de60d8f10b 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -98,8 +98,9 @@ class slvr_blk_1m_common : public slvr_common void hook_ante_step() { parent_t::hook_ante_step(); + condevap(); // treat saturation adjustment as pre-advection, post-half-rhs adjustment // store rl for buoyancy - this->r_l(this->ijk) = this->state(ix::rc)(this->ijk) + this->state(ix::rr)(this->ijk); + //this->r_l(this->ijk) = this->state(ix::rc)(this->ijk) + this->state(ix::rr)(this->ijk); } void update_rhs( @@ -107,6 +108,9 @@ class slvr_blk_1m_common : public slvr_common const typename parent_t::real_t &dt, const int &at ) { + // store rl for buoyancy + this->r_l(this->ijk) = this->state(ix::rc)(this->ijk) + this->state(ix::rr)(this->ijk); + parent_t::update_rhs(rhs, dt, at); // shouldnt forcings be after condensation to be consistent with lgrngn solver? // cell-wise @@ -124,7 +128,7 @@ class slvr_blk_1m_common : public slvr_common // void hook_post_step() { - condevap(); // treat saturation adjustment as post-advection, pre-rhs adjustment + //condevap(); // treat saturation adjustment as post-advection, pre-rhs adjustment parent_t::hook_post_step(); // includes the above forcings } From 4918726506e8b5432c95cbe492ec91ac0a036651 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 25 Apr 2018 12:26:14 +0200 Subject: [PATCH 07/74] add subsidence of rc and rr in blk_1m --- src/bicycles.cpp | 9 ++- src/calc_forces.hpp | 39 +++++++++++++ src/cases/CasesCommon.hpp | 2 +- src/cases/DYCOMS98.hpp | 2 + src/cases/DryThermalGMD2015.hpp | 2 + src/cases/LasherTrapp2001.hpp | 2 + src/cases/MoistThermalGrabowskiClark99.hpp | 2 + src/slvr_blk_1m.hpp | 64 ++++++++++++++++++++++ 8 files changed, 119 insertions(+), 3 deletions(-) diff --git a/src/bicycles.cpp b/src/bicycles.cpp index 0daebe3144..d9dea33770 100644 --- a/src/bicycles.cpp +++ b/src/bicycles.cpp @@ -103,7 +103,7 @@ void run(int nx, int nz, const user_params_t &user_params) setopts_micro(p, user_params, case_ptr); // reference profiles shared among threads - setup::arr_1D_t th_e(nz), p_e(nz), rv_e(nz), th_ref(nz), pre_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr_vctr(nz), hgt_fctr_sclr(nz); + setup::arr_1D_t th_e(nz), p_e(nz), rv_e(nz), th_ref(nz), pre_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr_vctr(nz), hgt_fctr_sclr(nz); // rhod needs to be bigger, cause it divides vertical courant number, TODO: should have a halo both up and down, not only up like now; then it should be interpolated in courant calculation // assign their values @@ -214,7 +214,7 @@ void run(int nx, int ny, int nz, const user_params_t &user_params) setopts_micro(p, user_params, case_ptr); // reference profiles shared among threads - setup::arr_1D_t th_e(nz), p_e(nz), rv_e(nz), th_ref(nz), pre_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr_vctr(nz), hgt_fctr_sclr(nz); + setup::arr_1D_t th_e(nz), p_e(nz), rv_e(nz), th_ref(nz), pre_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr_vctr(nz), hgt_fctr_sclr(nz); // rhod needs to be bigger, cause it divides vertical courant number, TODO: should have a halo both up and down, not only up like now; then it should be interpolated in courant calculation // assign their values @@ -388,6 +388,7 @@ void run_hlpr(bool piggy, std::string type, Args&&... args) struct ct_params_final : ct_params_piggy { enum { opts = opts::nug | opts::iga | opts::fct }; }; run>(args...); } + } else // piggybacking { @@ -432,6 +433,8 @@ int main(int argc, char** argv) ("serial", po::value()->default_value(false), "force advection and microphysics to be computed on single thread") ("th_src", po::value()->default_value(true) , "temp src") ("rv_src", po::value()->default_value(true) , "water vap source") + ("rc_src", po::value()->default_value(true) , "cloud water source (in blk_1m)") + ("rr_src", po::value()->default_value(true) , "rain water source (in blk_1m)") ("uv_src", po::value()->default_value(true) , "horizontal vel src") ("w_src", po::value()->default_value(true) , "vertical vel src") ("piggy", po::value()->default_value(false) , "is it a piggybacking run") @@ -490,6 +493,8 @@ int main(int argc, char** argv) // handling sources flags user_params.th_src = vm["th_src"].as(); user_params.rv_src = vm["rv_src"].as(); + user_params.rc_src = vm["rc_src"].as(); + user_params.rr_src = vm["rr_src"].as(); user_params.uv_src = vm["uv_src"].as(); user_params.w_src = vm["w_src"].as(); diff --git a/src/calc_forces.hpp b/src/calc_forces.hpp index 78d9222a17..0dbcfdc463 100644 --- a/src/calc_forces.hpp +++ b/src/calc_forces.hpp @@ -133,4 +133,43 @@ void slvr_common::w_src(typename parent_t::arr_t &th, typename pare alpha(ijk) += F(ijk); } +template +void slvr_blk_1m::rc_src() +{ + const auto &ijk = this->ijk; + if(params.rc_src) + { + // large-scale vertical wind + subsidence(ix::rc); + + alpha(ijk) += F(ijk); + } + else + alpha(ijk) = 0.; + + beta(ijk) = 0.; + // nudging, todo: use some other coeff than vab_coeff +// alpha(ijk).reindex(this->zero) += (*this->mem->vab_coeff)(ijk).reindex(this->zero) * (*params.rv_e)(this->vert_idx); // TODO: its a constant, cache it +// beta(ijk) = - (*this->mem->vab_coeff)(ijk); +} + +template +void slvr_blk_1m::rr_src() +{ + const auto &ijk = this->ijk; + if(params.rr_src) + { + // large-scale vertical wind + subsidence(ix::rr); + + alpha(ijk) += F(ijk); + } + else + alpha(ijk) = 0.; + + beta(ijk) = 0.; + // nudging, todo: use some other coeff than vab_coeff +// alpha(ijk).reindex(this->zero) += (*this->mem->vab_coeff)(ijk).reindex(this->zero) * (*params.rv_e)(this->vert_idx); // TODO: its a constant, cache it +// beta(ijk) = - (*this->mem->vab_coeff)(ijk); +} diff --git a/src/cases/CasesCommon.hpp b/src/cases/CasesCommon.hpp index 6534ecf5b9..2ed57978c5 100644 --- a/src/cases/CasesCommon.hpp +++ b/src/cases/CasesCommon.hpp @@ -17,7 +17,7 @@ struct user_params_t int nt, outfreq, spinup, rng_seed; setup::real_t dt, z_rlx_sclr; std::string outdir, model_case; - bool th_src, rv_src, uv_src, w_src; + bool th_src, rv_src, rc_src, rr_src, uv_src, w_src; }; namespace setup diff --git a/src/cases/DYCOMS98.hpp b/src/cases/DYCOMS98.hpp index 6ddc2b3924..cf0ba31fc6 100644 --- a/src/cases/DYCOMS98.hpp +++ b/src/cases/DYCOMS98.hpp @@ -123,6 +123,8 @@ namespace setup params.uv_src = user_params.uv_src; params.th_src = user_params.th_src; params.rv_src = user_params.rv_src; + params.rc_src = user_params.rc_src; + params.rr_src = user_params.rr_src; params.dt = user_params.dt; params.nt = user_params.nt; params.buoyancy_wet = true; diff --git a/src/cases/DryThermalGMD2015.hpp b/src/cases/DryThermalGMD2015.hpp index 6073a51743..6fbaf7d71f 100644 --- a/src/cases/DryThermalGMD2015.hpp +++ b/src/cases/DryThermalGMD2015.hpp @@ -30,6 +30,8 @@ namespace setup params.uv_src = false; params.th_src = false; params.rv_src = false; + params.rc_src = false; + params.rr_src = false; params.dt = user_params.dt; params.nt = user_params.nt; params.buoyancy_wet = false; diff --git a/src/cases/LasherTrapp2001.hpp b/src/cases/LasherTrapp2001.hpp index c9e956f775..38c1d841c7 100644 --- a/src/cases/LasherTrapp2001.hpp +++ b/src/cases/LasherTrapp2001.hpp @@ -60,6 +60,8 @@ namespace setup params.uv_src = false; // ? params.th_src = user_params.th_src; params.rv_src = user_params.rv_src; + params.rc_src = user_params.rc_src; + params.rr_src = user_params.rr_src; params.dt = user_params.dt; params.nt = user_params.nt; params.buoyancy_wet = true; diff --git a/src/cases/MoistThermalGrabowskiClark99.hpp b/src/cases/MoistThermalGrabowskiClark99.hpp index 7162d6d95a..313a942342 100644 --- a/src/cases/MoistThermalGrabowskiClark99.hpp +++ b/src/cases/MoistThermalGrabowskiClark99.hpp @@ -166,6 +166,8 @@ namespace setup params.uv_src = false; params.th_src = false; params.rv_src = false; + params.rc_src = false; + params.rr_src = false; params.dt = user_params.dt; params.nt = user_params.nt; params.buoyancy_wet = true; diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index de60d8f10b..c2ecbcfdeb 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -53,9 +53,68 @@ class slvr_blk_1m_common : public slvr_common protected: + void rc_src(); + void rr_src(); bool get_rain() { return opts.conv; } void set_rain(bool val) { opts.conv = val; }; + + void update_rhs( + arrvec_t &rhs, + const typename parent_t::real_t &dt, + const int &at + ) + { + parent_t::update_rhs(rhs, dt, at); + this->mem->barrier(); + if(this->rank == 0) + this->tbeg = clock::now(); + + using ix = typename ct_params_t::ix; + + const auto &ijk = this->ijk; + + // forcing + switch (at) + { + // for eulerian integration or used to init trapezoidal integration + case (0): + { + // ---- cloud water sources ---- + rc_src(); + rhs.at(ix::rc)(ijk) += alpha(ijk) + beta(ijk) * this->state(ix::rc)(ijk); + + // ---- rain water sources ---- + rr_src(); + rhs.at(ix::rr)(ijk) += alpha(ijk) + beta(ijk) * this->state(ix::rr)(ijk); + + break; + } + + case (1): + // trapezoidal rhs^n+1 + { + // ---- cloud water sources ---- + rc_src(); + rhs.at(ix::rc)(ijk) += alpha(ijk) + beta(ijk) * this->state(ix::rc)(ijk) / (1. - 0.5 * this->dt * beta(ijk)); + + // ---- rain water sources ---- + rr_src(); + rhs.at(ix::rr)(ijk) += alpha(ijk) + beta(ijk) * this->state(ix::rr)(ijk) / (1. - 0.5 * this->dt * beta(ijk)); + + break; + } + } + this->mem->barrier(); + if(this->rank == 0) + { + nancheck(rhs.at(ix::rc)(this->domain), "RHS of rc after rhs_update"); + nancheck(rhs.at(ix::rr)(this->domain), "RHS of rr after rhs_update"); + this->tend = clock::now(); + this->tupdate += std::chrono::duration_cast( this->tend - this->tbeg ); + } + } + // deals with initial supersaturation void hook_ante_loop(int nt) { @@ -92,6 +151,10 @@ class slvr_blk_1m_common : public slvr_common this->record_aux_const("r_c0", opts.r_c0); this->record_aux_const("k_acnv", opts.k_acnv); this->record_aux_const("r_eps", opts.r_eps); + this->record_aux_const("user_params rc_src", params.user_params.rc_src); + this->record_aux_const("user_params rr_src", params.user_params.rr_src); + this->record_aux_const("rt_params rc_src", params.rc_src); + this->record_aux_const("rt_params rr_src", params.rr_src); } } @@ -139,6 +202,7 @@ class slvr_blk_1m_common : public slvr_common struct rt_params_t : parent_t::rt_params_t { libcloudphxx::blk_1m::opts_t cloudph_opts; + bool rr_src, rc_src; }; // ctor From 11d0314e750bbd621f18789e969ec9eb77b39555 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 25 Apr 2018 14:15:55 +0200 Subject: [PATCH 08/74] fix suubsidence of rc and rr in blk_1m --- src/calc_forces_blk_1m.hpp | 45 +++++++++++++ src/calc_forces_common.hpp | 134 +++++++++++++++++++++++++++++++++++++ 2 files changed, 179 insertions(+) create mode 100644 src/calc_forces_blk_1m.hpp create mode 100644 src/calc_forces_common.hpp diff --git a/src/calc_forces_blk_1m.hpp b/src/calc_forces_blk_1m.hpp new file mode 100644 index 0000000000..0174ff2620 --- /dev/null +++ b/src/calc_forces_blk_1m.hpp @@ -0,0 +1,45 @@ +//TODO: move calc_forces and forcings to case class? +#pragma once +#include "forcings.hpp" + +// single-moment bulk forcing functions +template +void slvr_blk_1m_common::rc_src() +{ + const auto &ijk = this->ijk; + if(params.rc_src) + { + // large-scale vertical wind + parent_t::subsidence(ix::rc); + + this->alpha(ijk) += this->F(ijk); + } + else + this->alpha(ijk) = 0.; + + this->beta(ijk) = 0.; + // nudging, todo: use some other coeff than vab_coeff +// this->alpha(ijk).reindex(this->zero) += (*this->mem->vab_coeff)(ijk).reindex(this->zero) * (*params.rv_e)(this->vert_idx); // TODO: its a constant, cache it +// this->beta(ijk) = - (*this->mem->vab_coeff)(ijk); +} + + +template +void slvr_blk_1m_common::rr_src() +{ + const auto &ijk = this->ijk; + if(params.rr_src) + { + // large-scale vertical wind + parent_t::subsidence(ix::rr); + + this->alpha(ijk) += this->F(ijk); + } + else + this->alpha(ijk) = 0.; + + this->beta(ijk) = 0.; + // nudging, todo: use some other coeff than vab_coeff +// this->alpha(ijk).reindex(this->zero) += (*this->mem->vab_coeff)(ijk).reindex(this->zero) * (*params.rv_e)(this->vert_idx); // TODO: its a constant, cache it +// this->beta(ijk) = - (*this->mem->vab_coeff)(ijk); +} diff --git a/src/calc_forces_common.hpp b/src/calc_forces_common.hpp new file mode 100644 index 0000000000..e80a0d9263 --- /dev/null +++ b/src/calc_forces_common.hpp @@ -0,0 +1,134 @@ +//TODO: move calc_forces and forcings to case class? +#pragma once +#include "forcings.hpp" + +// helper functors +struct calc_c_p +{ + setup::real_t operator()(setup::real_t rv) const + {return libcloudphxx::common::moist_air::c_p(rv) * si::kilograms * si::kelvins / si::joules;} + BZ_DECLARE_FUNCTOR(calc_c_p) +}; + +struct calc_T +{ + setup::real_t operator()(setup::real_t th, setup::real_t rhod) const + {return libcloudphxx::common::theta_dry::T(th * si::kelvins, rhod * si::kilograms / si::metres / si::metres / si::metres) / si::kelvins;} + BZ_DECLARE_FUNCTOR2(calc_T) +}; + +struct calc_exner +{ + setup::real_t operator()(setup::real_t p) const + {return libcloudphxx::common::theta_dry::exner(p * si::pascals);} + BZ_DECLARE_FUNCTOR(calc_exner) +}; + +// common forcing functions +// TODO: make functions return blitz arrays to avoid unnecessary copies +template +void slvr_common::rv_src() +{ + const auto &ijk = this->ijk; + if(params.rv_src) + { + // surface flux + surf_latent(); + // sum of rv flux + this->vert_grad_fwd(F, alpha, params.dz); + + // change of rv[1/s] = latent heating[W/m^3] / lat_heat_of_evap[J/kg] / density[kg/m^3] + if(params.ForceParameters.surf_latent_flux_in_watts_per_square_meter) + alpha(ijk).reindex(this->zero) /= (libcloudphxx::common::const_cp::l_tri() * si::kilograms / si::joules) * (*params.rhod)(this->vert_idx); + + // large-scale vertical wind + subsidence(ix::rv); // TODO: in case 1, rv here should be in step n+1, calc it explicitly as rv + 0.5 * dt * rhs(rv); + // could also be calculated implicitly, but we would need implicit rv^n+1 in other cells + + alpha(ijk) += F(ijk); + } + else + alpha(ijk) = 0.; + + beta(ijk) = 0.; + // nudging, todo: use some other coeff than vab_coeff +// alpha(ijk).reindex(this->zero) += (*this->mem->vab_coeff)(ijk).reindex(this->zero) * (*params.rv_e)(this->vert_idx); // TODO: its a constant, cache it +// beta(ijk) = - (*this->mem->vab_coeff)(ijk); +} + +template +void slvr_common::th_src(typename parent_t::arr_t &rv) +{ + const auto &ijk = this->ijk; + if(params.th_src) + { + // -- heating -- + // surface flux + if(params.ForceParameters.surf_sensible_flux_in_watts_per_square_meter) + { + surf_sens(); + nancheck(F(ijk), "sensible surf forcing"); + // beta as tmp storage + beta(ijk) = F(ijk); + } + else + beta(ijk) = 0.; + + // radiation + radiation(rv); + nancheck(beta(ijk), "radiation"); + // add fluxes from radiation and surface + F(ijk) += beta(ijk); + nancheck(F(ijk), "sensible surf forcing + radiation"); + // sum of th flux, F(j) is upward flux through the bottom of the j-th cell + this->vert_grad_fwd(F, alpha, params.dz); + nancheck(alpha(ijk), "sum of th flux"); + + // change of theta[K/s] = heating[W/m^3] * exner / c_p[J/K/kg] / this->rhod[kg/m^3] + alpha(ijk).reindex(this->zero) *= calc_exner()((*params.p_e)(this->vert_idx)) / + calc_c_p()(rv(ijk).reindex(this->zero)) / + (*params.rhod)(this->vert_idx); + + nancheck2(alpha(ijk), this->state(ix::th)(ijk), "change of theta"); + + // surf flux if it is specified as mean(theta*w) + if(!params.ForceParameters.surf_sensible_flux_in_watts_per_square_meter) + { + surf_sens(); + nancheck(F(ijk), "sensible surf forcing"); + this->vert_grad_fwd(F, beta, params.dz); + nancheck(beta(ijk), "grad of sensible surf forcing"); + alpha(ijk) += beta(ijk); + } + + // large-scale vertical wind + subsidence(ix::th); // TODO: in case 1 th here should be in step n+1, calc it explicitly as th + 0.5 * dt * rhs(th); + // could also be calculated implicitly, but we would need implicit th^n+1 in other cells + nancheck(F(ijk), "subsidence"); + + alpha(ijk) += F(ijk); + nancheck(alpha(ijk), "alpha in th_src"); + } + else + alpha(ijk) = 0.; + + beta(ijk) = 0.; + // nudging, todo: use some other coeff than vab_coeff + //alpha(ijk).reindex(this->zero) += (*this->mem->vab_coeff)(ijk).reindex(this->zero) * (*params.th_e)(this->vert_idx); + //beta(ijk) = - (*this->mem->vab_coeff)(ijk); +} + +template +void slvr_common::w_src(typename parent_t::arr_t &th, typename parent_t::arr_t &rv) +{ + const auto &ijk = this->ijk; + // buoyancy + buoyancy(th, rv); + alpha(ijk) = F(ijk); + // large-scale vertical wind + subsidence(ix::w); // TODO: in case 1, w here should be in step n+1, calc it explicitly as w + 0.5 * dt * rhs(w); + // could also be calculated implicitly, but we would need implicit w^n+1 in other cells; + // also include absorber in w^n+1 estimate... + + alpha(ijk) += F(ijk); +} From 7c5002beeb5007599a0560d53aff53cdeaee1577 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 25 Apr 2018 14:19:44 +0200 Subject: [PATCH 09/74] add subsidence of rc and rr in blk_1m cd --- src/bicycles.cpp | 1 + src/calc_forces.hpp | 175 -------------------------------------------- src/opts_blk_1m.hpp | 2 + src/opts_lgrngn.hpp | 2 +- src/slvr_blk_1m.hpp | 145 +++++++++++++++++++++--------------- src/slvr_common.hpp | 4 + 6 files changed, 92 insertions(+), 237 deletions(-) delete mode 100644 src/calc_forces.hpp diff --git a/src/bicycles.cpp b/src/bicycles.cpp index d9dea33770..d9a8c05139 100644 --- a/src/bicycles.cpp +++ b/src/bicycles.cpp @@ -17,6 +17,7 @@ #include "opts_lgrngn.hpp" #include "opts_blk_1m.hpp" + #include "panic.hpp" #include diff --git a/src/calc_forces.hpp b/src/calc_forces.hpp deleted file mode 100644 index 0dbcfdc463..0000000000 --- a/src/calc_forces.hpp +++ /dev/null @@ -1,175 +0,0 @@ -//TODO: move calc_forces and forcings to case class? -#pragma once -#include "forcings.hpp" - -// helper functors -struct calc_c_p -{ - setup::real_t operator()(setup::real_t rv) const - {return libcloudphxx::common::moist_air::c_p(rv) * si::kilograms * si::kelvins / si::joules;} - BZ_DECLARE_FUNCTOR(calc_c_p) -}; - -struct calc_T -{ - setup::real_t operator()(setup::real_t th, setup::real_t rhod) const - {return libcloudphxx::common::theta_dry::T(th * si::kelvins, rhod * si::kilograms / si::metres / si::metres / si::metres) / si::kelvins;} - BZ_DECLARE_FUNCTOR2(calc_T) -}; - -struct calc_exner -{ - setup::real_t operator()(setup::real_t p) const - {return libcloudphxx::common::theta_dry::exner(p * si::pascals);} - BZ_DECLARE_FUNCTOR(calc_exner) -}; - -// forcing functions -// TODO: make functions return blitz arrays to avoid unnecessary copies -template -void slvr_common::rv_src() -{ - const auto &ijk = this->ijk; - if(params.rv_src) - { - // surface flux - surf_latent(); - // sum of rv flux - this->vert_grad_fwd(F, alpha, params.dz); - - // change of rv[1/s] = latent heating[W/m^3] / lat_heat_of_evap[J/kg] / density[kg/m^3] - if(params.ForceParameters.surf_latent_flux_in_watts_per_square_meter) - alpha(ijk).reindex(this->zero) /= (libcloudphxx::common::const_cp::l_tri() * si::kilograms / si::joules) * (*params.rhod)(this->vert_idx); - - // large-scale vertical wind - subsidence(ix::rv); // TODO: in case 1, rv here should be in step n+1, calc it explicitly as rv + 0.5 * dt * rhs(rv); - // could also be calculated implicitly, but we would need implicit rv^n+1 in other cells - - alpha(ijk) += F(ijk); - } - else - alpha(ijk) = 0.; - - beta(ijk) = 0.; - // nudging, todo: use some other coeff than vab_coeff -// alpha(ijk).reindex(this->zero) += (*this->mem->vab_coeff)(ijk).reindex(this->zero) * (*params.rv_e)(this->vert_idx); // TODO: its a constant, cache it -// beta(ijk) = - (*this->mem->vab_coeff)(ijk); -} - -template -void slvr_common::th_src(typename parent_t::arr_t &rv) -{ - const auto &ijk = this->ijk; - if(params.th_src) - { - // -- heating -- - // surface flux - if(params.ForceParameters.surf_sensible_flux_in_watts_per_square_meter) - { - surf_sens(); - nancheck(F(ijk), "sensible surf forcing"); - // beta as tmp storage - beta(ijk) = F(ijk); - } - else - beta(ijk) = 0.; - - // radiation - radiation(rv); - nancheck(beta(ijk), "radiation"); - // add fluxes from radiation and surface - F(ijk) += beta(ijk); - nancheck(F(ijk), "sensible surf forcing + radiation"); - // sum of th flux, F(j) is upward flux through the bottom of the j-th cell - this->vert_grad_fwd(F, alpha, params.dz); - nancheck(alpha(ijk), "sum of th flux"); - - // change of theta[K/s] = heating[W/m^3] * exner / c_p[J/K/kg] / this->rhod[kg/m^3] - alpha(ijk).reindex(this->zero) *= calc_exner()((*params.p_e)(this->vert_idx)) / - calc_c_p()(rv(ijk).reindex(this->zero)) / - (*params.rhod)(this->vert_idx); - - nancheck2(alpha(ijk), this->state(ix::th)(ijk), "change of theta"); - - // surf flux if it is specified as mean(theta*w) - if(!params.ForceParameters.surf_sensible_flux_in_watts_per_square_meter) - { - surf_sens(); - nancheck(F(ijk), "sensible surf forcing"); - this->vert_grad_fwd(F, beta, params.dz); - nancheck(beta(ijk), "grad of sensible surf forcing"); - alpha(ijk) += beta(ijk); - } - - // large-scale vertical wind - subsidence(ix::th); // TODO: in case 1 th here should be in step n+1, calc it explicitly as th + 0.5 * dt * rhs(th); - // could also be calculated implicitly, but we would need implicit th^n+1 in other cells - nancheck(F(ijk), "subsidence"); - - alpha(ijk) += F(ijk); - nancheck(alpha(ijk), "alpha in th_src"); - } - else - alpha(ijk) = 0.; - - beta(ijk) = 0.; - // nudging, todo: use some other coeff than vab_coeff - //alpha(ijk).reindex(this->zero) += (*this->mem->vab_coeff)(ijk).reindex(this->zero) * (*params.th_e)(this->vert_idx); - //beta(ijk) = - (*this->mem->vab_coeff)(ijk); -} - -template -void slvr_common::w_src(typename parent_t::arr_t &th, typename parent_t::arr_t &rv) -{ - const auto &ijk = this->ijk; - // buoyancy - buoyancy(th, rv); - alpha(ijk) = F(ijk); - // large-scale vertical wind - subsidence(ix::w); // TODO: in case 1, w here should be in step n+1, calc it explicitly as w + 0.5 * dt * rhs(w); - // could also be calculated implicitly, but we would need implicit w^n+1 in other cells; - // also include absorber in w^n+1 estimate... - - alpha(ijk) += F(ijk); -} - -template -void slvr_blk_1m::rc_src() -{ - const auto &ijk = this->ijk; - if(params.rc_src) - { - // large-scale vertical wind - subsidence(ix::rc); - - alpha(ijk) += F(ijk); - } - else - alpha(ijk) = 0.; - - beta(ijk) = 0.; - // nudging, todo: use some other coeff than vab_coeff -// alpha(ijk).reindex(this->zero) += (*this->mem->vab_coeff)(ijk).reindex(this->zero) * (*params.rv_e)(this->vert_idx); // TODO: its a constant, cache it -// beta(ijk) = - (*this->mem->vab_coeff)(ijk); -} - - -template -void slvr_blk_1m::rr_src() -{ - const auto &ijk = this->ijk; - if(params.rr_src) - { - // large-scale vertical wind - subsidence(ix::rr); - - alpha(ijk) += F(ijk); - } - else - alpha(ijk) = 0.; - - beta(ijk) = 0.; - // nudging, todo: use some other coeff than vab_coeff -// alpha(ijk).reindex(this->zero) += (*this->mem->vab_coeff)(ijk).reindex(this->zero) * (*params.rv_e)(this->vert_idx); // TODO: its a constant, cache it -// beta(ijk) = - (*this->mem->vab_coeff)(ijk); -} diff --git a/src/opts_blk_1m.hpp b/src/opts_blk_1m.hpp index 0adcfe12e2..4de8e61fbe 100644 --- a/src/opts_blk_1m.hpp +++ b/src/opts_blk_1m.hpp @@ -9,6 +9,8 @@ #include "opts_common.hpp" #include "slvr_blk_1m.hpp" +#include "calc_forces_common.hpp" +#include "calc_forces_blk_1m.hpp" // simulation and output parameters for micro=blk_1m template diff --git a/src/opts_lgrngn.hpp b/src/opts_lgrngn.hpp index ab5ff21070..a16f020c5a 100644 --- a/src/opts_lgrngn.hpp +++ b/src/opts_lgrngn.hpp @@ -11,7 +11,7 @@ #include "opts_common.hpp" #include "slvr_lgrngn.hpp" -#include "calc_forces.hpp" +#include "calc_forces_common.hpp" // string parsing #include diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index c2ecbcfdeb..9d802645ce 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -22,6 +22,7 @@ class slvr_blk_1m_common : public slvr_common using ix = typename ct_params_t::ix; // TODO: it's now in solver_common - is it needed here? using real_t = typename ct_params_t::real_t; using arr_sub_t = typename parent_t::arr_sub_t; + using clock = typename parent_t::clock; private: // a 2D/3D arrays with copies of the environmental total pressure/partial pressure of dry air, @@ -59,62 +60,6 @@ class slvr_blk_1m_common : public slvr_common void set_rain(bool val) { opts.conv = val; }; - void update_rhs( - arrvec_t &rhs, - const typename parent_t::real_t &dt, - const int &at - ) - { - parent_t::update_rhs(rhs, dt, at); - this->mem->barrier(); - if(this->rank == 0) - this->tbeg = clock::now(); - - using ix = typename ct_params_t::ix; - - const auto &ijk = this->ijk; - - // forcing - switch (at) - { - // for eulerian integration or used to init trapezoidal integration - case (0): - { - // ---- cloud water sources ---- - rc_src(); - rhs.at(ix::rc)(ijk) += alpha(ijk) + beta(ijk) * this->state(ix::rc)(ijk); - - // ---- rain water sources ---- - rr_src(); - rhs.at(ix::rr)(ijk) += alpha(ijk) + beta(ijk) * this->state(ix::rr)(ijk); - - break; - } - - case (1): - // trapezoidal rhs^n+1 - { - // ---- cloud water sources ---- - rc_src(); - rhs.at(ix::rc)(ijk) += alpha(ijk) + beta(ijk) * this->state(ix::rc)(ijk) / (1. - 0.5 * this->dt * beta(ijk)); - - // ---- rain water sources ---- - rr_src(); - rhs.at(ix::rr)(ijk) += alpha(ijk) + beta(ijk) * this->state(ix::rr)(ijk) / (1. - 0.5 * this->dt * beta(ijk)); - - break; - } - } - this->mem->barrier(); - if(this->rank == 0) - { - nancheck(rhs.at(ix::rc)(this->domain), "RHS of rc after rhs_update"); - nancheck(rhs.at(ix::rr)(this->domain), "RHS of rr after rhs_update"); - this->tend = clock::now(); - this->tupdate += std::chrono::duration_cast( this->tend - this->tbeg ); - } - } - // deals with initial supersaturation void hook_ante_loop(int nt) { @@ -151,10 +96,10 @@ class slvr_blk_1m_common : public slvr_common this->record_aux_const("r_c0", opts.r_c0); this->record_aux_const("k_acnv", opts.k_acnv); this->record_aux_const("r_eps", opts.r_eps); - this->record_aux_const("user_params rc_src", params.user_params.rc_src); - this->record_aux_const("user_params rr_src", params.user_params.rr_src); - this->record_aux_const("rt_params rc_src", params.rc_src); - this->record_aux_const("rt_params rr_src", params.rr_src); + this->record_aux_const("user_params rc_src", this->params.user_params.rc_src); + this->record_aux_const("user_params rr_src", this->params.user_params.rr_src); + this->record_aux_const("rt_params rc_src", this->params.rc_src); + this->record_aux_const("rt_params rr_src", this->params.rr_src); } } @@ -176,7 +121,12 @@ class slvr_blk_1m_common : public slvr_common parent_t::update_rhs(rhs, dt, at); // shouldnt forcings be after condensation to be consistent with lgrngn solver? + this->mem->barrier(); + if(this->rank == 0) + this->tbeg = clock::now(); + // cell-wise + // TODO: rozne cell-wise na n i n+1 ? { auto dot_rc = rhs.at(ix::rc)(this->ijk), @@ -186,6 +136,46 @@ class slvr_blk_1m_common : public slvr_common rr = this->state(ix::rr)(this->ijk); libcloudphxx::blk_1m::rhs_cellwise(opts, dot_rc, dot_rr, rc, rr); } + + // forcing + switch (at) + { + // for eulerian integration or used to init trapezoidal integration + case (0): + { + // ---- cloud water sources ---- + rc_src(); + rhs.at(ix::rc)(this->ijk) += this->alpha(this->ijk) + this->beta(this->ijk) * this->state(ix::rc)(this->ijk); + + // ---- rain water sources ---- + rr_src(); + rhs.at(ix::rr)(this->ijk) += this->alpha(this->ijk) + this->beta(this->ijk) * this->state(ix::rr)(this->ijk); + + break; + } + + case (1): + // trapezoidal rhs^n+1 + { + // ---- cloud water sources ---- + rc_src(); + rhs.at(ix::rc)(this->ijk) += this->alpha(this->ijk) + this->beta(this->ijk) * this->state(ix::rc)(this->ijk) / (1. - 0.5 * this->dt * this->beta(this->ijk)); + + // ---- rain water sources ---- + rr_src(); + rhs.at(ix::rr)(this->ijk) += this->alpha(this->ijk) + this->beta(this->ijk) * this->state(ix::rr)(this->ijk) / (1. - 0.5 * this->dt * this->beta(this->ijk)); + + break; + } + } + this->mem->barrier(); + if(this->rank == 0) + { + nancheck(rhs.at(ix::rc)(this->domain), "RHS of rc after rhs_update"); + nancheck(rhs.at(ix::rr)(this->domain), "RHS of rr after rhs_update"); + this->tend = clock::now(); + this->tupdate += std::chrono::duration_cast( this->tend - this->tbeg ); + } } // @@ -202,9 +192,16 @@ class slvr_blk_1m_common : public slvr_common struct rt_params_t : parent_t::rt_params_t { libcloudphxx::blk_1m::opts_t cloudph_opts; - bool rr_src, rc_src; }; + protected: + + // per-thread copy of params + // TODO: but slvr_common also has a copy of it's params.... + rt_params_t params; + + public: + // ctor slvr_blk_1m_common( typename parent_t::ctor_args_t args, @@ -232,6 +229,7 @@ class slvr_blk_1m< public: using parent_t = slvr_blk_1m_common; using real_t = typename ct_params_t::real_t; + using clock = typename parent_t::clock; // ctor slvr_blk_1m( @@ -249,6 +247,10 @@ class slvr_blk_1m< ) { parent_t::update_rhs(rhs, dt, at); // shouldnt forcings be after condensation to be consistent with lgrngn solver? + this->mem->barrier(); + if(this->rank == 0) + this->tbeg = clock::now(); + // column-wise for (int i = this->i.first(); i <= this->i.last(); ++i) { @@ -259,6 +261,14 @@ class slvr_blk_1m< rr = this->state(parent_t::ix::rr)(i, this->j); libcloudphxx::blk_1m::rhs_columnwise(this->opts, dot_rr, rhod, rr, this->params.dz); } + + this->mem->barrier(); + if(this->rank == 0) + { + nancheck(rhs.at(ix::rr)(this->domain), "RHS of rr after rhs_update"); + this->tend = clock::now(); + this->tupdate += std::chrono::duration_cast( this->tend - this->tbeg ); + } } }; @@ -272,6 +282,7 @@ class slvr_blk_1m< public: using parent_t = slvr_blk_1m_common; using real_t = typename ct_params_t::real_t; + using clock = typename parent_t::clock; // ctor slvr_blk_1m( @@ -289,6 +300,10 @@ class slvr_blk_1m< ) { parent_t::update_rhs(rhs, dt, at); // shouldnt forcings be after condensation to be consistent with lgrngn solver? + this->mem->barrier(); + if(this->rank == 0) + this->tbeg = clock::now(); + // column-wise for (int i = this->i.first(); i <= this->i.last(); ++i) for (int j = this->j.first(); j <= this->j.last(); ++j) @@ -300,5 +315,13 @@ class slvr_blk_1m< rr = this->state(parent_t::ix::rr)(i, j, this->k); libcloudphxx::blk_1m::rhs_columnwise(this->opts, dot_rr, rhod, rr, this->params.dz); } + + this->mem->barrier(); + if(this->rank == 0) + { + nancheck(rhs.at(ix::rr)(this->domain), "RHS of rr after rhs_update"); + this->tend = clock::now(); + this->tupdate += std::chrono::duration_cast( this->tend - this->tbeg ); + } } }; diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index 12e74195a0..b3ca4bcc34 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -20,6 +20,9 @@ class slvr_common : public slvr_dim using clock = std::chrono::high_resolution_clock; // TODO: option to disable timing, as it may affect performance a little? // timing fields + // TODO: timing slows down simulations + // either remove it and use profiling tools (e.g. vtune) + // or add some compile-time flag to turn it off clock::time_point tbeg, tend, tbeg1, tbeg_loop; std::chrono::milliseconds tdiag, tupdate, tsync, tasync, tasync_wait, tloop, tvip_rhs; @@ -311,6 +314,7 @@ class slvr_common : public slvr_dim int spinup = 0, // number of timesteps during which autoconversion is to be turned off nt; // total number of timesteps bool rv_src, th_src, uv_src, w_src, subsidence, friction, buoyancy_wet, radiation; + bool rc_src, rr_src; // these two are only relevant for blk_1m, but need to be here so that Cases can have access to it setup::arr_1D_t *th_e, *p_e, *rv_e, *th_ref, *pre_ref, *rhod, *w_LS, *hgt_fctr_sclr, *hgt_fctr_vctr; typename ct_params_t::real_t dz; // vertical grid size setup::ForceParameters_t ForceParameters; From 2108039ae5a70f9679a966b4ad4cfc84224883c3 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 26 Apr 2018 11:39:54 +0200 Subject: [PATCH 10/74] fix params copy in slvr_blk_1m --- src/slvr_blk_1m.hpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index 9d802645ce..8fd49cd46d 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -69,13 +69,13 @@ class slvr_blk_1m_common : public slvr_common // init the p_e array p_e.resize(this->shape(this->ijk)); - p_e = (*this->params.p_e)(this->vert_idx); + p_e = (*params.p_e)(this->vert_idx); p_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? // init the p_d_e array p_d_e.resize(this->shape(this->ijk)); // p_d_e = p_e - p_v_e - p_d_e = (*this->params.p_e)(this->vert_idx) - calc_p_v()((*this->params.p_e)(this->vert_idx), (*this->params.rv_e)(this->vert_idx)); + p_d_e = (*params.p_e)(this->vert_idx) - calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); p_d_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? // deal with initial supersaturation @@ -96,10 +96,10 @@ class slvr_blk_1m_common : public slvr_common this->record_aux_const("r_c0", opts.r_c0); this->record_aux_const("k_acnv", opts.k_acnv); this->record_aux_const("r_eps", opts.r_eps); - this->record_aux_const("user_params rc_src", this->params.user_params.rc_src); - this->record_aux_const("user_params rr_src", this->params.user_params.rr_src); - this->record_aux_const("rt_params rc_src", this->params.rc_src); - this->record_aux_const("rt_params rr_src", this->params.rr_src); + this->record_aux_const("user_params rc_src", params.user_params.rc_src); + this->record_aux_const("user_params rr_src", params.user_params.rr_src); + this->record_aux_const("rt_params rc_src", params.rc_src); + this->record_aux_const("rt_params rr_src", params.rr_src); } } @@ -208,6 +208,7 @@ class slvr_blk_1m_common : public slvr_common const rt_params_t &p ) : parent_t(args, p), + params(p), opts(p.cloudph_opts) {} }; @@ -265,7 +266,7 @@ class slvr_blk_1m< this->mem->barrier(); if(this->rank == 0) { - nancheck(rhs.at(ix::rr)(this->domain), "RHS of rr after rhs_update"); + nancheck(rhs.at(parent_t::ix::rr)(this->domain), "RHS of rr after rhs_update"); this->tend = clock::now(); this->tupdate += std::chrono::duration_cast( this->tend - this->tbeg ); } @@ -319,7 +320,7 @@ class slvr_blk_1m< this->mem->barrier(); if(this->rank == 0) { - nancheck(rhs.at(ix::rr)(this->domain), "RHS of rr after rhs_update"); + nancheck(rhs.at(parent_t::ix::rr)(this->domain), "RHS of rr after rhs_update"); this->tend = clock::now(); this->tupdate += std::chrono::duration_cast( this->tend - this->tbeg ); } From bb036c9f875606e24c80e416bce1338f1956f776 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 7 May 2018 15:31:18 +0200 Subject: [PATCH 11/74] constant pressure profile in lgrngn too --- src/calc_forces_common.hpp | 2 +- src/slvr_blk_1m.hpp | 24 +----------------------- src/slvr_common.hpp | 25 +++++++++++++++++++++++++ src/slvr_lgrngn.hpp | 4 +++- 4 files changed, 30 insertions(+), 25 deletions(-) diff --git a/src/calc_forces_common.hpp b/src/calc_forces_common.hpp index e80a0d9263..04876c1010 100644 --- a/src/calc_forces_common.hpp +++ b/src/calc_forces_common.hpp @@ -20,7 +20,7 @@ struct calc_T struct calc_exner { setup::real_t operator()(setup::real_t p) const - {return libcloudphxx::common::theta_dry::exner(p * si::pascals);} + {return libcloudphxx::common::theta_std::exner(p * si::pascals);} BZ_DECLARE_FUNCTOR(calc_exner) }; diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index 8fd49cd46d..934281ce73 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -6,12 +6,6 @@ #include #include -struct calc_p_v -{ - setup::real_t operator()(setup::real_t p, setup::real_t rv) const - {return libcloudphxx::common::moist_air::p_v(p * si::pascals, rv) / si::pascals;} - BZ_DECLARE_FUNCTOR2(calc_p_v) -}; template class slvr_blk_1m_common : public slvr_common @@ -25,11 +19,6 @@ class slvr_blk_1m_common : public slvr_common using clock = typename parent_t::clock; private: - // a 2D/3D arrays with copies of the environmental total pressure/partial pressure of dry air, - // it's needed by adj_cellwise_constp - blitz::Array p_e; - blitz::Array p_d_e; - void condevap() { auto @@ -41,7 +30,7 @@ class slvr_blk_1m_common : public slvr_common rhod = (*this->mem->G)(this->ijk); libcloudphxx::blk_1m::adj_cellwise_constp( - opts, rhod, p_e, p_d_e, th, rv, rc, rr, this->dt + opts, rhod, this->p_e, this->p_d_e, th, rv, rc, rr, this->dt ); this->mem->barrier(); } @@ -67,17 +56,6 @@ class slvr_blk_1m_common : public slvr_common zero_if_uninitialised(ix::rc); zero_if_uninitialised(ix::rr); - // init the p_e array - p_e.resize(this->shape(this->ijk)); - p_e = (*params.p_e)(this->vert_idx); - p_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? - - // init the p_d_e array - p_d_e.resize(this->shape(this->ijk)); - // p_d_e = p_e - p_v_e - p_d_e = (*params.p_e)(this->vert_idx) - calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); - p_d_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? - // deal with initial supersaturation condevap(); diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index b3ca4bcc34..76b76a2d1d 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -6,6 +6,15 @@ #include #include "../git_revision.h" +namespace detail +{ + struct calc_p_v + { + setup::real_t operator()(setup::real_t p, setup::real_t rv) const + {return libcloudphxx::common::moist_air::p_v(p * si::pascals, rv) / si::pascals;} + BZ_DECLARE_FUNCTOR2(calc_p_v) + }; +}; template class slvr_common : public slvr_dim @@ -43,6 +52,11 @@ class slvr_common : public slvr_dim &alpha, // 'explicit' rhs part - does not depend on the value at n+1 β // 'implicit' rhs part - coefficient of the value at n+1 + + // a 2D/3D arrays with copies of the environmental total pressure/partial pressure of dry air, + blitz::Array p_e; + blitz::Array p_d_e; + // surface precip stuff std::ofstream f_puddle; // output precipitation file @@ -347,6 +361,17 @@ class slvr_common : public slvr_dim surf_flux_sens.resize(this->shape(this->hrzntl_domain)); // TODO: resize to hrzntl_subdomain surf_flux_lat.resize(this->shape(this->hrzntl_domain)); // TODO: resize to hrzntl_subdomain r_l = 0.; + + // init the p_e array + p_e.resize(this->shape(this->ijk)); + p_e = (*params.p_e)(this->vert_idx); + p_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? + + // init the p_d_e array + p_d_e.resize(this->shape(this->ijk)); + // p_d_e = p_e - p_v_e + p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); + p_d_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? } static void alloc(typename parent_t::mem_t *mem, const int &n_iters) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index afa078567d..5a7cb3cf3a 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -284,7 +284,9 @@ class slvr_lgrngn : public slvr_common prtcls->init( make_arrinfo(this->mem->advectee(ix::th)), make_arrinfo(this->mem->advectee(ix::rv)), - make_arrinfo(rhod) + make_arrinfo(rhod), + make_arrinfo(this->p_e), + make_arrinfo(this->p_d_e) ); // writing diagnostic data for the initial condition From 3bd98b98753826849579e454ea6cc89febdf99a2 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 17 May 2018 11:47:04 +0200 Subject: [PATCH 12/74] pressure profile 3d array in lgrngn --- src/slvr_blk_1m.hpp | 17 ++++++++++++++++- src/slvr_common.hpp | 16 ---------------- src/slvr_lgrngn.hpp | 12 ++++++++++-- 3 files changed, 26 insertions(+), 19 deletions(-) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index 934281ce73..bfa463b2ab 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -19,6 +19,10 @@ class slvr_blk_1m_common : public slvr_common using clock = typename parent_t::clock; private: + // a 2D/3D arrays with copies of the environmental total pressure/partial pressure of dry air, + blitz::Array p_e; + blitz::Array p_d_e; + void condevap() { auto @@ -30,7 +34,7 @@ class slvr_blk_1m_common : public slvr_common rhod = (*this->mem->G)(this->ijk); libcloudphxx::blk_1m::adj_cellwise_constp( - opts, rhod, this->p_e, this->p_d_e, th, rv, rc, rr, this->dt + opts, rhod, p_e, p_d_e, th, rv, rc, rr, this->dt ); this->mem->barrier(); } @@ -55,6 +59,17 @@ class slvr_blk_1m_common : public slvr_common // if uninitialised fill with zeros zero_if_uninitialised(ix::rc); zero_if_uninitialised(ix::rr); + + // init the p_e array + p_e.resize(this->shape(this->ijk)); + p_e = (*params.p_e)(this->vert_idx); + p_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? + + // init the p_d_e array + p_d_e.resize(this->shape(this->ijk)); + // p_d_e = p_e - p_v_e + p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); + p_d_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? // deal with initial supersaturation condevap(); diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index 27ec9ee7df..75c6cc43f1 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -52,11 +52,6 @@ class slvr_common : public slvr_dim &alpha, // 'explicit' rhs part - does not depend on the value at n+1 β // 'implicit' rhs part - coefficient of the value at n+1 - - // a 2D/3D arrays with copies of the environmental total pressure/partial pressure of dry air, - blitz::Array p_e; - blitz::Array p_d_e; - // surface precip stuff std::ofstream f_puddle; // output precipitation file @@ -385,17 +380,6 @@ class slvr_common : public slvr_dim surf_flux_sens.resize(this->shape(this->hrzntl_domain)); // TODO: resize to hrzntl_subdomain surf_flux_lat.resize(this->shape(this->hrzntl_domain)); // TODO: resize to hrzntl_subdomain r_l = 0.; - - // init the p_e array - p_e.resize(this->shape(this->ijk)); - p_e = (*params.p_e)(this->vert_idx); - p_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? - - // init the p_d_e array - p_d_e.resize(this->shape(this->ijk)); - // p_d_e = p_e - p_v_e - p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); - p_d_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? } static void alloc(typename parent_t::mem_t *mem, const int &n_iters) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 5a7cb3cf3a..7b3fcacb83 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -281,12 +281,20 @@ class slvr_lgrngn : public slvr_common typename parent_t::arr_t rhod(this->mem->advectee(ix::th).shape()); // TODO: replace all rhod arrays with this->mem->G rhod = (*params.rhod)(this->vert_idx); + // temporary array of pressure - prtcls cant be init'd with 1D profile + typename parent_t::arr_t p_e(this->mem->advectee(ix::th).shape()); + p_e = (*params.p_e)(this->vert_idx); + + // temporary array of partial pressure of dry air - prtcls cant be init'd with 1D profile + typename parent_t::arr_t p_d_e(this->mem->advectee(ix::th).shape()); + p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); + prtcls->init( make_arrinfo(this->mem->advectee(ix::th)), make_arrinfo(this->mem->advectee(ix::rv)), make_arrinfo(rhod), - make_arrinfo(this->p_e), - make_arrinfo(this->p_d_e) + make_arrinfo(p_e), + make_arrinfo(p_d_e) ); // writing diagnostic data for the initial condition From 510b8e7593e1d04e6b67efacd07a982a3dd28b9b Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 17 May 2018 16:10:46 +0200 Subject: [PATCH 13/74] add pressure plotting --- drawbicyc/plot_fields.hpp | 10 ++++++++++ drawbicyc/plots.hpp | 3 ++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/drawbicyc/plot_fields.hpp b/drawbicyc/plot_fields.hpp index 8b4b524158..a993148263 100644 --- a/drawbicyc/plot_fields.hpp +++ b/drawbicyc/plot_fields.hpp @@ -205,6 +205,16 @@ void plot_fields(Plotter_t plotter, Plots plots) } catch(...){} } + else if (plt == "lib_pres") + { + try{ + std::string title = "libcloud pressure [Pa]"; + gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n"; + auto tmp = plotter.h5load_timestep("libcloud_pressure", at * n["outfreq"]); + plotter.plot(gp, tmp); + } + catch(...){} + } else if (plt == "supersat") { try{ diff --git a/drawbicyc/plots.hpp b/drawbicyc/plots.hpp index aa49af167f..86697184be 100644 --- a/drawbicyc/plots.hpp +++ b/drawbicyc/plots.hpp @@ -52,7 +52,8 @@ std::vector fields_dycoms({ "th", "rv", "u", "w", "sd_conc",//, "r_dry", -"RH", "supersat" +"RH", "supersat", +"lib_pres" }); std::vector fields_moist_thermal({ From 02c843b383fafab34b4163a233fc01e1f7547364 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 17 May 2018 16:11:21 +0200 Subject: [PATCH 14/74] bring back smoothing of buoyancy - otherwise occasional nans --- src/forcings.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/forcings.hpp b/src/forcings.hpp index 7e5670e2f8..d4dd86194d 100644 --- a/src/forcings.hpp +++ b/src/forcings.hpp @@ -24,8 +24,8 @@ void slvr_common::buoyancy(typename parent_t::arr_t &th, typename p (th(ijk).reindex(this->zero) - (*params.th_e)(this->vert_idx)) / (*params.th_ref)(this->vert_idx) ); -// this->smooth(tmp1, F); - F(ijk) = tmp1(ijk); + this->smooth(tmp1, F); +// F(ijk) = tmp1(ijk); } template From 56298c7eab52e89165de2a56b9c1324dbedc00cc Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 17 May 2018 16:11:58 +0200 Subject: [PATCH 15/74] save pressure + pass only env prof of p, not p_d --- src/slvr_lgrngn.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 7b3fcacb83..050a2ab02a 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -42,6 +42,10 @@ class slvr_lgrngn : public slvr_common prtcls->diag_RH(); this->record_aux("RH", prtcls->outbuf()); + // recording pressure + prtcls->diag_pressure(); + this->record_aux("libcloud_pressure", prtcls->outbuf()); + // recording precipitation rate per grid cel prtcls->diag_all(); prtcls->diag_precip_rate(); @@ -293,8 +297,7 @@ class slvr_lgrngn : public slvr_common make_arrinfo(this->mem->advectee(ix::th)), make_arrinfo(this->mem->advectee(ix::rv)), make_arrinfo(rhod), - make_arrinfo(p_e), - make_arrinfo(p_d_e) + make_arrinfo(p_e) ); // writing diagnostic data for the initial condition From e670848aef1965f6ab5dff622a7dfca152a0a1ba Mon Sep 17 00:00:00 2001 From: Maciej Waruszewski Date: Fri, 18 May 2018 10:33:05 +0200 Subject: [PATCH 16/74] fix pressure profile at the ground in moist thermal --- src/cases/MoistThermalGrabowskiClark99.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cases/MoistThermalGrabowskiClark99.hpp b/src/cases/MoistThermalGrabowskiClark99.hpp index eef4f362d0..98f58f4338 100644 --- a/src/cases/MoistThermalGrabowskiClark99.hpp +++ b/src/cases/MoistThermalGrabowskiClark99.hpp @@ -252,6 +252,7 @@ namespace setup th_e = th_std_fctr(th_e_surf)(k * dz); pre_ref(0.) = p_0 / si::pascals; + p_e(0) = pre_ref(0); T(0.) = T_0 / si::kelvins; for(int k=1; k Date: Fri, 18 May 2018 13:01:24 +0200 Subject: [PATCH 17/74] add temperature plotting --- drawbicyc/plot_fields.hpp | 10 ++++++++++ drawbicyc/plots.hpp | 2 +- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/drawbicyc/plot_fields.hpp b/drawbicyc/plot_fields.hpp index a993148263..1a4de8305b 100644 --- a/drawbicyc/plot_fields.hpp +++ b/drawbicyc/plot_fields.hpp @@ -215,6 +215,16 @@ void plot_fields(Plotter_t plotter, Plots plots) } catch(...){} } + else if (plt == "lib_temp") + { + try{ + std::string title = "libcloud temperature [K]"; + gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n"; + auto tmp = plotter.h5load_timestep("libcloud_temperature", at * n["outfreq"]); + plotter.plot(gp, tmp); + } + catch(...){} + } else if (plt == "supersat") { try{ diff --git a/drawbicyc/plots.hpp b/drawbicyc/plots.hpp index 86697184be..8d53484a80 100644 --- a/drawbicyc/plots.hpp +++ b/drawbicyc/plots.hpp @@ -53,7 +53,7 @@ std::vector fields_dycoms({ "u", "w", "sd_conc",//, "r_dry", "RH", "supersat", -"lib_pres" +"lib_pres", "lib_temp" }); std::vector fields_moist_thermal({ From fb012421df8828194ce9a471fc7ee48cad849fa9 Mon Sep 17 00:00:00 2001 From: Maciej Waruszewski Date: Fri, 18 May 2018 15:43:39 +0200 Subject: [PATCH 18/74] remove dry pressure in one-moment bulk solver --- src/slvr_blk_1m.hpp | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index bfa463b2ab..b5c6c5c817 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -19,9 +19,8 @@ class slvr_blk_1m_common : public slvr_common using clock = typename parent_t::clock; private: - // a 2D/3D arrays with copies of the environmental total pressure/partial pressure of dry air, - blitz::Array p_e; - blitz::Array p_d_e; + // a 2D/3D array with copy of the environmental total pressure of dry air + typename parent_t::arr_t &p_e; void condevap() { @@ -31,10 +30,11 @@ class slvr_blk_1m_common : public slvr_common rc = this->state(ix::rc)(this->ijk), // cloud water mixing ratio rr = this->state(ix::rr)(this->ijk); // rain water mixing ratio auto const - rhod = (*this->mem->G)(this->ijk); + rhod = (*this->mem->G)(this->ijk), + &p_e_arg = p_e(this->ijk); libcloudphxx::blk_1m::adj_cellwise_constp( - opts, rhod, p_e, p_d_e, th, rv, rc, rr, this->dt + opts, rhod, p_e_arg, th, rv, rc, rr, this->dt ); this->mem->barrier(); } @@ -61,15 +61,7 @@ class slvr_blk_1m_common : public slvr_common zero_if_uninitialised(ix::rr); // init the p_e array - p_e.resize(this->shape(this->ijk)); p_e = (*params.p_e)(this->vert_idx); - p_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? - - // init the p_d_e array - p_d_e.resize(this->shape(this->ijk)); - // p_d_e = p_e - p_v_e - p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); - p_d_e.reindexSelf(this->state(ix::rv).base()); // TODO: reindex not necessary? // deal with initial supersaturation condevap(); @@ -195,6 +187,12 @@ class slvr_blk_1m_common : public slvr_common public: + static void alloc(typename parent_t::mem_t *mem, const int &n_iters) + { + parent_t::alloc(mem, n_iters); + parent_t::alloc_tmp_sclr(mem, __FILE__, 1); // p_e + } + // ctor slvr_blk_1m_common( typename parent_t::ctor_args_t args, @@ -202,7 +200,8 @@ class slvr_blk_1m_common : public slvr_common ) : parent_t(args, p), params(p), - opts(p.cloudph_opts) + opts(p.cloudph_opts), + p_e(args.mem->tmp[__FILE__][0][0]) {} }; From 876a52e5ef8f044f442da4582b8a6dc38152addd Mon Sep 17 00:00:00 2001 From: Maciej Waruszewski Date: Fri, 18 May 2018 15:44:44 +0200 Subject: [PATCH 19/74] convert full tht to dry tht using initial rv not rv_e in moist thermal --- src/cases/MoistThermalGrabowskiClark99.hpp | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/cases/MoistThermalGrabowskiClark99.hpp b/src/cases/MoistThermalGrabowskiClark99.hpp index 98f58f4338..37b6be1792 100644 --- a/src/cases/MoistThermalGrabowskiClark99.hpp +++ b/src/cases/MoistThermalGrabowskiClark99.hpp @@ -394,7 +394,19 @@ namespace setup ), blitz::tensor::j * dz ); - + + // reinitialise th with dry conversion based on perturbed rv TODO: without loops, also in 3D + for (int i = 0; i < nx; ++i) + { + for(int k = 0; k < nz; ++k) + { + using libcloudphxx::common::theta_dry::std2dry; + using libcloudphxx::common::theta_dry::dry2std; + quantity si_rv_e = rv_e(k); + quantity si_rv = solver.advectee(ix::rv)(i, k); + solver.advectee(ix::th)(i, k) = std2dry(dry2std(th_e(k) * si::kelvins, si_rv_e), si_rv) / si::kelvins; + } + } } }; From 1cfba98a6fe0ba6e79e080d1eb3b80ac13b0729a Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 21 May 2018 16:18:05 +0200 Subject: [PATCH 20/74] DYCOMS: reference state dry density and potential temp at z=0 equal to env state's dry density and temp, not total density and temp --- src/cases/DYCOMS98.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/cases/DYCOMS98.hpp b/src/cases/DYCOMS98.hpp index acbdcc9bcb..9a0aa21b8f 100644 --- a/src/cases/DYCOMS98.hpp +++ b/src/cases/DYCOMS98.hpp @@ -229,7 +229,9 @@ namespace setup st(notopbot) = (th_e(notopbot+1) - th_e(notopbot-1)) / th_e(notopbot); real_t st_avg = blitz::sum(st) / (nz-2) / (2.*dz); // reference theta - th_ref = th_e(0) * exp(st_avg * k * dz); + // th_ref = th_e(0) * exp(st_avg * k * dz); + th_ref = th_e(0) * pow(1 + rv_e(0) / a, f) + * exp(st_avg * k * dz); // virtual temp at surface using libcloudphxx::common::moist_air::R_d_over_c_pd; using libcloudphxx::common::moist_air::c_pd; @@ -238,7 +240,8 @@ namespace setup real_t T_surf = th_e(0) * pow(p_0 / p_1000(), R_d_over_c_pd()); real_t T_virt_surf = T_surf * (1. + 0.608 * rv_e(0)); - real_t rho_surf = (p_0 / si::pascals) / T_virt_surf / 287. ; // TODO: R_d instead of 287 + real_t rho_surf = (p_0 / si::pascals) / T_virt_surf / 287. ; // TODO: R_d instead of 287, its the total, not dry density! + rho_surf /= (1 + rv_e(0)); // turn it into dry air density! TODO: is this correct? TODO2: approp change in the paper real_t cs = 9.81 / (c_pd() / si::joules * si::kilograms * si::kelvins) / st_avg / T_surf; // real_t cs = 9.81 / (c_pd() / si::joules * si::kilograms * si::kelvins) / st_avg / th_e(0); // this is correct? // rhod profile From 7e05bf3f93bfb61e6553dd08a5274abf526bb59d Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 21 May 2018 16:18:34 +0200 Subject: [PATCH 21/74] turn off smoothing in buoy and subsidence (again) --- src/forcings.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/forcings.hpp b/src/forcings.hpp index d4dd86194d..65a8c15a16 100644 --- a/src/forcings.hpp +++ b/src/forcings.hpp @@ -24,8 +24,8 @@ void slvr_common::buoyancy(typename parent_t::arr_t &th, typename p (th(ijk).reindex(this->zero) - (*params.th_e)(this->vert_idx)) / (*params.th_ref)(this->vert_idx) ); - this->smooth(tmp1, F); -// F(ijk) = tmp1(ijk); +// this->smooth(tmp1, F); + F(ijk) = tmp1(ijk); } template @@ -134,8 +134,8 @@ void slvr_common::subsidence(const int &type) // large-scale vertic this->vert_grad_cnt(tmp1, F, params.dz); F(ijk).reindex(this->zero) *= - (*params.w_LS)(this->vert_idx); - tmp1(ijk)=F(ijk); //TODO: unnecessary copy - this->smooth(tmp1, F); +// tmp1(ijk)=F(ijk); //TODO: unnecessary copy + // this->smooth(tmp1, F); } else F(ijk)=0.; From 32a3e563bb8645890211194d9d2a40e3f9aabefa Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 21 May 2018 16:19:07 +0200 Subject: [PATCH 22/74] record temperaturwe + remove old p_d_e --- src/slvr_lgrngn.hpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 050a2ab02a..17921a92bb 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -46,6 +46,10 @@ class slvr_lgrngn : public slvr_common prtcls->diag_pressure(); this->record_aux("libcloud_pressure", prtcls->outbuf()); + // recording temperature + prtcls->diag_temperature(); + this->record_aux("libcloud_temperature", prtcls->outbuf()); + // recording precipitation rate per grid cel prtcls->diag_all(); prtcls->diag_precip_rate(); @@ -289,15 +293,11 @@ class slvr_lgrngn : public slvr_common typename parent_t::arr_t p_e(this->mem->advectee(ix::th).shape()); p_e = (*params.p_e)(this->vert_idx); - // temporary array of partial pressure of dry air - prtcls cant be init'd with 1D profile - typename parent_t::arr_t p_d_e(this->mem->advectee(ix::th).shape()); - p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); - prtcls->init( make_arrinfo(this->mem->advectee(ix::th)), make_arrinfo(this->mem->advectee(ix::rv)), - make_arrinfo(rhod), - make_arrinfo(p_e) + make_arrinfo(rhod) + ,make_arrinfo(p_e) ); // writing diagnostic data for the initial condition From 4248ac54926d269e871d10c9db3bc940cd2666f8 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Tue, 22 May 2018 11:29:49 +0200 Subject: [PATCH 23/74] comment in dycoms case --- src/cases/DYCOMS98.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cases/DYCOMS98.hpp b/src/cases/DYCOMS98.hpp index 9a0aa21b8f..8901b9c34c 100644 --- a/src/cases/DYCOMS98.hpp +++ b/src/cases/DYCOMS98.hpp @@ -230,7 +230,7 @@ namespace setup real_t st_avg = blitz::sum(st) / (nz-2) / (2.*dz); // reference theta // th_ref = th_e(0) * exp(st_avg * k * dz); - th_ref = th_e(0) * pow(1 + rv_e(0) / a, f) + th_ref = th_e(0) * pow(1 + rv_e(0) / a, f) // calc dry theta at z=0 * exp(st_avg * k * dz); // virtual temp at surface using libcloudphxx::common::moist_air::R_d_over_c_pd; From 5274dd00c7e7a15ffbbcf973174b1ce778535255 Mon Sep 17 00:00:00 2001 From: Maciej Waruszewski Date: Tue, 22 May 2018 12:08:24 +0200 Subject: [PATCH 24/74] fix one-moment bulk model forces --- src/calc_forces_blk_1m.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/calc_forces_blk_1m.hpp b/src/calc_forces_blk_1m.hpp index 0174ff2620..a8d06c1263 100644 --- a/src/calc_forces_blk_1m.hpp +++ b/src/calc_forces_blk_1m.hpp @@ -12,7 +12,7 @@ void slvr_blk_1m_common::rc_src() // large-scale vertical wind parent_t::subsidence(ix::rc); - this->alpha(ijk) += this->F(ijk); + this->alpha(ijk) = this->F(ijk); } else this->alpha(ijk) = 0.; @@ -33,7 +33,7 @@ void slvr_blk_1m_common::rr_src() // large-scale vertical wind parent_t::subsidence(ix::rr); - this->alpha(ijk) += this->F(ijk); + this->alpha(ijk) = this->F(ijk); } else this->alpha(ijk) = 0.; From a935b86f0e3514b8bcac5103d3e0341d457023c0 Mon Sep 17 00:00:00 2001 From: Maciej Waruszewski Date: Tue, 22 May 2018 12:12:02 +0200 Subject: [PATCH 25/74] add cheating to one-moment bulk solver (negtozero) --- src/slvr_blk_1m.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index b5c6c5c817..913e2fec11 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -91,6 +91,11 @@ class slvr_blk_1m_common : public slvr_common void hook_ante_step() { parent_t::hook_ante_step(); + + negtozero(this->mem->advectee(ix::rv), "rv after first half of rhs"); + negtozero(this->mem->advectee(ix::rc), "rc after first half of rhs"); + negtozero(this->mem->advectee(ix::rr), "rr after first half of rhs"); + condevap(); // treat saturation adjustment as pre-advection, post-half-rhs adjustment // store rl for buoyancy //this->r_l(this->ijk) = this->state(ix::rc)(this->ijk) + this->state(ix::rr)(this->ijk); From ac4e6c2bd950461fc72511e3088655bbdc5c4d38 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Tue, 22 May 2018 16:10:21 +0200 Subject: [PATCH 26/74] fix dycoms plots for blk_1m micro --- drawbicyc/PlotterMicro.hpp | 56 +++++++++++++++++++++++++++++++- drawbicyc/plot_prof.hpp | 66 ++++++++++---------------------------- 2 files changed, 72 insertions(+), 50 deletions(-) diff --git a/drawbicyc/PlotterMicro.hpp b/drawbicyc/PlotterMicro.hpp index 7b05bc565f..3431b5ad11 100644 --- a/drawbicyc/PlotterMicro.hpp +++ b/drawbicyc/PlotterMicro.hpp @@ -20,6 +20,18 @@ class PlotterMicro_t : public Plotter_t public: // functions for diagnosing fields // + // aerosol droplets mixing ratio + auto h5load_ra_timestep( + int at + ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + { + if(this->micro == "lgrngn") + res = this->h5load_timestep("aerosol_rw_mom3", at) * 4./3. * 3.1416 * 1e3; + else if(this->micro == "blk_1m") + res = 0; + return blitz::safeToReturn(res + 0); + } + // cloud droplets mixing ratio auto h5load_rc_timestep( int at @@ -40,7 +52,7 @@ class PlotterMicro_t : public Plotter_t if(this->micro == "lgrngn") res = this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3; else if(this->micro == "blk_1m") - res = this->h5load_timestep("rc", at); + res = this->h5load_timestep("rr", at); return blitz::safeToReturn(res + 0); } @@ -59,6 +71,48 @@ class PlotterMicro_t : public Plotter_t return blitz::safeToReturn(res + 0); } + // cloud droplets concentration [1/kg] + auto h5load_nc_timestep( + int at + ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + { + if(this->micro == "lgrngn") + res = this->h5load_timestep("cloud_rw_mom0", at); + else if(this->micro == "blk_1m") + res = 0; + return blitz::safeToReturn(res + 0); + } + + // precipitation flux [W/m2] + auto h5load_prflux_timestep( + int at + ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + { + if(this->micro == "lgrngn") + { + res = this->h5load_timestep("precip_rate", at) + * 4./3 * 3.14 * 1e3 // to get mass + / this->CellVol // averaged over cell volume, TODO: make precip rate return specific moment? wouldnt need the dx and dy + * 2264.76e3; // latent heat of evaporation [J/kg] + } + else if(this->micro == "blk_1m") + res = 0; + return blitz::safeToReturn(res + 0); + } + + // RH + auto h5load_RH_timestep( + int at + ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + { + if(this->micro == "lgrngn") + res = this->h5load_timestep("RH", at); + else if(this->micro == "blk_1m") + res = 0; + return blitz::safeToReturn(res + 0); + } + + // functions for diagnosing statistics // helper function that calculates staistics (mean and std_dev) of a field only in cloudy cells diff --git a/drawbicyc/plot_prof.hpp b/drawbicyc/plot_prof.hpp index 7b541df6a3..a709eaee79 100644 --- a/drawbicyc/plot_prof.hpp +++ b/drawbicyc/plot_prof.hpp @@ -71,10 +71,10 @@ void plot_profiles(Plotter_t plotter, Plots plots) std::cout << at * n["outfreq"] << std::endl; if (plt == "rliq") { - // liquid water content (cloud + rain, missing droplets with r<0.5um!) - //res += plotter.h5load_timestep("aerosol_rw_mom3", at * n["outfreq"]) * 4./3 * 3.1416 * 1e3 * 1e3; - res += plotter.h5load_timestep("cloud_rw_mom3", at * n["outfreq"]) * 4./3 * 3.1416 * 1e3 * 1e3; - res += plotter.h5load_timestep("rain_rw_mom3", at * n["outfreq"]) * 4./3 * 3.1416 * 1e3 * 1e3; + // liquid water content + res += plotter.h5load_ra_timestep(at * n["outfreq"]) * 1e3; // aerosol + res += plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud + res += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; // rain gp << "set title 'liquid water [g/kg]'\n"; } if (plt == "gccn_rw") @@ -332,7 +332,7 @@ void plot_profiles(Plotter_t plotter, Plots plots) else if (plt == "sat_RH") { { - auto tmp = plotter.h5load_timestep("RH", at * n["outfreq"]); + auto tmp = plotter.h5load_RH_timestep(at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); res_tmp = (snap -1) * 100; } @@ -360,26 +360,11 @@ void plot_profiles(Plotter_t plotter, Plots plots) } else if (plt == "00rtot") { - { - auto tmp = plotter.h5load_timestep("aerosol_rw_mom3", at * n["outfreq"]) * 4./3 * 3.1416 * 1e3 * 1e3; - typename Plotter_t::arr_t snap(tmp); - res_tmp = snap; - } - { - auto tmp = plotter.h5load_timestep("cloud_rw_mom3", at * n["outfreq"]) * 4./3 * 3.1416 * 1e3 * 1e3; - typename Plotter_t::arr_t snap(tmp); - res_tmp += snap; - } - { - auto tmp = plotter.h5load_timestep("rain_rw_mom3", at * n["outfreq"]) * 4./3 * 3.1416 * 1e3 * 1e3; - typename Plotter_t::arr_t snap(tmp); - res_tmp += snap; - } - { - auto tmp = plotter.h5load_timestep("rv", at * n["outfreq"]) * 1e3; - typename Plotter_t::arr_t snap(tmp); - res_tmp += snap; - } + res_tmp = plotter.h5load_ra_timestep(at * n["outfreq"]) * 1e3; // aerosol + res_tmp += plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud + res_tmp += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; // rain + res_tmp += plotter.h5load_timestep("rv", at * n["outfreq"]) * 1e3; // vapour + res += res_tmp; res_prof = plotter.horizontal_mean(res_tmp); // average in x // find instantaneous inversion height @@ -390,7 +375,7 @@ void plot_profiles(Plotter_t plotter, Plots plots) { // cloud drops concentration [1/cm^3] { - auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); + auto tmp = plotter.h5load_nc_timestep(at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); snap *= rhod; // b4 it was specific moment snap /= 1e6; // per cm^3 @@ -404,7 +389,7 @@ void plot_profiles(Plotter_t plotter, Plots plots) try { // cloud fraction (cloudy if N_c > 20/cm^3) - auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); + auto tmp = plotter.h5load_nc_timestep(at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); snap *= rhod; // b4 it was specific moment snap /= 1e6; // per cm^3 @@ -426,21 +411,9 @@ void plot_profiles(Plotter_t plotter, Plots plots) { auto &ql(res_tmp2); auto &T(res_tmp); - { - auto tmp = plotter.h5load_timestep("aerosol_rw_mom3", at * n["outfreq"]) * 4./3 * 3.1416 * 1e3; - typename Plotter_t::arr_t snap(tmp); - ql = snap; - } - { - auto tmp = plotter.h5load_timestep("cloud_rw_mom3", at * n["outfreq"]) * 4./3 * 3.14 * 1e3; - typename Plotter_t::arr_t snap(tmp); - ql += snap; - } - { - auto tmp = plotter.h5load_timestep("rain_rw_mom3", at * n["outfreq"]) * 4./3 * 3.14 * 1e3; - typename Plotter_t::arr_t snap(tmp); - ql += snap; - } + ql = plotter.h5load_ra_timestep(at * n["outfreq"]); // aerosol + ql += plotter.h5load_rc_timestep(at * n["outfreq"]); // cloud + ql += plotter.h5load_rr_timestep(at * n["outfreq"]); // rain // ql is now q_l (liq water content) // auto tmp = plotter.h5load_timestep("th", at * n["outfreq"]); // typename Plotter_t::arr_t th_d(tmp); @@ -461,7 +434,7 @@ void plot_profiles(Plotter_t plotter, Plots plots) { // cloud fraction (cloudy if N_c > 20/cm^3) { - auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); + auto tmp = plotter.h5load_nc_timestep(at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); snap *= rhod; // b4 it was specific moment snap /= 1e6; // per cm^3 @@ -474,12 +447,7 @@ void plot_profiles(Plotter_t plotter, Plots plots) { // precipitation flux(doesnt include vertical volicty w!) { - auto tmp = plotter.h5load_timestep("precip_rate", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap = snap * 4./3 * 3.14 * 1e3 // to get mass - / plotter.CellVol // averaged over cell volume, TODO: make precip rate return specific moment? wouldnt need the dx and dy - * 2264.76e3; // latent heat of evaporation [J/kg] - res += snap; + res += plotter.h5load_prflux_timestep(at * n["outfreq"]); } // add vertical velocity to precipitation flux (3rd mom of cloud drops * w) /* From 673088e63d2e818289c08373444859ecc82894b3 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 23 May 2018 11:56:48 +0200 Subject: [PATCH 27/74] blk_1m: save puddle + use newton raphson condensation --- src/slvr_blk_1m.hpp | 50 +++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 46 insertions(+), 4 deletions(-) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index 913e2fec11..eb9ec2365d 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -33,9 +33,14 @@ class slvr_blk_1m_common : public slvr_common rhod = (*this->mem->G)(this->ijk), &p_e_arg = p_e(this->ijk); +/* libcloudphxx::blk_1m::adj_cellwise_constp( opts, rhod, p_e_arg, th, rv, rc, rr, this->dt ); +*/ + libcloudphxx::blk_1m::adj_cellwise_nwtrph( + opts, p_e_arg, th, rv, rc, this->dt + ); this->mem->barrier(); } @@ -46,12 +51,42 @@ class slvr_blk_1m_common : public slvr_common } protected: + + // accumulated water falling out of domain + real_t puddle; + + void diag() + { + assert(this->rank == 0); + parent_t::tbeg = parent_t::clock::now(); + + // recording puddle + for(int i=0; i < 10; ++i) + { + this->f_puddle << i << " " << (i == 8 ? this->puddle : 0) << "\n"; + } + this->f_puddle << "\n"; + + parent_t::tend = parent_t::clock::now(); + parent_t::tdiag += std::chrono::duration_cast( parent_t::tend - parent_t::tbeg ); + } + void rc_src(); void rr_src(); bool get_rain() { return opts.conv; } void set_rain(bool val) { opts.conv = val; }; + void record_all() + { + // plain (no xdmf) hdf5 output + parent_t::output_t::parent_t::record_all(); + // UWLCM output + diag(); + // xmf markup + this->write_xmfs(); + } + // deals with initial supersaturation void hook_ante_loop(int nt) @@ -119,12 +154,18 @@ class slvr_blk_1m_common : public slvr_common // TODO: rozne cell-wise na n i n+1 ? { auto + dot_th = rhs.at(ix::th)(this->ijk), + dot_rv = rhs.at(ix::rv)(this->ijk), dot_rc = rhs.at(ix::rc)(this->ijk), dot_rr = rhs.at(ix::rr)(this->ijk); const auto + th = this->state(ix::th)(this->ijk), + rv = this->state(ix::rv)(this->ijk), rc = this->state(ix::rc)(this->ijk), - rr = this->state(ix::rr)(this->ijk); - libcloudphxx::blk_1m::rhs_cellwise(opts, dot_rc, dot_rr, rc, rr); + rr = this->state(ix::rr)(this->ijk), + rhod = (*this->mem->G)(this->ijk), + &p_e_arg = p_e(this->ijk); + libcloudphxx::blk_1m::rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, rhod, p_e_arg, th, rv, rc, rr); } // forcing @@ -206,6 +247,7 @@ class slvr_blk_1m_common : public slvr_common parent_t(args, p), params(p), opts(p.cloudph_opts), + puddle(0), p_e(args.mem->tmp[__FILE__][0][0]) {} }; @@ -257,7 +299,7 @@ class slvr_blk_1m< const auto rhod = (*this->mem->G)(i, this->j), rr = this->state(parent_t::ix::rr)(i, this->j); - libcloudphxx::blk_1m::rhs_columnwise(this->opts, dot_rr, rhod, rr, this->params.dz); + this->puddle += - libcloudphxx::blk_1m::rhs_columnwise(this->opts, dot_rr, rhod, rr, this->params.dz); } this->mem->barrier(); @@ -311,7 +353,7 @@ class slvr_blk_1m< const auto rhod = (*this->mem->G)(i, j, this->k), rr = this->state(parent_t::ix::rr)(i, j, this->k); - libcloudphxx::blk_1m::rhs_columnwise(this->opts, dot_rr, rhod, rr, this->params.dz); + this->puddle += - libcloudphxx::blk_1m::rhs_columnwise(this->opts, dot_rr, rhod, rr, this->params.dz); } this->mem->barrier(); From 55d2cd1c72164ba99cdd399bec4a89f92b7264c8 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 23 May 2018 12:01:45 +0200 Subject: [PATCH 28/74] negtozero: output suppressed in production runs --- src/detail/checknan.cpp | 39 ++++++++++++++++++++++++++------------- src/slvr_blk_1m.hpp | 6 +++--- 2 files changed, 29 insertions(+), 16 deletions(-) diff --git a/src/detail/checknan.cpp b/src/detail/checknan.cpp index 10f2781f7a..24840a86ab 100644 --- a/src/detail/checknan.cpp +++ b/src/detail/checknan.cpp @@ -18,14 +18,12 @@ #define negcheck(arr, name) {nancheck_hlprs::negcheck_hlpr(arr, name);} #endif +#ifdef NDEBUG // actually not to zero, but to 1e-10 (we need rv>0 in libcloud and cond substepping numerical errors colud lead to rv<0 if we would set it here to 0) -// we don't make it a critical section, because it is also used in production runs -#define negtozero(arr, name) {if(min(arr) < 0.) {\ - std::cout << "A negative number detected in: " << name << std::endl;\ - std::cout << arr;\ - std::cout << "CHEATING: turning negative values to small positive values" << std::endl;\ - arr = where(arr <= 0., 1e-10, arr);\ - }} +#define negtozero(arr, name) {arr = where(arr <= 0., 1e-10, arr);} +#else +#define negtozero(arr, name) {nancheck_hlprs::negtozero_hlpr(arr, name);} +#endif #ifndef NDEBUG namespace nancheck_hlprs @@ -33,9 +31,9 @@ namespace nancheck_hlprs template void nancheck_hlpr(const arr_t &arr, std::string name) { - #pragma omp critical + if(!std::isfinite(sum(arr))) { - if(!std::isfinite(sum(arr))) + #pragma omp critical { std::cout << "A not-finite number detected in: " << name << std::endl; std::cout << arr; @@ -47,9 +45,9 @@ namespace nancheck_hlprs template void nancheck2_hlpr(const arr_t &arrcheck, const arr_t &arrout, std::string name) { - #pragma omp critical + if(!std::isfinite(sum(arrcheck))) { - if(!std::isfinite(sum(arrcheck))) + #pragma omp critical { std::cout << "A not-finite number detected in: " << name << std::endl; std::cout << arrcheck; @@ -62,9 +60,9 @@ namespace nancheck_hlprs template void negcheck_hlpr(const arr_t &arr, std::string name) { - #pragma omp critical + if(min(arr) < 0.) { - if(min(arr) < 0.) + #pragma omp critical { std::cout << "A negative number detected in: " << name << std::endl; std::cout << arr; @@ -72,5 +70,20 @@ namespace nancheck_hlprs } } } + + template + void negtozero_hlpr(arr_t arr, std::string name) + { + if(min(arr) < 0.) + { + #pragma omp critical + { + std::cout << "A negative number detected in: " << name << std::endl; + std::cout << arr; + std::cout << "CHEATING: turning negative values to small positive values" << std::endl; + } + arr = where(arr <= 0., 1e-10, arr); + } + } }; #endif diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index eb9ec2365d..e651b64eb8 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -127,9 +127,9 @@ class slvr_blk_1m_common : public slvr_common { parent_t::hook_ante_step(); - negtozero(this->mem->advectee(ix::rv), "rv after first half of rhs"); - negtozero(this->mem->advectee(ix::rc), "rc after first half of rhs"); - negtozero(this->mem->advectee(ix::rr), "rr after first half of rhs"); + negtozero(this->mem->advectee(ix::rv)(this->ijk), "rv after first half of rhs"); + negtozero(this->mem->advectee(ix::rc)(this->ijk), "rc after first half of rhs"); + negtozero(this->mem->advectee(ix::rr)(this->ijk), "rr after first half of rhs"); condevap(); // treat saturation adjustment as pre-advection, post-half-rhs adjustment // store rl for buoyancy From 735ee62a74096666261116b428c63a73fe8059f0 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 24 May 2018 14:42:33 +0200 Subject: [PATCH 29/74] dycoms case output lwp --- src/cases/DYCOMS98.hpp | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/cases/DYCOMS98.hpp b/src/cases/DYCOMS98.hpp index 8901b9c34c..cfb3a6c574 100644 --- a/src/cases/DYCOMS98.hpp +++ b/src/cases/DYCOMS98.hpp @@ -199,6 +199,8 @@ namespace setup real_t c = l_tri() / c_pd() / si::kelvins; real_t d = l_tri() / si::joules * si::kilograms / rv; real_t f = R_d_over_c_pd(); + +real_t lwp_env = 0; for(int k=1; k() / si::pascals), f); + + bottom = R_d() / si::joules * si::kelvins * si::kilograms * T(k) * (1 + 0.61 * rv_e(k)); // (p / rho) of moist air at k-1 + rho1 = p_e(k) / bottom; // rho at k-1 + lwp_env += delta * rho1; +std::cout << k << " env_prof temp: " << T(k) << " env prof delta: " << delta << std::endl; } + lwp_env = lwp_env * 5 * 1e3; +std::cout << "lwp env: " << lwp_env << std::endl; // compute reference state theta and rhod blitz::firstIndex k; @@ -239,18 +248,22 @@ namespace setup using libcloudphxx::common::theta_std::p_1000; real_t T_surf = th_e(0) * pow(p_0 / p_1000(), R_d_over_c_pd()); +/* real_t T_virt_surf = T_surf * (1. + 0.608 * rv_e(0)); real_t rho_surf = (p_0 / si::pascals) / T_virt_surf / 287. ; // TODO: R_d instead of 287, its the total, not dry density! rho_surf /= (1 + rv_e(0)); // turn it into dry air density! TODO: is this correct? TODO2: approp change in the paper - real_t cs = 9.81 / (c_pd() / si::joules * si::kilograms * si::kelvins) / st_avg / T_surf; - // real_t cs = 9.81 / (c_pd() / si::joules * si::kilograms * si::kelvins) / st_avg / th_e(0); // this is correct? - // rhod profile - rhod = rho_surf * exp(- st_avg * k * dz) * pow( - 1. - cs * (1 - exp(- st_avg * k * dz)), (1. / R_d_over_c_pd()) - 1); +*/ + real_t rho_surf = (p_0 / si::pascals) / T_surf / (1. + 29. / 18. * rv_e(0)) / 287. ; // dry air density at the surface TODO: R_d instead of 287 // theta_std env prof to theta_dry_e for(int k=1; k(th_e(k) * si::kelvins, quantity(rv_e(k))) / si::kelvins; + + // real_t cs = 9.81 / (c_pd() / si::joules * si::kilograms * si::kelvins) / st_avg / T_surf; // this is from Wojtek + real_t cs = 9.81 / (c_pd() / si::joules * si::kilograms * si::kelvins) / st_avg / th_e(0); // this is correct? or total, not dry th_e(0) should be here? + // rhod profile + rhod = rho_surf * exp(- st_avg * k * dz) * pow( + 1. - cs * (1 - exp(- st_avg * k * dz)), (1. / R_d_over_c_pd()) - 1); // subsidence rate w_LS = w_LS_fctr()(k * dz); From e733c5ce8737f84a296fc3286bc6098e852b687b Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Fri, 8 Jun 2018 13:24:12 +0200 Subject: [PATCH 30/74] pass p_d_e to lgrngn micro init --- src/slvr_lgrngn.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 17921a92bb..d321217225 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -292,12 +292,16 @@ class slvr_lgrngn : public slvr_common // temporary array of pressure - prtcls cant be init'd with 1D profile typename parent_t::arr_t p_e(this->mem->advectee(ix::th).shape()); p_e = (*params.p_e)(this->vert_idx); - + + // temporary array of partial pressure of dry air - prtcls cant be init'd with 1D profile + typename parent_t::arr_t p_d_e(this->mem->advectee(ix::th).shape()); + p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); prtcls->init( make_arrinfo(this->mem->advectee(ix::th)), make_arrinfo(this->mem->advectee(ix::rv)), make_arrinfo(rhod) ,make_arrinfo(p_e) + ,make_arrinfo(p_d_e) ); // writing diagnostic data for the initial condition From 1ac61eb638789c562bedcd4a1b8a1a18e63940a5 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Fri, 8 Jun 2018 13:27:31 +0200 Subject: [PATCH 31/74] adapt moist thermal to the constant pressure profile --- src/bicycles.cpp | 9 +- src/cases/CasesCommon.hpp | 2 +- src/cases/DYCOMS98.hpp | 4 +- src/cases/DryThermalGMD2015.hpp | 4 +- src/cases/LasherTrapp2001.hpp | 4 +- src/cases/MoistThermalGrabowskiClark99.hpp | 164 ++++++--------------- src/slvr_common.hpp | 10 -- 7 files changed, 52 insertions(+), 145 deletions(-) diff --git a/src/bicycles.cpp b/src/bicycles.cpp index 8bcf12bcf6..55b005403c 100644 --- a/src/bicycles.cpp +++ b/src/bicycles.cpp @@ -124,17 +124,17 @@ void run(int nx, int nz, const user_params_t &user_params) if(user_params.model_case == "dry_thermal") { concurr.reset(new concurr_openmp_cyclic_t(p)); - case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, user_params.rng_seed); // works only by chance? + case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, p_e, user_params.rng_seed); // works only by chance? } else if(user_params.model_case == "lasher_trapp") { concurr.reset(new concurr_openmp_rigid_rigid_t(p)); - case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, user_params.rng_seed); // works only by chance? + case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, p_e, user_params.rng_seed); // works only by chance? } else { concurr.reset(new concurr_openmp_rigid_t(p)); - case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, user_params.rng_seed); + case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, p_e, user_params.rng_seed); } // setup panic pointer and the signal handler @@ -389,7 +389,6 @@ void run_hlpr(bool piggy, std::string type, Args&&... args) struct ct_params_final : ct_params_piggy { enum { opts = opts::nug | opts::iga | opts::fct }; }; run>(args...); } - } else // piggybacking { @@ -520,7 +519,7 @@ int main(int argc, char** argv) else if (micro == "blk_1m" && ny > 0) // 3D one-moment run_hlpr(piggy, user_params.model_case, nx, ny, nz, user_params); - // TODO: not only micro can be wrong + // TODO: not only micro can be wrong else throw po::validation_error( po::validation_error::invalid_option_value, micro, "micro" diff --git a/src/cases/CasesCommon.hpp b/src/cases/CasesCommon.hpp index 2ed57978c5..b518a5ca8a 100644 --- a/src/cases/CasesCommon.hpp +++ b/src/cases/CasesCommon.hpp @@ -63,7 +63,7 @@ namespace setup virtual void setopts(typename concurr_t::solver_t::rt_params_t ¶ms, int nx, int nz, const user_params_t &user_params) {assert(false);}; virtual void setopts(typename concurr_t::solver_t::rt_params_t ¶ms, int nx, int ny, int nz, const user_params_t &user_params) {assert(false);}; - virtual void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, int rng_seed) =0; + virtual void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &p_e, int rng_seed) =0; virtual void env_prof(arr_1D_t &th_e, arr_1D_t &p_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) =0; virtual void update_surf_flux_sens(typename concurr_t::solver_t::arr_sub_t &surf_flux_sens, int timestep, real_t dt) {if(timestep==0) surf_flux_sens = 0.;}; virtual void update_surf_flux_lat(typename concurr_t::solver_t::arr_sub_t &surf_flux_lat, int timestep, real_t dt) {if(timestep==0) surf_flux_lat = 0.;}; diff --git a/src/cases/DYCOMS98.hpp b/src/cases/DYCOMS98.hpp index cfb3a6c574..e51a1e2110 100644 --- a/src/cases/DYCOMS98.hpp +++ b/src/cases/DYCOMS98.hpp @@ -325,7 +325,7 @@ std::cout << "lwp env: " << lwp_env << std::endl; params.dz = params.dj; } - void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, int rng_seed) + void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &p_e, int rng_seed) { blitz::secondIndex k; this->intcond_hlpr(solver, rhod, rng_seed, k); @@ -346,7 +346,7 @@ std::cout << "lwp env: " << lwp_env << std::endl; params.dz = params.dk; } - void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, int rng_seed) + void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &p_e, int rng_seed) { blitz::thirdIndex k; this->intcond_hlpr(solver, rhod, rng_seed, k); diff --git a/src/cases/DryThermalGMD2015.hpp b/src/cases/DryThermalGMD2015.hpp index 6fbaf7d71f..83672b1426 100644 --- a/src/cases/DryThermalGMD2015.hpp +++ b/src/cases/DryThermalGMD2015.hpp @@ -110,7 +110,7 @@ namespace setup } // function expecting a libmpdata++ solver as argument - void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, int rng_seed) + void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &p_e, int rng_seed) { blitz::secondIndex k; this->intcond_hlpr(solver, rhod, rng_seed, k); @@ -133,7 +133,7 @@ namespace setup } // function expecting a libmpdata++ solver as argument - void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, int rng_seed) + void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &p_e, int rng_seed) { blitz::thirdIndex k; this->intcond_hlpr(solver, rhod, rng_seed, k); diff --git a/src/cases/LasherTrapp2001.hpp b/src/cases/LasherTrapp2001.hpp index 42142f1cc6..9adbf8adfb 100644 --- a/src/cases/LasherTrapp2001.hpp +++ b/src/cases/LasherTrapp2001.hpp @@ -304,7 +304,7 @@ namespace setup params.dz = params.dj; } - void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, int rng_seed) + void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &p_e, int rng_seed) { blitz::secondIndex k; this->intcond_hlpr(solver, rhod, rng_seed, k); @@ -325,7 +325,7 @@ namespace setup params.dz = params.dk; } - void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, int rng_seed) + void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &p_e, int rng_seed) { blitz::thirdIndex k; this->intcond_hlpr(solver, rhod, rng_seed, k); diff --git a/src/cases/MoistThermalGrabowskiClark99.hpp b/src/cases/MoistThermalGrabowskiClark99.hpp index 37b6be1792..592fea8adf 100644 --- a/src/cases/MoistThermalGrabowskiClark99.hpp +++ b/src/cases/MoistThermalGrabowskiClark99.hpp @@ -4,6 +4,16 @@ #pragma once #include "CasesCommon.hpp" +namespace detail +{ + struct calc_p_v + { + setup::real_t operator()(setup::real_t p, setup::real_t rv) const + {return libcloudphxx::common::moist_air::p_v(p * si::pascals, rv) / si::pascals;} + BZ_DECLARE_FUNCTOR2(calc_p_v) + }; +}; + namespace setup { namespace moist_thermal @@ -27,6 +37,21 @@ namespace setup { return moist_air::eps() * RH * const_cp::p_vs(T) / (p - RH * const_cp::p_vs(T)); } + +/* + // rv(RH, th_dry, rhod) + real_t RH_th_rhod_to_rv(const real_t &RH, const real_t &th_dry, const real_t &rhod) + { + real_t T = theta_dry::T(th_dry * si::kelvins, (rhod * si::kilograms / si::cubic_metres)) / si::kelvins; + return (p_vs(T * si::kelvins) / si::pascals) * RH / (rhod * T * (R_v() / si::joules * si::kilograms * si::kelvins)); + } + + // rv(RH, T, rhod) + real_t RH_T_rhod_to_rv(const real_t &RH, const real_t &T, const real_t &rhod) + { + return (p_vs(T * si::kelvins) / si::pascals) * RH / (rhod * T * (R_v() / si::joules * si::kilograms * si::kelvins)); + } +*/ const quantity T_0(283. * si::kelvins); // surface temperature const quantity p_0(85000 * si::pascals); // total surface temperature @@ -37,15 +62,11 @@ namespace setup const quantity rv_0 = RH_T_p_to_rv(env_RH, T_0, p_0); const quantity qv_0 = rv_0 / (1. + rv_0); // specific humidity at surface const quantity -// z_0 ( 0 * si::metres), - Z ( 2400 * si::metres), // DYCOMS: 1500 - X ( 3600 * si::metres), // DYCOMS: 6400 - Y ( 3600 * si::metres), // DYCOMS: 6400 + Z ( 2400 * si::metres), + X ( 3600 * si::metres), + Y ( 3600 * si::metres), z_prtrb ( 800 * si::metres); -// const setup::real_t rhod_surf = (p_0 / si::pascals) / (T_0 / si::kelvins) /( R_d() / si::joules * si::kelvins * si::kilograms + rv_0 * R_v() / si::joules * si::kelvins * si::kilograms); - // const setup::real_t rhod_surf = (p_0 / si::pascals) / (T_0 / si::kelvins) /( R_d() / si::joules * si::kelvins * si::kilograms) / (1+0.608 * qv_0); const setup::real_t rhod_surf = theta_std::rhod(p_0, th_std_0, rv_0) * si::cubic_metres / si::kilograms; - //const setup::real_t rhod_surfW = (p_0 / si::pascals) / (T_0 / si::kelvins) /( R_d() / si::joules * si::kelvins * si::kilograms); const setup::real_t cs = (libcloudphxx::common::earth::g() / si::metres_per_second_squared) / (c_pd() / si::joules * si::kilograms * si::kelvins) / stab / (T_0 / si::kelvins); const real_t z_abs = 125000; // [m] height above which absorber works, no absorber @@ -65,21 +86,6 @@ namespace setup BZ_DECLARE_FUNCTOR(th_std_fctr); }; -/* - // temperature profile (constant stability atmosphere, Clark Farley 1984) - real_t T(const real_t &z) - { - return (T_0 / si::kelvins) / exp(- stab * z) * ( - 1. - cs * (1 - exp(- stab * z))); - } - - // pressure profile (constant stability atmosphere, Clark Farley 1984), dry... - real_t p(const real_t &z) - { - return p_0 / si::pascals * pow( 1. - cs * (1 - exp(- stab * z)), 1. / R_d_over_c_pd() ); - } - -*/ // air density profile (constant stability atmosphere, Clark Farley 1984), dry... struct rho_fctr { @@ -95,18 +101,6 @@ namespace setup BZ_DECLARE_FUNCTOR(rho_fctr); }; - // rv(RH, th_dry, rhod) - real_t RH_th_rhod_to_rv(const real_t &RH, const real_t &th_dry, const real_t &rhod) - { - real_t T = theta_dry::T(th_dry * si::kelvins, (rhod * si::kilograms / si::cubic_metres)) / si::kelvins; - return (p_vs(T * si::kelvins) / si::pascals) * RH / (rhod * T * (R_v() / si::joules * si::kilograms * si::kelvins)); - } - - // rv(RH, T, rhod) - real_t RH_T_rhod_to_rv(const real_t &RH, const real_t &T, const real_t &rhod) - { - return (p_vs(T * si::kelvins) / si::pascals) * RH / (rhod * T * (R_v() / si::joules * si::kilograms * si::kelvins)); - } struct RH { @@ -122,32 +116,17 @@ namespace setup BZ_DECLARE_FUNCTOR(RH); }; -/* - struct env_rv - { - quantity operator()(const real_t &z) const - { - // return RH_T_rhod_to_rv(env_RH, T(z), rhod_fctr()(z)); - return RH_th_rhod_to_rv(env_RH, th_std_fctr(th_std_0 / si::kelvins)(z) , rho_fctr(rhod_surf)(z)); - // return RH_T_p_to_rv(env_RH, T(z) * si::kelvins, p(z) * si::pascals); - } - BZ_DECLARE_FUNCTOR(env_rv); - }; - */ - struct prtrb_rv { - arr_1D_t &_th_std, &_rhod; + arr_1D_t &_T, &_p; real_t dz; - prtrb_rv(arr_1D_t _th_std, arr_1D_t _rhod, real_t dz): _th_std(_th_std), _rhod(_rhod), dz(dz) {} + prtrb_rv(arr_1D_t _T, arr_1D_t _p, real_t dz): _T(_T), _p(_p), dz(dz) {} quantity operator()(const real_t &r, const real_t &z) const { - return RH_th_rhod_to_rv(RH()(r), this->_th_std(z/this->dz) , this->_rhod(z/this->dz)); - // return RH_T_p_to_rv(env_RH, T(z) * si::kelvins, p(z) * si::pascals); - // return RH_T_rhod_to_rv(env_RH, T(z), rhod_fctr()(z)); + return RH_T_p_to_rv(RH()(r), this->_T(z/this->dz) * si::kelvins , this->_p(z/this->dz) * si::pascals); } - BZ_DECLARE_FUNCTOR2(prtrb_rv); + BZ_DECLARE_FUNCTOR2(prtrb_rv); }; // its in fact the moist thermal from our 2017 GMD paper on Twomey SDs? differences: kappa=1.28, i.e. sea salt aerosol @@ -186,8 +165,6 @@ namespace setup int nx = solver.advectee().extent(0); // ix::w is the index of vertical domension both in 2D and 3D real_t dx = (X / si::metres) / (nx-1); - // solver.advectee(ix::rv) = rv_e(index); - solver.advectee(ix::u) = 0; solver.advectee(ix::w) = 0; @@ -201,7 +178,6 @@ namespace setup // initial potential temperature solver.advectee(ix::th) = th_e(index); -//libcloudphxx::common::theta_dry::std2dry(th_std_0, rv_0) / si::kelvins; } @@ -248,7 +224,7 @@ namespace setup //rv_e(0) = env_RH * qvs; rv_e(0) = rv_0;// env_RH * qvs; real_t th_e_surf = th_std_0 / si::kelvins * (1 + a * rv_e(0)); // virtual potential temp - + th_e = th_std_fctr(th_e_surf)(k * dz); pre_ref(0.) = p_0 / si::pascals; @@ -294,7 +270,7 @@ namespace setup T(k)=th_e(k)* pow(pre_ref(k)/1.e5, cap); T(k)=T(k)/(1.+a*rv_e(k)); } - rv_e(k) = RH_T_p_to_rv(env_RH, T(k) * si::kelvins, pre_ref(k) * si::pascals); // cheating! + //rv_e(k) = RH_T_p_to_rv(env_RH, T(k) * si::kelvins, pre_ref(k) * si::pascals); // cheating! p_e(k) = pre_ref(k); } @@ -309,6 +285,7 @@ namespace setup { quantity si_rv_e = rv_e(k); th_e(k) = libcloudphxx::common::theta_dry::std2dry(th_e(k) * si::kelvins, si_rv_e) / si::kelvins; + real_t p_d = pre_ref(k) - libcloudphxx::common::moist_air::p_v(pre_ref(k) * si::pascals, rv_e(k)) / si::pascals; } th_ref = th_e;//th_std_fctr(th_std_0 / si::kelvins)(k * dz); } @@ -318,48 +295,6 @@ namespace setup { this->kappa = 1.28; // NaCl aerosol } - - // calculate the initial environmental theta and rv profiles - /* - template - void env_prof(arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &th_ref, arr_1D_t &pre_ref, arr_1D_t &rhod, arr_1D_t &w_LS, arr_1D_t &hgt_fctr_vctr, arr_1D_t &hgt_fctr_sclr, int nz, const user_params_t &user_params) - { - setup::real_t dz = (Z / si::metres) / (nz-1); - blitz::firstIndex k; - th_ref = th_std_fctr()(k * dz); - th_e = th_ref; - rv_e = env_rv()(k * dz); - rhod = rho_fctr()(k * dz); - - std::cout << "th ref: " << th_ref << std::endl; - std::cout << "th e: " << th_e << std::endl; - std::cout << "rv e: " << rv_e << std::endl; - std::cout << "rho_ref: " << rhod << std::endl; - - // subsidence rate - //w_LS = setup::w_LS_fctr()(k * dz); - - // surface sources relaxation factors - // for vectors - - // real_t z_0 = setup::z_rlx_vctr / si::metres; - // hgt_fctr_vctr = exp(- (k-0.5) * dz / z_0); // z=0 at k=1/2 - // hgt_fctr_vctr(0) = 1; - // // for scalars - // z_0 = user_params.z_rlx_sclr; - // hgt_fctr_sclr = exp(- (k-0.5) * dz / z_0); - // hgt_fctr_sclr(0) = 1; - - - // calc divergence directly - real_t z_0 = z_rlx_vctr / si::metres; - hgt_fctr_vctr = exp(- k * dz / z_0) / z_0; - // for scalars - z_0 = user_params.z_rlx_sclr; - hgt_fctr_sclr = exp(- k * dz / z_0) / z_0; - } - */ - }; // 2d/3d children @@ -377,17 +312,20 @@ namespace setup } // function expecting a libmpdata++ solver as argument - void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, int rng_seed) + void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &p_e, int rng_seed) { blitz::secondIndex k; this->intcond_hlpr(solver, rhod, th_e, rv_e, rng_seed, k); + arr_1D_t p_d_e(p_e - detail::calc_p_v()(p_e, rv_e)); + arr_1D_t T(th_e * pow(p_d_e / 1.e5, R_d_over_c_pd())); + using ix = typename concurr_t::solver_t::ix; int nz = solver.advectee().extent(ix::w); real_t dz = (Z / si::metres) / (nz-1); int nx = solver.advectee().extent(0); real_t dx = (X / si::metres) / (nx-1); - solver.advectee(ix::rv) = prtrb_rv(th_e, rhod, dz)( + solver.advectee(ix::rv) = prtrb_rv(T, p_e, dz)( sqrt( pow(blitz::tensor::i * dx - (X / si::metres / 2.), 2) + pow(blitz::tensor::j * dz - (z_prtrb / si::metres), 2) @@ -395,18 +333,6 @@ namespace setup blitz::tensor::j * dz ); - // reinitialise th with dry conversion based on perturbed rv TODO: without loops, also in 3D - for (int i = 0; i < nx; ++i) - { - for(int k = 0; k < nz; ++k) - { - using libcloudphxx::common::theta_dry::std2dry; - using libcloudphxx::common::theta_dry::dry2std; - quantity si_rv_e = rv_e(k); - quantity si_rv = solver.advectee(ix::rv)(i, k); - solver.advectee(ix::th)(i, k) = std2dry(dry2std(th_e(k) * si::kelvins, si_rv_e), si_rv) / si::kelvins; - } - } } }; @@ -425,7 +351,7 @@ namespace setup } // function expecting a libmpdata++ solver as argument - void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, int rng_seed) + void intcond(concurr_t &solver, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &p_e, int rng_seed) { blitz::thirdIndex k; this->intcond_hlpr(solver, rhod, th_e, rv_e, rng_seed, k); @@ -448,14 +374,6 @@ namespace setup solver.advectee(ix::v) = 0; solver.vab_relaxed_state(1) = 0; } - -/* - // TODO: make it work in 3d? - MoistThermalGrabowskiClark99_3d() - { - throw std::runtime_error("Moist Thermal doesn't work in 3d yet"); - } -*/ }; }; }; diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index 75c6cc43f1..5c0dd5147d 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -6,16 +6,6 @@ #include #include "../git_revision.h" -namespace detail -{ - struct calc_p_v - { - setup::real_t operator()(setup::real_t p, setup::real_t rv) const - {return libcloudphxx::common::moist_air::p_v(p * si::pascals, rv) / si::pascals;} - BZ_DECLARE_FUNCTOR2(calc_p_v) - }; -}; - template class slvr_common : public slvr_dim { From 6f53fb36ce36fc3f8b0e10ce67a7476c684e1347 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Fri, 8 Jun 2018 14:17:56 +0200 Subject: [PATCH 32/74] use p_d profile in blk_1m: --- src/slvr_blk_1m.hpp | 26 ++++++++++++++++---------- src/slvr_lgrngn.hpp | 17 ++++++++++------- 2 files changed, 26 insertions(+), 17 deletions(-) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index e651b64eb8..1b2ee04d61 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -19,8 +19,8 @@ class slvr_blk_1m_common : public slvr_common using clock = typename parent_t::clock; private: - // a 2D/3D array with copy of the environmental total pressure of dry air - typename parent_t::arr_t &p_e; + // a 2D/3D array with copy of the environmental total pressure and env partial pressure of dry air + typename parent_t::arr_t &p_e, &p_d_e; void condevap() { @@ -31,15 +31,16 @@ class slvr_blk_1m_common : public slvr_common rr = this->state(ix::rr)(this->ijk); // rain water mixing ratio auto const rhod = (*this->mem->G)(this->ijk), + &p_d_e_arg = p_e(this->ijk), &p_e_arg = p_e(this->ijk); /* libcloudphxx::blk_1m::adj_cellwise_constp( - opts, rhod, p_e_arg, th, rv, rc, rr, this->dt + opts, rhod, p_e_arg, p_d_e_arg, th, rv, rc, rr, this->dt ); */ libcloudphxx::blk_1m::adj_cellwise_nwtrph( - opts, p_e_arg, th, rv, rc, this->dt + opts, p_e_arg, p_d_e_arg, th, rv, rc, this->dt ); this->mem->barrier(); } @@ -98,6 +99,9 @@ class slvr_blk_1m_common : public slvr_common // init the p_e array p_e = (*params.p_e)(this->vert_idx); + // init the p_d_e array + p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); + // deal with initial supersaturation condevap(); @@ -164,8 +168,9 @@ class slvr_blk_1m_common : public slvr_common rc = this->state(ix::rc)(this->ijk), rr = this->state(ix::rr)(this->ijk), rhod = (*this->mem->G)(this->ijk), - &p_e_arg = p_e(this->ijk); - libcloudphxx::blk_1m::rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, rhod, p_e_arg, th, rv, rc, rr); + &p_e_arg = p_e(this->ijk), + &p_d_e_arg = p_d_e(this->ijk); + libcloudphxx::blk_1m::rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, rhod, p_e_arg, p_d_e_arg, th, rv, rc, rr); } // forcing @@ -202,8 +207,8 @@ class slvr_blk_1m_common : public slvr_common this->mem->barrier(); if(this->rank == 0) { - nancheck(rhs.at(ix::rc)(this->domain), "RHS of rc after rhs_update"); - nancheck(rhs.at(ix::rr)(this->domain), "RHS of rr after rhs_update"); + nancheck2(rhs.at(ix::rc)(this->domain), this->state(ix::rc)(this->domain), "RHS of rc after rhs_update"); + nancheck2(rhs.at(ix::rr)(this->domain), this->state(ix::rr)(this->domain), "RHS of rr after rhs_update"); this->tend = clock::now(); this->tupdate += std::chrono::duration_cast( this->tend - this->tbeg ); } @@ -236,7 +241,7 @@ class slvr_blk_1m_common : public slvr_common static void alloc(typename parent_t::mem_t *mem, const int &n_iters) { parent_t::alloc(mem, n_iters); - parent_t::alloc_tmp_sclr(mem, __FILE__, 1); // p_e + parent_t::alloc_tmp_sclr(mem, __FILE__, 2); // p_e, p_d_e } // ctor @@ -248,7 +253,8 @@ class slvr_blk_1m_common : public slvr_common params(p), opts(p.cloudph_opts), puddle(0), - p_e(args.mem->tmp[__FILE__][0][0]) + p_e(args.mem->tmp[__FILE__][0][0]), + p_d_e(args.mem->tmp[__FILE__][0][1]) {} }; diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index d321217225..68e719ad55 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -285,10 +285,12 @@ class slvr_lgrngn : public slvr_common params.cloudph_opts_init )); + // temporary array of densities - prtcls cant be init'd with 1D profile typename parent_t::arr_t rhod(this->mem->advectee(ix::th).shape()); // TODO: replace all rhod arrays with this->mem->G rhod = (*params.rhod)(this->vert_idx); + // TODO: below arrays of pressures are constant, init them and store them once as it is done in slvr_blk_1m // temporary array of pressure - prtcls cant be init'd with 1D profile typename parent_t::arr_t p_e(this->mem->advectee(ix::th).shape()); p_e = (*params.p_e)(this->vert_idx); @@ -296,13 +298,14 @@ class slvr_lgrngn : public slvr_common // temporary array of partial pressure of dry air - prtcls cant be init'd with 1D profile typename parent_t::arr_t p_d_e(this->mem->advectee(ix::th).shape()); p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); - prtcls->init( - make_arrinfo(this->mem->advectee(ix::th)), - make_arrinfo(this->mem->advectee(ix::rv)), - make_arrinfo(rhod) - ,make_arrinfo(p_e) - ,make_arrinfo(p_d_e) - ); + + prtcls->init( + make_arrinfo(this->mem->advectee(ix::th)), + make_arrinfo(this->mem->advectee(ix::rv)), + make_arrinfo(rhod) + ,make_arrinfo(p_e) + ,make_arrinfo(p_d_e) + ); // writing diagnostic data for the initial condition parent_t::tbeg_loop = parent_t::clock::now(); From 643d73207d479836c4cf7ce54212b216ea80be39 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 11 Jun 2018 13:12:33 +0200 Subject: [PATCH 33/74] pass p_e to intcond also in 3d --- src/bicycles.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/bicycles.cpp b/src/bicycles.cpp index 55b005403c..421ccb8ce5 100644 --- a/src/bicycles.cpp +++ b/src/bicycles.cpp @@ -237,17 +237,17 @@ void run(int nx, int ny, int nz, const user_params_t &user_params) if(user_params.model_case == "dry_thermal") { concurr.reset(new concurr_openmp_cyclic_t(p)); - case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, user_params.rng_seed); // works only by chance? + case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, p_e, user_params.rng_seed); // works only by chance? } else if(user_params.model_case == "lasher_trapp") { concurr.reset(new concurr_openmp_rigid_rigid_t(p)); - case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, user_params.rng_seed); // works only by chance? + case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, p_e, user_params.rng_seed); // works only by chance? } else { concurr.reset(new concurr_openmp_rigid_t(p)); - case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, user_params.rng_seed); + case_ptr->intcond(*static_cast(concurr.get()), rhod, th_e, rv_e, p_e, user_params.rng_seed); } // setup panic pointer and the signal handler From 6870c8d2165de4857eddae23d5296791860ec0ca Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Tue, 12 Jun 2018 14:44:54 +0200 Subject: [PATCH 34/74] blk_1m: fix p_d value --- src/slvr_blk_1m.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index 1b2ee04d61..72123465fa 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -31,10 +31,13 @@ class slvr_blk_1m_common : public slvr_common rr = this->state(ix::rr)(this->ijk); // rain water mixing ratio auto const rhod = (*this->mem->G)(this->ijk), - &p_d_e_arg = p_e(this->ijk), + &p_d_e_arg = p_d_e(this->ijk), &p_e_arg = p_e(this->ijk); /* + libcloudphxx::blk_1m::adj_cellwise( + opts, rhod, th, rv, rc, rr, this->dt + ); libcloudphxx::blk_1m::adj_cellwise_constp( opts, rhod, p_e_arg, p_d_e_arg, th, rv, rc, rr, this->dt ); @@ -171,6 +174,7 @@ class slvr_blk_1m_common : public slvr_common &p_e_arg = p_e(this->ijk), &p_d_e_arg = p_d_e(this->ijk); libcloudphxx::blk_1m::rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, rhod, p_e_arg, p_d_e_arg, th, rv, rc, rr); + //libcloudphxx::blk_1m::rhs_cellwise(opts, dot_rc, dot_rr, rc, rr); } // forcing From 01404220703ea7b2cafa6e4defdb2fb3c806204d Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Tue, 12 Jun 2018 17:04:23 +0200 Subject: [PATCH 35/74] timing of nondelayed step --- src/slvr_common.hpp | 14 ++++++++++++-- src/slvr_lgrngn.hpp | 6 +++++- 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index 5c0dd5147d..7efd1747e9 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -22,8 +22,8 @@ class slvr_common : public slvr_dim // TODO: timing slows down simulations // either remove it and use profiling tools (e.g. vtune) // or add some compile-time flag to turn it off - clock::time_point tbeg, tend, tbeg1, tbeg_loop; - std::chrono::milliseconds tdiag, tupdate, tsync, tasync, tasync_wait, tloop, tvip_rhs; + clock::time_point tbeg, tend, tbeg_loop; + std::chrono::milliseconds tdiag, tupdate, tsync, tasync, tasync_wait, tloop, tvip_rhs, tnondelayed_step; int spinup; // number of timesteps @@ -116,6 +116,15 @@ class slvr_common : public slvr_dim parent_t::hook_ante_step(); } + void hook_ante_delayed_step() + { + if(this->rank == 0) + { + tend = clock::now(); + tnondelayed_step += std::chrono::duration_cast( tend - tbeg ); + } + } + // get shape from a rng_t or an idx_t inline int shape(const rng_t &rng) { return rng.length();} @@ -324,6 +333,7 @@ class slvr_common : public slvr_dim << "custom vip_rhs: " << tvip_rhs.count() << " ("<< setup::real_t(tvip_rhs.count())/tloop.count()*100 <<"%)" << std::endl << "diag: " << tdiag.count() << " ("<< setup::real_t(tdiag.count())/tloop.count()*100 <<"%)" << std::endl << "sync: " << tsync.count() << " ("<< setup::real_t(tsync.count())/tloop.count()*100 <<"%)" << std::endl + << "nondelayed step: " << tnondelayed_step.count() << " ("<< setup::real_t(tnondelayed_step.count())/tloop.count()*100 <<"%)" << std::endl << "async: " << tasync.count() << " ("<< setup::real_t(tasync.count())/tloop.count()*100 <<"%)" << std::endl << "async_wait: " << tasync_wait.count() << " ("<< setup::real_t(tasync_wait.count())/tloop.count()*100 <<"%)" << std::endl; } diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 68e719ad55..c6a314ca8e 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -383,7 +383,6 @@ class slvr_lgrngn : public slvr_common this->mem->barrier(); if (this->rank == 0) { - parent_t::tbeg1 = parent_t::clock::now(); // assuring previous async step finished ... #if defined(STD_FUTURE_WORKS) if ( @@ -489,6 +488,11 @@ class slvr_lgrngn : public slvr_common } } this->mem->barrier(); + // start recording time of the non-delayed advection step + if (this->rank == 0) + { + parent_t::tbeg = parent_t::clock::now(); + } } void record_all() From 98fcbcad9e2d7b2c048f4758ef7f2878ba774fb4 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 14 Jun 2018 17:16:34 +0200 Subject: [PATCH 36/74] fix slvr_blk_1m pressures --- src/slvr_blk_1m.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index 72123465fa..ca72432f89 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -100,10 +100,10 @@ class slvr_blk_1m_common : public slvr_common zero_if_uninitialised(ix::rr); // init the p_e array - p_e = (*params.p_e)(this->vert_idx); + p_e(this->ijk).reindex(this->zero) = (*params.p_e)(this->vert_idx); // init the p_d_e array - p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); + p_d_e(this->ijk).reindex(this->zero) = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); // deal with initial supersaturation condevap(); From e39477572baab92b0244aac2cc945123ce6b2c67 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 18 Jun 2018 12:34:04 +0200 Subject: [PATCH 37/74] delay advection of th and rv, async calc of step_sync during advection of velocities --- src/bicycles.cpp | 2 + src/slvr_common.hpp | 14 ++---- src/slvr_lgrngn.hpp | 110 ++++++++++++++++++++++++++++++++++++-------- 3 files changed, 97 insertions(+), 29 deletions(-) diff --git a/src/bicycles.cpp b/src/bicycles.cpp index 421ccb8ce5..b641df996f 100644 --- a/src/bicycles.cpp +++ b/src/bicycles.cpp @@ -340,6 +340,7 @@ struct ct_params_2D_sd : ct_params_common u, w, th, rv, vip_i=u, vip_j=w, vip_den=-1 }; }; + enum { delayed_step = opts::bit(ix::th) | opts::bit(ix::rv) }; }; struct ct_params_3D_sd : ct_params_common @@ -350,6 +351,7 @@ struct ct_params_3D_sd : ct_params_common u, v, w, th, rv, vip_i=u, vip_j=v, vip_k=w, vip_den=-1 }; }; + enum { delayed_step = opts::bit(ix::th) | opts::bit(ix::rv) }; }; struct ct_params_2D_blk_1m : ct_params_common diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index 7efd1747e9..f690073ce9 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -23,7 +23,7 @@ class slvr_common : public slvr_dim // either remove it and use profiling tools (e.g. vtune) // or add some compile-time flag to turn it off clock::time_point tbeg, tend, tbeg_loop; - std::chrono::milliseconds tdiag, tupdate, tsync, tasync, tasync_wait, tloop, tvip_rhs, tnondelayed_step; + std::chrono::milliseconds tdiag, tupdate, tsync, tsync_wait, tasync, tasync_wait, tloop, tvip_rhs, tnondelayed_step; int spinup; // number of timesteps @@ -116,15 +116,6 @@ class slvr_common : public slvr_dim parent_t::hook_ante_step(); } - void hook_ante_delayed_step() - { - if(this->rank == 0) - { - tend = clock::now(); - tnondelayed_step += std::chrono::duration_cast( tend - tbeg ); - } - } - // get shape from a rng_t or an idx_t inline int shape(const rng_t &rng) { return rng.length();} @@ -335,7 +326,8 @@ class slvr_common : public slvr_dim << "sync: " << tsync.count() << " ("<< setup::real_t(tsync.count())/tloop.count()*100 <<"%)" << std::endl << "nondelayed step: " << tnondelayed_step.count() << " ("<< setup::real_t(tnondelayed_step.count())/tloop.count()*100 <<"%)" << std::endl << "async: " << tasync.count() << " ("<< setup::real_t(tasync.count())/tloop.count()*100 <<"%)" << std::endl - << "async_wait: " << tasync_wait.count() << " ("<< setup::real_t(tasync_wait.count())/tloop.count()*100 <<"%)" << std::endl; + << "async_wait: " << tasync_wait.count() << " ("<< setup::real_t(tasync_wait.count())/tloop.count()*100 <<"%)" << std::endl + << "sync_wait: " << tsync_wait.count() << " ("<< setup::real_t(tsync_wait.count())/tloop.count()*100 <<"%)" << std::endl; } } } diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index c6a314ca8e..c0b6b5ad43 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -434,27 +434,106 @@ class slvr_lgrngn : public slvr_common nancheck(this->mem->advectee(ix::rv), "rv before step sync"); negcheck(this->mem->advectee(ix::th), "th before step sync"); negcheck(this->mem->advectee(ix::rv), "rv before step sync"); - prtcls->step_sync( - params.cloudph_opts, - make_arrinfo(this->mem->advectee(ix::th)), - make_arrinfo(this->mem->advectee(ix::rv)), - libcloudphxx::lgrngn::arrinfo_t(), - make_arrinfo(Cx), - this->n_dims == 2 ? libcloudphxx::lgrngn::arrinfo_t() : make_arrinfo(Cy), - make_arrinfo(Cz) - ); - + + using libcloudphxx::lgrngn::particles_t; + using libcloudphxx::lgrngn::CUDA; + using libcloudphxx::lgrngn::multi_CUDA; + +#if defined(STD_FUTURE_WORKS) + if (params.async) + { + + prtcls->sync_in( + make_arrinfo(this->mem->advectee(ix::th)), + make_arrinfo(this->mem->advectee(ix::rv)), + libcloudphxx::lgrngn::arrinfo_t(), + make_arrinfo(Cx), + this->n_dims == 2 ? libcloudphxx::lgrngn::arrinfo_t() : make_arrinfo(Cy), + make_arrinfo(Cz) + ); + + assert(!ftr.valid()); + if(params.backend == CUDA) + ftr = std::async( + std::launch::async, + &particles_t::step_cond, + dynamic_cast*>(prtcls.get()), + params.cloudph_opts, + make_arrinfo(this->mem->advectee(ix::th)), + make_arrinfo(this->mem->advectee(ix::rv)), + std::map >() + ); + else if(params.backend == multi_CUDA) + ftr = std::async( + std::launch::async, + &particles_t::step_cond, + dynamic_cast*>(prtcls.get()), + params.cloudph_opts, + make_arrinfo(this->mem->advectee(ix::th)), + make_arrinfo(this->mem->advectee(ix::rv)), + std::map >() + ); + assert(ftr.valid()); + } else +#endif + { + prtcls->step_sync( + params.cloudph_opts, + make_arrinfo(this->mem->advectee(ix::th)), + make_arrinfo(this->mem->advectee(ix::rv)), + libcloudphxx::lgrngn::arrinfo_t(), + make_arrinfo(Cx), + this->n_dims == 2 ? libcloudphxx::lgrngn::arrinfo_t() : make_arrinfo(Cy), + make_arrinfo(Cz) + ); + // microphysics could have led to rv < 0 ? + negtozero(this->mem->advectee(ix::rv), "rv after step sync"); + + nancheck(this->mem->advectee(ix::th), "th after step sync"); + nancheck(this->mem->advectee(ix::rv), "rv after step sync"); + negcheck(this->mem->advectee(ix::th), "th after step sync"); + negcheck(this->mem->advectee(ix::rv), "rv after step sync"); + } + + parent_t::tend = parent_t::clock::now(); + parent_t::tsync += std::chrono::duration_cast( parent_t::tend - parent_t::tbeg ); + } + // start recording time of the non-delayed advection step + parent_t::tbeg = parent_t::clock::now(); + } + this->mem->barrier(); + } + + + void hook_ante_delayed_step() + { + parent_t::hook_ante_delayed_step(); + if (this->rank == 0) + { + parent_t::tend = parent_t::clock::now(); + parent_t::tnondelayed_step += std::chrono::duration_cast( parent_t::tend - parent_t::tbeg ); + + // assuring previous sync step finished ... +#if defined(STD_FUTURE_WORKS) + if ( + params.async + ) { + assert(ftr.valid()); + parent_t::tbeg = parent_t::clock::now(); + ftr.get(); // microphysics could have led to rv < 0 ? + // TODO: repeated above, unify negtozero(this->mem->advectee(ix::rv), "rv after step sync"); nancheck(this->mem->advectee(ix::th), "th after step sync"); nancheck(this->mem->advectee(ix::rv), "rv after step sync"); negcheck(this->mem->advectee(ix::th), "th after step sync"); negcheck(this->mem->advectee(ix::rv), "rv after step sync"); - parent_t::tend = parent_t::clock::now(); - parent_t::tsync += std::chrono::duration_cast( parent_t::tend - parent_t::tbeg ); - } + parent_t::tsync_wait += std::chrono::duration_cast( parent_t::tend - parent_t::tbeg ); + } else assert(!ftr.valid()); +#endif + // running asynchronous stuff { using libcloudphxx::lgrngn::particles_t; @@ -488,11 +567,6 @@ class slvr_lgrngn : public slvr_common } } this->mem->barrier(); - // start recording time of the non-delayed advection step - if (this->rank == 0) - { - parent_t::tbeg = parent_t::clock::now(); - } } void record_all() From 3df4499241100bcd8c4c96f8c9b45b3d1a5c19cc Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 14 Jun 2018 17:16:34 +0200 Subject: [PATCH 38/74] fix slvr_blk_1m pressures --- src/slvr_blk_1m.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index 72123465fa..ca72432f89 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -100,10 +100,10 @@ class slvr_blk_1m_common : public slvr_common zero_if_uninitialised(ix::rr); // init the p_e array - p_e = (*params.p_e)(this->vert_idx); + p_e(this->ijk).reindex(this->zero) = (*params.p_e)(this->vert_idx); // init the p_d_e array - p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); + p_d_e(this->ijk).reindex(this->zero) = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); // deal with initial supersaturation condevap(); From 8453f720aab7715211823abb77705daf9ce7b2f2 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 10 Jul 2018 03:07:10 +0200 Subject: [PATCH 39/74] Revert "pass p_d_e to lgrngn micro init" This reverts commit e733c5ce8737f84a296fc3286bc6098e852b687b. Conflicts: src/slvr_lgrngn.hpp --- src/slvr_lgrngn.hpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 68e719ad55..7506e78472 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -294,17 +294,12 @@ class slvr_lgrngn : public slvr_common // temporary array of pressure - prtcls cant be init'd with 1D profile typename parent_t::arr_t p_e(this->mem->advectee(ix::th).shape()); p_e = (*params.p_e)(this->vert_idx); - - // temporary array of partial pressure of dry air - prtcls cant be init'd with 1D profile - typename parent_t::arr_t p_d_e(this->mem->advectee(ix::th).shape()); - p_d_e = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); prtcls->init( make_arrinfo(this->mem->advectee(ix::th)), make_arrinfo(this->mem->advectee(ix::rv)), make_arrinfo(rhod) ,make_arrinfo(p_e) - ,make_arrinfo(p_d_e) ); // writing diagnostic data for the initial condition From cbfe44b7c9a2a285bc52e3c683e03bcc4b8e2820 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 10 Jul 2018 03:10:21 +0200 Subject: [PATCH 40/74] Revert "use p_d profile in blk_1m:" This reverts commit 6f53fb36ce36fc3f8b0e10ce67a7476c684e1347. Conflicts: src/slvr_blk_1m.hpp src/slvr_lgrngn.hpp --- src/slvr_blk_1m.hpp | 27 ++++++++++----------------- src/slvr_lgrngn.hpp | 2 -- 2 files changed, 10 insertions(+), 19 deletions(-) diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index ca72432f89..49ef5273cd 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -19,8 +19,8 @@ class slvr_blk_1m_common : public slvr_common using clock = typename parent_t::clock; private: - // a 2D/3D array with copy of the environmental total pressure and env partial pressure of dry air - typename parent_t::arr_t &p_e, &p_d_e; + // a 2D/3D array with copy of the environmental total pressure of dry air + typename parent_t::arr_t &p_e; void condevap() { @@ -31,7 +31,6 @@ class slvr_blk_1m_common : public slvr_common rr = this->state(ix::rr)(this->ijk); // rain water mixing ratio auto const rhod = (*this->mem->G)(this->ijk), - &p_d_e_arg = p_d_e(this->ijk), &p_e_arg = p_e(this->ijk); /* @@ -39,11 +38,11 @@ class slvr_blk_1m_common : public slvr_common opts, rhod, th, rv, rc, rr, this->dt ); libcloudphxx::blk_1m::adj_cellwise_constp( - opts, rhod, p_e_arg, p_d_e_arg, th, rv, rc, rr, this->dt + opts, rhod, p_e_arg, th, rv, rc, rr, this->dt ); */ libcloudphxx::blk_1m::adj_cellwise_nwtrph( - opts, p_e_arg, p_d_e_arg, th, rv, rc, this->dt + opts, p_e_arg, th, rv, rc, this->dt ); this->mem->barrier(); } @@ -102,9 +101,6 @@ class slvr_blk_1m_common : public slvr_common // init the p_e array p_e(this->ijk).reindex(this->zero) = (*params.p_e)(this->vert_idx); - // init the p_d_e array - p_d_e(this->ijk).reindex(this->zero) = (*params.p_e)(this->vert_idx) - detail::calc_p_v()((*params.p_e)(this->vert_idx), (*params.rv_e)(this->vert_idx)); - // deal with initial supersaturation condevap(); @@ -171,10 +167,8 @@ class slvr_blk_1m_common : public slvr_common rc = this->state(ix::rc)(this->ijk), rr = this->state(ix::rr)(this->ijk), rhod = (*this->mem->G)(this->ijk), - &p_e_arg = p_e(this->ijk), - &p_d_e_arg = p_d_e(this->ijk); - libcloudphxx::blk_1m::rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, rhod, p_e_arg, p_d_e_arg, th, rv, rc, rr); - //libcloudphxx::blk_1m::rhs_cellwise(opts, dot_rc, dot_rr, rc, rr); + &p_e_arg = p_e(this->ijk); + libcloudphxx::blk_1m::rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, rhod, p_e_arg, th, rv, rc, rr); } // forcing @@ -211,8 +205,8 @@ class slvr_blk_1m_common : public slvr_common this->mem->barrier(); if(this->rank == 0) { - nancheck2(rhs.at(ix::rc)(this->domain), this->state(ix::rc)(this->domain), "RHS of rc after rhs_update"); - nancheck2(rhs.at(ix::rr)(this->domain), this->state(ix::rr)(this->domain), "RHS of rr after rhs_update"); + nancheck(rhs.at(ix::rc)(this->domain), "RHS of rc after rhs_update"); + nancheck(rhs.at(ix::rr)(this->domain), "RHS of rr after rhs_update"); this->tend = clock::now(); this->tupdate += std::chrono::duration_cast( this->tend - this->tbeg ); } @@ -245,7 +239,7 @@ class slvr_blk_1m_common : public slvr_common static void alloc(typename parent_t::mem_t *mem, const int &n_iters) { parent_t::alloc(mem, n_iters); - parent_t::alloc_tmp_sclr(mem, __FILE__, 2); // p_e, p_d_e + parent_t::alloc_tmp_sclr(mem, __FILE__, 1); // p_e } // ctor @@ -257,8 +251,7 @@ class slvr_blk_1m_common : public slvr_common params(p), opts(p.cloudph_opts), puddle(0), - p_e(args.mem->tmp[__FILE__][0][0]), - p_d_e(args.mem->tmp[__FILE__][0][1]) + p_e(args.mem->tmp[__FILE__][0][0]) {} }; diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 7506e78472..d1295d6687 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -285,12 +285,10 @@ class slvr_lgrngn : public slvr_common params.cloudph_opts_init )); - // temporary array of densities - prtcls cant be init'd with 1D profile typename parent_t::arr_t rhod(this->mem->advectee(ix::th).shape()); // TODO: replace all rhod arrays with this->mem->G rhod = (*params.rhod)(this->vert_idx); - // TODO: below arrays of pressures are constant, init them and store them once as it is done in slvr_blk_1m // temporary array of pressure - prtcls cant be init'd with 1D profile typename parent_t::arr_t p_e(this->mem->advectee(ix::th).shape()); p_e = (*params.p_e)(this->vert_idx); From b3bed3dcaa8a6bdfeb7465c0068db4587176604e Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 11 Jul 2018 01:05:20 +0200 Subject: [PATCH 41/74] DYCOMS: revert to reference state dry density and potential temp at z=0 equal to env state's total density and standard potential temp --- src/cases/DYCOMS98.hpp | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/cases/DYCOMS98.hpp b/src/cases/DYCOMS98.hpp index e51a1e2110..3f1f250411 100644 --- a/src/cases/DYCOMS98.hpp +++ b/src/cases/DYCOMS98.hpp @@ -238,9 +238,9 @@ std::cout << "lwp env: " << lwp_env << std::endl; st(notopbot) = (th_e(notopbot+1) - th_e(notopbot-1)) / th_e(notopbot); real_t st_avg = blitz::sum(st) / (nz-2) / (2.*dz); // reference theta - // th_ref = th_e(0) * exp(st_avg * k * dz); - th_ref = th_e(0) * pow(1 + rv_e(0) / a, f) // calc dry theta at z=0 - * exp(st_avg * k * dz); + th_ref = th_e(0) * exp(st_avg * k * dz); + // th_ref = th_e(0) * pow(1 + rv_e(0) / a, f) // calc dry theta at z=0 + // * exp(st_avg * k * dz); // virtual temp at surface using libcloudphxx::common::moist_air::R_d_over_c_pd; using libcloudphxx::common::moist_air::c_pd; @@ -248,22 +248,23 @@ std::cout << "lwp env: " << lwp_env << std::endl; using libcloudphxx::common::theta_std::p_1000; real_t T_surf = th_e(0) * pow(p_0 / p_1000(), R_d_over_c_pd()); -/* + real_t T_virt_surf = T_surf * (1. + 0.608 * rv_e(0)); real_t rho_surf = (p_0 / si::pascals) / T_virt_surf / 287. ; // TODO: R_d instead of 287, its the total, not dry density! - rho_surf /= (1 + rv_e(0)); // turn it into dry air density! TODO: is this correct? TODO2: approp change in the paper -*/ - real_t rho_surf = (p_0 / si::pascals) / T_surf / (1. + 29. / 18. * rv_e(0)) / 287. ; // dry air density at the surface TODO: R_d instead of 287 +// rho_surf /= (1 + rv_e(0)); // turn it into dry air density! TODO: is this correct? TODO2: approp change in the paper - // theta_std env prof to theta_dry_e - for(int k=1; k(th_e(k) * si::kelvins, quantity(rv_e(k))) / si::kelvins; + // real_t rho_surf = (p_0 / si::pascals) / T_surf / (1. + 29. / 18. * rv_e(0)) / 287. ; // dry air density at the surface TODO: R_d instead of 287 // real_t cs = 9.81 / (c_pd() / si::joules * si::kilograms * si::kelvins) / st_avg / T_surf; // this is from Wojtek real_t cs = 9.81 / (c_pd() / si::joules * si::kilograms * si::kelvins) / st_avg / th_e(0); // this is correct? or total, not dry th_e(0) should be here? // rhod profile rhod = rho_surf * exp(- st_avg * k * dz) * pow( 1. - cs * (1 - exp(- st_avg * k * dz)), (1. / R_d_over_c_pd()) - 1); + + + // theta_std env prof to theta_dry_e + for(int k=1; k(th_e(k) * si::kelvins, quantity(rv_e(k))) / si::kelvins; // subsidence rate w_LS = w_LS_fctr()(k * dz); From faceda31121b33c2e9e4527b79d5165ded00025c Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 18 Jul 2018 12:27:18 +0200 Subject: [PATCH 42/74] record RH_formula used --- src/slvr_lgrngn.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 302abca2ed..5de21d549a 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -342,6 +342,7 @@ class slvr_lgrngn : public slvr_common this->record_aux_const("rng_seed", params.cloudph_opts_init.rng_seed); this->record_aux_const("kernel", params.cloudph_opts_init.kernel); this->record_aux_const("terminal_velocity", params.cloudph_opts_init.terminal_velocity); + this->record_aux_const("RH_formula", params.cloudph_opts_init.RH_formula); this->record_aux_const("backend", params.backend); this->record_aux_const("async", params.async); this->record_aux_const("adve", params.cloudph_opts.adve); From 928fc59d0da66dd0bb072c87bb01b4c2756642e2 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 19 Jul 2018 14:13:43 +0200 Subject: [PATCH 43/74] sd_conc_large_tail option --- src/opts_lgrngn.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/opts_lgrngn.hpp b/src/opts_lgrngn.hpp index a16f020c5a..9ae6f986ff 100644 --- a/src/opts_lgrngn.hpp +++ b/src/opts_lgrngn.hpp @@ -53,6 +53,7 @@ void setopts_micro( ("dev_id", po::value()->default_value(-1), "CUDA backend - id of device to be used") // free parameters ("exact_sstp_cond", po::value()->default_value(rt_params.cloudph_opts_init.exact_sstp_cond), "exact(per-particle) logic for substeps for condensation") + ("sd_conc_large_tail", po::value()->default_value(rt_params.cloudph_opts_init.sd_conc_large_tail), "add SDs to better represent the large tail") ("sstp_cond", po::value()->default_value(rt_params.cloudph_opts_init.sstp_cond), "no. of substeps for condensation") ("sstp_coal", po::value()->default_value(rt_params.cloudph_opts_init.sstp_coal), "no. of substeps for coalescence") ("sstp_chem", po::value()->default_value(rt_params.cloudph_opts_init.sstp_chem), "no. of substeps for chemistry") @@ -146,6 +147,7 @@ void setopts_micro( rt_params.cloudph_opts_init.sstp_coal = vm["sstp_coal"].as(); rt_params.cloudph_opts_init.sstp_chem = vm["sstp_chem"].as(); rt_params.cloudph_opts_init.exact_sstp_cond = vm["exact_sstp_cond"].as(); + rt_params.cloudph_opts_init.sd_conc_large_tail = vm["sd_conc_large_tail"].as(); rt_params.cloudph_opts_init.rng_seed = user_params.rng_seed; From e40b2ae77c42b7232c18891e50db2015dece9a07 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 25 Jul 2018 15:24:07 +0200 Subject: [PATCH 44/74] mixed explicit/implicit integration; breaks blk_1m --- src/bicycles.cpp | 2 +- src/calc_forces_common.hpp | 17 ++++++++++------- src/slvr_common.hpp | 25 ++++++++++++++----------- src/slvr_lgrngn.hpp | 17 +++++++++++++++++ 4 files changed, 42 insertions(+), 19 deletions(-) diff --git a/src/bicycles.cpp b/src/bicycles.cpp index b641df996f..2e555c14fd 100644 --- a/src/bicycles.cpp +++ b/src/bicycles.cpp @@ -327,7 +327,7 @@ void run(int nx, int ny, int nz, const user_params_t &user_params) struct ct_params_common : ct_params_default_t { using real_t = setup::real_t; - enum { rhs_scheme = solvers::trapez }; + enum { rhs_scheme = solvers::mixed }; enum { prs_scheme = solvers::cr }; enum { vip_vab = solvers::expl }; }; diff --git a/src/calc_forces_common.hpp b/src/calc_forces_common.hpp index bc73995799..9f234edd73 100644 --- a/src/calc_forces_common.hpp +++ b/src/calc_forces_common.hpp @@ -134,16 +134,19 @@ void slvr_common::th_src(typename parent_t::arr_t &rv) } template -void slvr_common::w_src(typename parent_t::arr_t &th, typename parent_t::arr_t &rv) +void slvr_common::w_src(typename parent_t::arr_t &th, typename parent_t::arr_t &rv, const int at) { const auto &ijk = this->ijk; // buoyancy buoyancy(th, rv); - alpha(ijk) = F(ijk); - // large-scale vertical wind - subsidence(ix::w); // TODO: in case 1, w here should be in step n+1, calc it explicitly as w + 0.5 * dt * rhs(w); - // could also be calculated implicitly, but we would need implicit w^n+1 in other cells; - // also include absorber in w^n+1 estimate... + alpha(ijk) = 0.5 * F(ijk); // halved, because it is applied trapezoidaly + if(at == 0) // subsidence added explicitly, so updated only at n + { + // large-scale vertical wind + subsidence(ix::w); // TODO: in case 1, w here should be in step n+1, calc it explicitly as w + 0.5 * dt * rhs(w); + // could also be calculated implicitly, but we would need implicit w^n+1 in other cells; + // also include absorber in w^n+1 estimate... - alpha(ijk) += F(ijk); + alpha(ijk) += F(ijk); + } } diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index f690073ce9..5c2b525cdd 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -126,7 +126,7 @@ class slvr_common : public slvr_dim void radiation(typename parent_t::arr_t &rv); void rv_src(); void th_src(typename parent_t::arr_t &rv); - void w_src(typename parent_t::arr_t &th, typename parent_t::arr_t &rv); + void w_src(typename parent_t::arr_t &th, typename parent_t::arr_t &rv, const int at); void surf_sens(); void surf_latent(); void subsidence(const int&); @@ -137,7 +137,7 @@ class slvr_common : public slvr_dim const int &at ) { - parent_t::update_rhs(rhs, dt, at); + parent_t::update_rhs(rhs, dt, at); // zero-out rhs this->mem->barrier(); if(this->rank == 0) tbeg = clock::now(); @@ -164,7 +164,7 @@ class slvr_common : public slvr_dim // vertical velocity sources if(params.w_src && (!ct_params_t::piggy)) { - w_src(this->state(ix::th), this->state(ix::rv)); + w_src(this->state(ix::th), this->state(ix::rv), 0); rhs.at(ix_w)(ijk) += alpha(ijk); } @@ -185,17 +185,17 @@ class slvr_common : public slvr_dim // trapezoidal rhs^n+1 { // ---- water vapor sources ---- - rv_src(); - rhs.at(ix::rv)(ijk) += (alpha(ijk) + beta(ijk) * this->state(ix::rv)(ijk)) / (1. - 0.5 * this->dt * beta(ijk)); +// rv_src(); +// rhs.at(ix::rv)(ijk) += (alpha(ijk) + beta(ijk) * this->state(ix::rv)(ijk)) / (1. - 0.5 * this->dt * beta(ijk)); // TODO: alpha should also take (possibly impolicit) estimate of rv^n+1 too // becomes important when nudging is introduced? // ---- potential temp sources ---- - tmp2(ijk) = this->state(ix::rv)(ijk) + 0.5 * this->dt * rhs.at(ix::rv)(ijk); +// tmp2(ijk) = this->state(ix::rv)(ijk) + 0.5 * this->dt * rhs.at(ix::rv)(ijk); // todo: once rv_src beta!=0 (e.g. nudging), rv^n+1 estimate should be implicit here - th_src(tmp2); - rhs.at(ix::th)(ijk) += (alpha(ijk) + beta(ijk) * this->state(ix::th)(ijk)) / (1. - 0.5 * this->dt * beta(ijk)); +// th_src(tmp2); +// rhs.at(ix::th)(ijk) += (alpha(ijk) + beta(ijk) * this->state(ix::th)(ijk)) / (1. - 0.5 * this->dt * beta(ijk)); // TODO: alpha should also take (possibly impolicit) estimate of th^n+1 too // becomes important when nudging is introduced? @@ -203,19 +203,21 @@ class slvr_common : public slvr_dim if(params.w_src && (!ct_params_t::piggy)) { // temporarily use beta to store the th^n+1 estimate - beta(ijk) = this->state(ix::th)(ijk) + 0.5 * this->dt * rhs.at(ix::th)(ijk); +// beta(ijk) = this->state(ix::th)(ijk) + 0.5 * this->dt * rhs.at(ix::th)(ijk); // todo: oncethv_src beta!=0 (e.g. nudging), th^n+1 estimate should be implicit here // temporarily use F to store the rv^n+1 estimate - F(ijk) = this->state(ix::rv)(ijk) + 0.5 * this->dt * rhs.at(ix::rv)(ijk); + // F(ijk) = this->state(ix::rv)(ijk) + 0.5 * this->dt * rhs.at(ix::rv)(ijk); // todo: once rv_src beta!=0 (e.g. nudging), rv^n+1 estimate should be implicit here - w_src(beta, F); + // w_src(beta, F, 1); + w_src(this->state(ix::th)(ijk), this->state(ix::rv)(ijk), 1); rhs.at(ix::w)(ijk) += alpha(ijk); } // horizontal velocity sources // large-scale vertical wind +/* if(params.uv_src) { for(auto type : this->hori_vel) @@ -226,6 +228,7 @@ class slvr_common : public slvr_dim rhs.at(type)(ijk) += F(ijk); } } +*/ break; } } diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 5de21d549a..95537ca768 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -367,6 +367,9 @@ class slvr_lgrngn : public slvr_common // TODO: barrier? } + void hook_mixed_rhs_ante_loop() + {} // empty, because update_rhs called before each step; defined, to avoid assert from default hook_mixed_rhs_ante_loop + #if defined(STD_FUTURE_WORKS) std::future ftr; #endif @@ -499,6 +502,8 @@ class slvr_lgrngn : public slvr_common } + + void hook_ante_delayed_step() { parent_t::hook_ante_delayed_step(); @@ -562,6 +567,18 @@ class slvr_lgrngn : public slvr_common } this->mem->barrier(); } + + void hook_mixed_rhs_ante_step() + { + update_rhs(this->rhs, this->dt, this->n); + apply_rhs(this->dt); + } + + void hook_mixed_rhs_post_step() + { + update_rhs(this->rhs, this->dt, this->n+1); + apply_rhs(this->dt); + } void record_all() { From 4d1c2054941d0e680e54c937670bd0daae3c8b1b Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 25 Jul 2018 15:56:29 +0200 Subject: [PATCH 45/74] lgrngn: sync in before rhs_apply, cond after rhs_apply --- src/slvr_common.hpp | 2 +- src/slvr_lgrngn.hpp | 211 ++++++++++++++++++++++++-------------------- 2 files changed, 114 insertions(+), 99 deletions(-) diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index 5c2b525cdd..26bd12a658 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -211,7 +211,7 @@ class slvr_common : public slvr_dim // todo: once rv_src beta!=0 (e.g. nudging), rv^n+1 estimate should be implicit here // w_src(beta, F, 1); - w_src(this->state(ix::th)(ijk), this->state(ix::rv)(ijk), 1); + w_src(this->state(ix::th), this->state(ix::rv), 1); rhs.at(ix::w)(ijk) += alpha(ijk); } diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 95537ca768..070d9298bc 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -403,100 +403,6 @@ class slvr_lgrngn : public slvr_common nancheck(rl, "rl after copying from diag_wet_mom(3)"); rl = rl * 4./3. * 1000. * 3.14159; // get mixing ratio [kg/kg] - { - // temporarily Cx & Cz are multiplied by this->rhod ... - auto - Cx = this->mem->GC[0](this->Cx_domain).reindex(this->zero).copy(), - Cy = this->mem->GC[1](this->Cy_domain).reindex(this->zero).copy(), - Cz = this->mem->GC[ix::w](this->Cz_domain).reindex(this->zero).copy(); - nancheck(Cx, "Cx after copying from mpdata"); - nancheck(Cy, "Cy after copying from mpdata"); - nancheck(Cz, "Cz after copying from mpdata"); - - // ... and now dividing them by this->rhod (TODO: z=0 is located at k=1/2) - { - Cx.reindex(this->zero) /= (*params.rhod)(this->vert_idx); - Cy.reindex(this->zero) /= (*params.rhod)(this->vert_idx); - Cz.reindex(this->zero) /= (*params.rhod)(this->vert_idx); // TODO: should be interpolated, since theres a shift between positions of rhod and Cz - } - - // running synchronous stuff - parent_t::tbeg = parent_t::clock::now(); - - // rv might be negative due to large negative RHS from SD fluctuations + large-scale subsidence? - // turn all negative rv into rv = 0... CHEATING - negtozero(this->mem->advectee(ix::rv), "rv before step sync"); - - nancheck(this->mem->advectee(ix::th), "th before step sync"); - nancheck(this->mem->advectee(ix::rv), "rv before step sync"); - negcheck(this->mem->advectee(ix::th), "th before step sync"); - negcheck(this->mem->advectee(ix::rv), "rv before step sync"); - - using libcloudphxx::lgrngn::particles_t; - using libcloudphxx::lgrngn::CUDA; - using libcloudphxx::lgrngn::multi_CUDA; - -#if defined(STD_FUTURE_WORKS) - if (params.async) - { - - prtcls->sync_in( - make_arrinfo(this->mem->advectee(ix::th)), - make_arrinfo(this->mem->advectee(ix::rv)), - libcloudphxx::lgrngn::arrinfo_t(), - make_arrinfo(Cx), - this->n_dims == 2 ? libcloudphxx::lgrngn::arrinfo_t() : make_arrinfo(Cy), - make_arrinfo(Cz) - ); - - assert(!ftr.valid()); - if(params.backend == CUDA) - ftr = std::async( - std::launch::async, - &particles_t::step_cond, - dynamic_cast*>(prtcls.get()), - params.cloudph_opts, - make_arrinfo(this->mem->advectee(ix::th)), - make_arrinfo(this->mem->advectee(ix::rv)), - std::map >() - ); - else if(params.backend == multi_CUDA) - ftr = std::async( - std::launch::async, - &particles_t::step_cond, - dynamic_cast*>(prtcls.get()), - params.cloudph_opts, - make_arrinfo(this->mem->advectee(ix::th)), - make_arrinfo(this->mem->advectee(ix::rv)), - std::map >() - ); - assert(ftr.valid()); - } else -#endif - { - prtcls->step_sync( - params.cloudph_opts, - make_arrinfo(this->mem->advectee(ix::th)), - make_arrinfo(this->mem->advectee(ix::rv)), - libcloudphxx::lgrngn::arrinfo_t(), - make_arrinfo(Cx), - this->n_dims == 2 ? libcloudphxx::lgrngn::arrinfo_t() : make_arrinfo(Cy), - make_arrinfo(Cz) - ); - // microphysics could have led to rv < 0 ? - negtozero(this->mem->advectee(ix::rv), "rv after step sync"); - - nancheck(this->mem->advectee(ix::th), "th after step sync"); - nancheck(this->mem->advectee(ix::rv), "rv after step sync"); - negcheck(this->mem->advectee(ix::th), "th after step sync"); - negcheck(this->mem->advectee(ix::rv), "rv after step sync"); - } - - parent_t::tend = parent_t::clock::now(); - parent_t::tsync += std::chrono::duration_cast( parent_t::tend - parent_t::tbeg ); - } - // start recording time of the non-delayed advection step - parent_t::tbeg = parent_t::clock::now(); } this->mem->barrier(); } @@ -570,14 +476,123 @@ class slvr_lgrngn : public slvr_common void hook_mixed_rhs_ante_step() { - update_rhs(this->rhs, this->dt, this->n); - apply_rhs(this->dt); + this->update_rhs(this->rhs, this->dt, 0); + + // pass Eulerian fields to microphysics + if (this->rank == 0) + { + // temporarily Cx & Cz are multiplied by this->rhod ... + auto + Cx = this->mem->GC[0](this->Cx_domain).reindex(this->zero).copy(), + Cy = this->mem->GC[1](this->Cy_domain).reindex(this->zero).copy(), + Cz = this->mem->GC[ix::w](this->Cz_domain).reindex(this->zero).copy(); + nancheck(Cx, "Cx after copying from mpdata"); + nancheck(Cy, "Cy after copying from mpdata"); + nancheck(Cz, "Cz after copying from mpdata"); + + // ... and now dividing them by this->rhod (TODO: z=0 is located at k=1/2) + { + Cx.reindex(this->zero) /= (*params.rhod)(this->vert_idx); + Cy.reindex(this->zero) /= (*params.rhod)(this->vert_idx); + Cz.reindex(this->zero) /= (*params.rhod)(this->vert_idx); // TODO: should be interpolated, since theres a shift between positions of rhod and Cz + } + + // start synchronous stuff timer + parent_t::tbeg = parent_t::clock::now(); + + using libcloudphxx::lgrngn::particles_t; + using libcloudphxx::lgrngn::CUDA; + using libcloudphxx::lgrngn::multi_CUDA; + + prtcls->sync_in( + make_arrinfo(this->mem->advectee(ix::th)), + make_arrinfo(this->mem->advectee(ix::rv)), + libcloudphxx::lgrngn::arrinfo_t(), + make_arrinfo(Cx), + this->n_dims == 2 ? libcloudphxx::lgrngn::arrinfo_t() : make_arrinfo(Cy), + make_arrinfo(Cz) + ); + } + this->mem->barrier(); + + // apply rhs (except condensation) + this->apply_rhs(this->dt); + // rv might be negative due to large negative RHS from SD fluctuations + large-scale subsidence? + // turn all negative rv into rv = 0... CHEATING + negtozero(this->mem->advectee(ix::rv), "rv before step sync"); + nancheck(this->mem->advectee(ix::th), "th before step sync"); + nancheck(this->mem->advectee(ix::rv), "rv before step sync"); + negcheck(this->mem->advectee(ix::th), "th before step sync"); + negcheck(this->mem->advectee(ix::rv), "rv before step sync"); + + // start sync/async run of step_cond + // step_cond takes th and rv only for sync_out purposes - the values of th and rv before condensation come from sync_in, i.e. before apply_rhs + if (this->rank == 0) + { + using libcloudphxx::lgrngn::particles_t; + using libcloudphxx::lgrngn::CUDA; + using libcloudphxx::lgrngn::multi_CUDA; + +#if defined(STD_FUTURE_WORKS) + if (params.async) + { + assert(!ftr.valid()); + if(params.backend == CUDA) + ftr = std::async( + std::launch::async, + &particles_t::step_cond, + dynamic_cast*>(prtcls.get()), + params.cloudph_opts, + make_arrinfo(this->mem->advectee(ix::th)), + make_arrinfo(this->mem->advectee(ix::rv)), + std::map >() + ); + else if(params.backend == multi_CUDA) + ftr = std::async( + std::launch::async, + &particles_t::step_cond, + dynamic_cast*>(prtcls.get()), + params.cloudph_opts, + make_arrinfo(this->mem->advectee(ix::th)), + make_arrinfo(this->mem->advectee(ix::rv)), + std::map >() + ); + assert(ftr.valid()); + } else +#endif + { + prtcls->step_cond( + params.cloudph_opts, + make_arrinfo(this->mem->advectee(ix::th)), + make_arrinfo(this->mem->advectee(ix::rv)) + ); + // microphysics could have led to rv < 0 ? + negtozero(this->mem->advectee(ix::rv), "rv after step sync"); + + nancheck(this->mem->advectee(ix::th), "th after step sync"); + nancheck(this->mem->advectee(ix::rv), "rv after step sync"); + negcheck(this->mem->advectee(ix::th), "th after step sync"); + negcheck(this->mem->advectee(ix::rv), "rv after step sync"); + } + + parent_t::tend = parent_t::clock::now(); + parent_t::tsync += std::chrono::duration_cast( parent_t::tend - parent_t::tbeg ); + // start recording time of the non-delayed advection step + parent_t::tbeg = parent_t::clock::now(); + } + this->mem->barrier(); } void hook_mixed_rhs_post_step() { - update_rhs(this->rhs, this->dt, this->n+1); - apply_rhs(this->dt); + this->update_rhs(this->rhs, this->dt, 1); + this->apply_rhs(this->dt); + // no negtozero after apply, because only w changed here + // TODO: add these nanchecks/negchecks to apply_rhs, since they are repeated twice now + nancheck(this->mem->advectee(ix::th), "th after step sync"); + nancheck(this->mem->advectee(ix::rv), "rv after step sync"); + negcheck(this->mem->advectee(ix::th), "th after step sync"); + negcheck(this->mem->advectee(ix::rv), "rv after step sync"); } void record_all() From 9c40e58b755aa0882f24a0f4b534523091ea696d Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 26 Jul 2018 12:39:14 +0200 Subject: [PATCH 46/74] w8 for async before mixed_rhs_ante --- src/slvr_lgrngn.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 070d9298bc..9aaa2cc8d0 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -376,8 +376,6 @@ class slvr_lgrngn : public slvr_common void hook_ante_step() { - parent_t::hook_ante_step(); // includes output - this->mem->barrier(); if (this->rank == 0) { // assuring previous async step finished ... @@ -405,6 +403,7 @@ class slvr_lgrngn : public slvr_common } this->mem->barrier(); + parent_t::hook_ante_step(); // includes output } @@ -525,6 +524,7 @@ class slvr_lgrngn : public slvr_common negcheck(this->mem->advectee(ix::th), "th before step sync"); negcheck(this->mem->advectee(ix::rv), "rv before step sync"); + this->mem->barrier(); // start sync/async run of step_cond // step_cond takes th and rv only for sync_out purposes - the values of th and rv before condensation come from sync_in, i.e. before apply_rhs if (this->rank == 0) From 349f0eb867de301fb2bb663d1bde57e58014fcad Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 27 Jul 2018 14:47:19 +0200 Subject: [PATCH 47/74] add a todo in buoyancy --- src/calc_forces_common.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/calc_forces_common.hpp b/src/calc_forces_common.hpp index 9f234edd73..0cf4ebec12 100644 --- a/src/calc_forces_common.hpp +++ b/src/calc_forces_common.hpp @@ -138,6 +138,7 @@ void slvr_common::w_src(typename parent_t::arr_t &th, typename pare { const auto &ijk = this->ijk; // buoyancy + // TODO: buoyancy is now calculated twice, at n and at n+1, make it so that it is calculated once (will need to remove zeroing-out of w rhs in parent:update-rhs) buoyancy(th, rv); alpha(ijk) = 0.5 * F(ijk); // halved, because it is applied trapezoidaly if(at == 0) // subsidence added explicitly, so updated only at n From c1da6c27ffcf3001e80c0eec7d33ab76bd3dd1e8 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 27 Jul 2018 14:49:12 +0200 Subject: [PATCH 48/74] mixed_rhs_ante_step: update rhs after sync_in for code readability --- src/slvr_lgrngn.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 9aaa2cc8d0..c9d52525dd 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -475,7 +475,6 @@ class slvr_lgrngn : public slvr_common void hook_mixed_rhs_ante_step() { - this->update_rhs(this->rhs, this->dt, 0); // pass Eulerian fields to microphysics if (this->rank == 0) @@ -514,8 +513,9 @@ class slvr_lgrngn : public slvr_common } this->mem->barrier(); - // apply rhs (except condensation) + this->update_rhs(this->rhs, this->dt, 0); this->apply_rhs(this->dt); + // rv might be negative due to large negative RHS from SD fluctuations + large-scale subsidence? // turn all negative rv into rv = 0... CHEATING negtozero(this->mem->advectee(ix::rv), "rv before step sync"); From e9b1a4b76ea9bbb6ffe29854b644afccff59d34b Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 27 Jul 2018 15:13:34 +0200 Subject: [PATCH 49/74] todo on rhs:' --- src/slvr_lgrngn.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index c9d52525dd..6e5661f9b1 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -513,7 +513,7 @@ class slvr_lgrngn : public slvr_common } this->mem->barrier(); - this->update_rhs(this->rhs, this->dt, 0); + this->update_rhs(this->rhs, this->dt, 0); // TODO: update_rhs called twice per step causes halo filling twice (done by parent_t::update_rhs), probably not needed - we just need to set rhs to zero this->apply_rhs(this->dt); // rv might be negative due to large negative RHS from SD fluctuations + large-scale subsidence? From fa84db747b6c1a1d7df7c00353b035659683a033 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 27 Jul 2018 15:40:53 +0200 Subject: [PATCH 50/74] theta_std stuff: no std2dry and use rv_cc RH formula --- src/cases/DYCOMS98.hpp | 6 +++--- src/opts_lgrngn.hpp | 2 ++ 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/cases/DYCOMS98.hpp b/src/cases/DYCOMS98.hpp index 3f1f250411..462299662e 100644 --- a/src/cases/DYCOMS98.hpp +++ b/src/cases/DYCOMS98.hpp @@ -51,7 +51,7 @@ namespace setup { real_t operator()(const real_t &z) const { - return theta_dry::std2dry(th_l(z), r_t()(z)) / si::kelvins; + return th_l(z) / si::kelvins; } BZ_DECLARE_FUNCTOR(th_dry_fctr); }; @@ -263,8 +263,8 @@ std::cout << "lwp env: " << lwp_env << std::endl; // theta_std env prof to theta_dry_e - for(int k=1; k(th_e(k) * si::kelvins, quantity(rv_e(k))) / si::kelvins; +// for(int k=1; k(th_e(k) * si::kelvins, quantity(rv_e(k))) / si::kelvins; // subsidence rate w_LS = w_LS_fctr()(k * dz); diff --git a/src/opts_lgrngn.hpp b/src/opts_lgrngn.hpp index 9ae6f986ff..c39728ee60 100644 --- a/src/opts_lgrngn.hpp +++ b/src/opts_lgrngn.hpp @@ -163,6 +163,8 @@ void setopts_micro( // terminal velocity choice rt_params.cloudph_opts_init.terminal_velocity = libcloudphxx::lgrngn::vt_t::beard77fast; + rt_params.cloudph_opts_init.RH_formula = libcloudphxx::lgrngn::RH_formula_t::rv_cc; // use rv to be consistent with Lipps Hemler + // parsing --out_dry and --out_wet options values // the format is: "rmin:rmax|0,1,2;rmin:rmax|3;..." for (auto &opt : std::set({"out_dry", "out_wet"})) From 53dd5158e643ea7a2b8b2cc9c0be3a33971ba59d Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 27 Jul 2018 17:03:50 +0200 Subject: [PATCH 51/74] thl profile: follow that th is no th_std --- drawbicyc/plot_prof.hpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/drawbicyc/plot_prof.hpp b/drawbicyc/plot_prof.hpp index a709eaee79..453fb4e002 100644 --- a/drawbicyc/plot_prof.hpp +++ b/drawbicyc/plot_prof.hpp @@ -410,22 +410,23 @@ void plot_profiles(Plotter_t plotter, Plots plots) // liquid potential temp [K] { auto &ql(res_tmp2); - auto &T(res_tmp); ql = plotter.h5load_ra_timestep(at * n["outfreq"]); // aerosol ql += plotter.h5load_rc_timestep(at * n["outfreq"]); // cloud ql += plotter.h5load_rr_timestep(at * n["outfreq"]); // rain // ql is now q_l (liq water content) // auto tmp = plotter.h5load_timestep("th", at * n["outfreq"]); // typename Plotter_t::arr_t th_d(tmp); - typename Plotter_t::arr_t th_d(plotter.h5load_timestep("th", at * n["outfreq"])); + typename Plotter_t::arr_t th(plotter.h5load_timestep("th", at * n["outfreq"])); ///auto tmp = plotter.h5load_timestep("rv", at * n["outfreq"]); typename Plotter_t::arr_t rv(plotter.h5load_timestep("rv", at * n["outfreq"])); + + typename Plotter_t::arr_t T(plotter.h5load_timestep("libcloud_temperature", at * n["outfreq"])); // init pressure, from rv just to get correct size - typename Plotter_t::arr_t p(rv); - T = pow(th_d * pow(rhod * R_d / (p_1000), R_d / c_pd), c_pd / (c_pd - R_d)); +// typename Plotter_t::arr_t p(rv); + // T = pow(th_d * pow(rhod * R_d / (p_1000), R_d / c_pd), c_pd / (c_pd - R_d)); // TODO: env pressure should be used below! - p = rhod * R_d * (1 + 29./18. * rv) * T; // Rv/Rd = 29/18 - res += pow(p_1000 / p, R_d / c_pd) * (T - ql * L / c_p); + // p = rhod * R_d * (1 + 29./18. * rv) * T; // Rv/Rd = 29/18 + res += th * (T - ql * L / c_p); // res += ql; } gp << "set title 'liquid potential temp [K]'\n"; From 2d119129f48d8511a811d580e58ba7e3d89f1054 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 27 Jul 2018 17:27:08 +0200 Subject: [PATCH 52/74] fix thl profile with th_std --- drawbicyc/plot_prof.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/drawbicyc/plot_prof.hpp b/drawbicyc/plot_prof.hpp index 453fb4e002..2ee9221501 100644 --- a/drawbicyc/plot_prof.hpp +++ b/drawbicyc/plot_prof.hpp @@ -426,7 +426,7 @@ void plot_profiles(Plotter_t plotter, Plots plots) // T = pow(th_d * pow(rhod * R_d / (p_1000), R_d / c_pd), c_pd / (c_pd - R_d)); // TODO: env pressure should be used below! // p = rhod * R_d * (1 + 29./18. * rv) * T; // Rv/Rd = 29/18 - res += th * (T - ql * L / c_p); + res += th / T * (T - ql * L / c_p); // res += ql; } gp << "set title 'liquid potential temp [K]'\n"; From b635e0dd9d265985921cf9ef5c5b78bb4728b518 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 30 Jul 2018 14:11:51 +0200 Subject: [PATCH 53/74] dycoms: theta_ref(0) = theta_virt_env(0) --- src/cases/DYCOMS98.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cases/DYCOMS98.hpp b/src/cases/DYCOMS98.hpp index 462299662e..64be7471e9 100644 --- a/src/cases/DYCOMS98.hpp +++ b/src/cases/DYCOMS98.hpp @@ -238,7 +238,7 @@ std::cout << "lwp env: " << lwp_env << std::endl; st(notopbot) = (th_e(notopbot+1) - th_e(notopbot-1)) / th_e(notopbot); real_t st_avg = blitz::sum(st) / (nz-2) / (2.*dz); // reference theta - th_ref = th_e(0) * exp(st_avg * k * dz); + th_ref = th_e(0) * (1. + 0.608 * rv_e(0)) * exp(st_avg * k * dz); // th_ref = th_e(0) * pow(1 + rv_e(0) / a, f) // calc dry theta at z=0 // * exp(st_avg * k * dz); // virtual temp at surface From 90398c3ffc139ce2e80a5aafd7eafdbe8cd58a27 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 30 Jul 2018 16:27:07 +0200 Subject: [PATCH 54/74] use hall kernel instead of hall_davis_no_waals --- src/opts_lgrngn.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/opts_lgrngn.hpp b/src/opts_lgrngn.hpp index c39728ee60..63ff7e024f 100644 --- a/src/opts_lgrngn.hpp +++ b/src/opts_lgrngn.hpp @@ -153,7 +153,7 @@ void setopts_micro( // coalescence kernel choice if(!onishi) - rt_params.cloudph_opts_init.kernel = libcloudphxx::lgrngn::kernel_t::hall_davis_no_waals; + rt_params.cloudph_opts_init.kernel = libcloudphxx::lgrngn::kernel_t::hall;//_davis_no_waals; else { rt_params.cloudph_opts_init.kernel = libcloudphxx::lgrngn::kernel_t::onishi_hall_davis_no_waals; From db95f2cb00ac076cad7674cd5de3fde10a46f532 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 30 Jul 2018 16:27:39 +0200 Subject: [PATCH 55/74] use the khvorostyanow_nonsppherical terminal velocity --- src/opts_lgrngn.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/opts_lgrngn.hpp b/src/opts_lgrngn.hpp index 63ff7e024f..6b2f4ba3a2 100644 --- a/src/opts_lgrngn.hpp +++ b/src/opts_lgrngn.hpp @@ -161,7 +161,7 @@ void setopts_micro( rt_params.cloudph_opts_init.kernel_parameters.push_back(ReL); } // terminal velocity choice - rt_params.cloudph_opts_init.terminal_velocity = libcloudphxx::lgrngn::vt_t::beard77fast; + rt_params.cloudph_opts_init.terminal_velocity = libcloudphxx::lgrngn::vt_t::khvorostyanov_nonspherical; rt_params.cloudph_opts_init.RH_formula = libcloudphxx::lgrngn::RH_formula_t::rv_cc; // use rv to be consistent with Lipps Hemler From 5111a601f3c732bd6e096305d721c46475d30a52 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Mon, 30 Jul 2018 16:28:19 +0200 Subject: [PATCH 56/74] use pv_cc instead of rv_cc - a little inconsistency with the lipps-hemler equations --- src/opts_lgrngn.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/opts_lgrngn.hpp b/src/opts_lgrngn.hpp index 6b2f4ba3a2..0798a3a24d 100644 --- a/src/opts_lgrngn.hpp +++ b/src/opts_lgrngn.hpp @@ -163,7 +163,7 @@ void setopts_micro( // terminal velocity choice rt_params.cloudph_opts_init.terminal_velocity = libcloudphxx::lgrngn::vt_t::khvorostyanov_nonspherical; - rt_params.cloudph_opts_init.RH_formula = libcloudphxx::lgrngn::RH_formula_t::rv_cc; // use rv to be consistent with Lipps Hemler + rt_params.cloudph_opts_init.RH_formula = libcloudphxx::lgrngn::RH_formula_t::pv_cc; // rv_cc would be more consistent with Lipps Hemler // parsing --out_dry and --out_wet options values // the format is: "rmin:rmax|0,1,2;rmin:rmax|3;..." From ad098dc827a1b9da579af8748dc50b23d21b5096 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 31 Jul 2018 11:12:57 +0200 Subject: [PATCH 57/74] Revert "use hall kernel instead of hall_davis_no_waals" This reverts commit 90398c3ffc139ce2e80a5aafd7eafdbe8cd58a27. --- src/opts_lgrngn.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/opts_lgrngn.hpp b/src/opts_lgrngn.hpp index 0798a3a24d..6fac1052d7 100644 --- a/src/opts_lgrngn.hpp +++ b/src/opts_lgrngn.hpp @@ -153,7 +153,7 @@ void setopts_micro( // coalescence kernel choice if(!onishi) - rt_params.cloudph_opts_init.kernel = libcloudphxx::lgrngn::kernel_t::hall;//_davis_no_waals; + rt_params.cloudph_opts_init.kernel = libcloudphxx::lgrngn::kernel_t::hall_davis_no_waals; else { rt_params.cloudph_opts_init.kernel = libcloudphxx::lgrngn::kernel_t::onishi_hall_davis_no_waals; From 3f26c43a7a3a1f35474b4bed4016926d054dba01 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 31 Jul 2018 11:13:21 +0200 Subject: [PATCH 58/74] Revert "use pv_cc instead of rv_cc - a little inconsistency with the lipps-hemler equations" This reverts commit 5111a601f3c732bd6e096305d721c46475d30a52. --- src/opts_lgrngn.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/opts_lgrngn.hpp b/src/opts_lgrngn.hpp index 6fac1052d7..e57cfe7ee8 100644 --- a/src/opts_lgrngn.hpp +++ b/src/opts_lgrngn.hpp @@ -163,7 +163,7 @@ void setopts_micro( // terminal velocity choice rt_params.cloudph_opts_init.terminal_velocity = libcloudphxx::lgrngn::vt_t::khvorostyanov_nonspherical; - rt_params.cloudph_opts_init.RH_formula = libcloudphxx::lgrngn::RH_formula_t::pv_cc; // rv_cc would be more consistent with Lipps Hemler + rt_params.cloudph_opts_init.RH_formula = libcloudphxx::lgrngn::RH_formula_t::rv_cc; // use rv to be consistent with Lipps Hemler // parsing --out_dry and --out_wet options values // the format is: "rmin:rmax|0,1,2;rmin:rmax|3;..." From b21bb82cf111a663fa6f51a72e5adb8e5ce4b783 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Tue, 31 Jul 2018 16:16:18 +0200 Subject: [PATCH 59/74] condensation after apply rhs at n --- src/slvr_lgrngn.hpp | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 6e5661f9b1..098c0e3519 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -403,7 +403,7 @@ class slvr_lgrngn : public slvr_common } this->mem->barrier(); - parent_t::hook_ante_step(); // includes output + parent_t::hook_ante_step(); // includes RHS, which in turn launches sync_in and step_cond } @@ -476,6 +476,19 @@ class slvr_lgrngn : public slvr_common void hook_mixed_rhs_ante_step() { + this->update_rhs(this->rhs, this->dt, 0); // TODO: update_rhs called twice per step causes halo filling twice (done by parent_t::update_rhs), probably not needed - we just need to set rhs to zero + this->apply_rhs(this->dt); + + // rv might be negative due to large negative RHS from SD fluctuations + large-scale subsidence? + // turn all negative rv into rv = 0... CHEATING + negtozero(this->mem->advectee(ix::rv), "rv before step sync"); + nancheck(this->mem->advectee(ix::th), "th before step sync"); + nancheck(this->mem->advectee(ix::rv), "rv before step sync"); + negcheck(this->mem->advectee(ix::th), "th before step sync"); + negcheck(this->mem->advectee(ix::rv), "rv before step sync"); + + this->mem->barrier(); + // pass Eulerian fields to microphysics if (this->rank == 0) { @@ -513,18 +526,6 @@ class slvr_lgrngn : public slvr_common } this->mem->barrier(); - this->update_rhs(this->rhs, this->dt, 0); // TODO: update_rhs called twice per step causes halo filling twice (done by parent_t::update_rhs), probably not needed - we just need to set rhs to zero - this->apply_rhs(this->dt); - - // rv might be negative due to large negative RHS from SD fluctuations + large-scale subsidence? - // turn all negative rv into rv = 0... CHEATING - negtozero(this->mem->advectee(ix::rv), "rv before step sync"); - nancheck(this->mem->advectee(ix::th), "th before step sync"); - nancheck(this->mem->advectee(ix::rv), "rv before step sync"); - negcheck(this->mem->advectee(ix::th), "th before step sync"); - negcheck(this->mem->advectee(ix::rv), "rv before step sync"); - - this->mem->barrier(); // start sync/async run of step_cond // step_cond takes th and rv only for sync_out purposes - the values of th and rv before condensation come from sync_in, i.e. before apply_rhs if (this->rank == 0) From b5ff2ff4731e29b7448faefa9c1fe0f186848e6a Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 1 Aug 2018 13:04:48 +0200 Subject: [PATCH 60/74] add rv/th pre/post cond arrays to slvr_lgrngn --- src/slvr_common.hpp | 6 ++---- src/slvr_lgrngn.hpp | 19 ++++++++++++++++++- 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index 26bd12a658..22d855e2ac 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -36,7 +36,6 @@ class slvr_common : public slvr_dim // global arrays, shared among threads, TODO: in fact no need to share them? typename parent_t::arr_t &tmp1, - &tmp2, &r_l, &F, // forcings helper &alpha, // 'explicit' rhs part - does not depend on the value at n+1 @@ -365,11 +364,10 @@ class slvr_common : public slvr_dim params(p), spinup(p.spinup), tmp1(args.mem->tmp[__FILE__][0][0]), - tmp2(args.mem->tmp[__FILE__][0][5]), r_l(args.mem->tmp[__FILE__][0][2]), alpha(args.mem->tmp[__FILE__][0][3]), beta(args.mem->tmp[__FILE__][0][4]), - F(args.mem->tmp[__FILE__][0][1]) + F(args.mem->tmp[__FILE__][0][1]), { k_i.resize(this->shape(this->hrzntl_domain)); // TODO: resize to hrzntl_subdomain surf_flux_sens.resize(this->shape(this->hrzntl_domain)); // TODO: resize to hrzntl_subdomain @@ -380,6 +378,6 @@ class slvr_common : public slvr_dim static void alloc(typename parent_t::mem_t *mem, const int &n_iters) { parent_t::alloc(mem, n_iters); - parent_t::alloc_tmp_sclr(mem, __FILE__, 6); // tmp1, tmp2, r_l, alpha, beta, F + parent_t::alloc_tmp_sclr(mem, __FILE__, 5); } }; diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 098c0e3519..8130768d38 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -22,6 +22,13 @@ class slvr_lgrngn : public slvr_common // member fields std::unique_ptr> prtcls; + // helpers for calculating RHS from condensation, probably some of the could be avoided e.g. if step_cond returnd deltas and not changed fields + // or if change in theta was calculated from change in rv + typename parent_t::arr_t &rv_pre_cond, + &rv_post_cond, + &th_pre_cond, + &th_post_cond; + // helper methods void diag() { @@ -639,10 +646,20 @@ class slvr_lgrngn : public slvr_common const rt_params_t &p ) : parent_t(args, p), - params(p) + params(p), + rv_pre_cond(args.mem->tmp[__FILE__][0][0]), + rv_post_cond(args.mem->tmp[__FILE__][0][1]), + th_pre_cond(args.mem->tmp[__FILE__][0][2]), + th_post_cond(args.mem->tmp[__FILE__][0][3]) { // TODO: equip rank() in libmpdata with an assert() checking if not in serial block } + static void alloc(typename parent_t::mem_t *mem, const int &n_iters) + { + parent_t::alloc(mem, n_iters); + parent_t::alloc_tmp_sclr(mem, __FILE__, 4); + } + }; From 609dc0b166cea1d4155bdb3971a23681035145bb Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 1 Aug 2018 14:24:11 +0200 Subject: [PATCH 61/74] slvr_lgrngn: store cond rhs and apply --- src/slvr_common.hpp | 2 +- src/slvr_lgrngn.hpp | 60 +++++++++++++++++++++++---------------------- 2 files changed, 32 insertions(+), 30 deletions(-) diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index 22d855e2ac..f258fcdc1f 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -367,7 +367,7 @@ class slvr_common : public slvr_dim r_l(args.mem->tmp[__FILE__][0][2]), alpha(args.mem->tmp[__FILE__][0][3]), beta(args.mem->tmp[__FILE__][0][4]), - F(args.mem->tmp[__FILE__][0][1]), + F(args.mem->tmp[__FILE__][0][1]) { k_i.resize(this->shape(this->hrzntl_domain)); // TODO: resize to hrzntl_subdomain surf_flux_sens.resize(this->shape(this->hrzntl_domain)); // TODO: resize to hrzntl_subdomain diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 8130768d38..a8b0a2836a 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -432,18 +432,24 @@ class slvr_lgrngn : public slvr_common assert(ftr.valid()); parent_t::tbeg = parent_t::clock::now(); ftr.get(); - // microphysics could have led to rv < 0 ? - // TODO: repeated above, unify - negtozero(this->mem->advectee(ix::rv), "rv after step sync"); - - nancheck(this->mem->advectee(ix::th), "th after step sync"); - nancheck(this->mem->advectee(ix::rv), "rv after step sync"); - negcheck(this->mem->advectee(ix::th), "th after step sync"); - negcheck(this->mem->advectee(ix::rv), "rv after step sync"); parent_t::tend = parent_t::clock::now(); parent_t::tsync_wait += std::chrono::duration_cast( parent_t::tend - parent_t::tbeg ); } else assert(!ftr.valid()); #endif + } + + // add microphysics contribution to th and rv + this->state(ix::rv)(this->ijk) += rv_post_cond(this->ijk) - rv_pre_cond(this->ijk); + this->state(ix::th)(this->ijk) += th_post_cond(this->ijk) - th_pre_cond(this->ijk); + + if (this->rank == 0) + { + // microphysics could have led to rv < 0 ? + negtozero(this->mem->advectee(ix::rv), "rv after step sync"); + nancheck(this->mem->advectee(ix::th), "th after step sync"); + nancheck(this->mem->advectee(ix::rv), "rv after step sync"); + negcheck(this->mem->advectee(ix::th), "th after step sync"); + negcheck(this->mem->advectee(ix::rv), "rv after step sync"); // running asynchronous stuff { @@ -494,6 +500,9 @@ class slvr_lgrngn : public slvr_common negcheck(this->mem->advectee(ix::th), "th before step sync"); negcheck(this->mem->advectee(ix::rv), "rv before step sync"); + rv_pre_cond(this->ijk) = this->state(ix::rv)(this->ijk); + th_pre_cond(this->ijk) = this->state(ix::th)(this->ijk); + this->mem->barrier(); // pass Eulerian fields to microphysics @@ -522,14 +531,14 @@ class slvr_lgrngn : public slvr_common using libcloudphxx::lgrngn::CUDA; using libcloudphxx::lgrngn::multi_CUDA; - prtcls->sync_in( - make_arrinfo(this->mem->advectee(ix::th)), - make_arrinfo(this->mem->advectee(ix::rv)), - libcloudphxx::lgrngn::arrinfo_t(), - make_arrinfo(Cx), - this->n_dims == 2 ? libcloudphxx::lgrngn::arrinfo_t() : make_arrinfo(Cy), - make_arrinfo(Cz) - ); + prtcls->sync_in( + make_arrinfo(this->mem->advectee(ix::th)), + make_arrinfo(this->mem->advectee(ix::rv)), + libcloudphxx::lgrngn::arrinfo_t(), + make_arrinfo(Cx), + this->n_dims == 2 ? libcloudphxx::lgrngn::arrinfo_t() : make_arrinfo(Cy), + make_arrinfo(Cz) + ); } this->mem->barrier(); @@ -551,8 +560,8 @@ class slvr_lgrngn : public slvr_common &particles_t::step_cond, dynamic_cast*>(prtcls.get()), params.cloudph_opts, - make_arrinfo(this->mem->advectee(ix::th)), - make_arrinfo(this->mem->advectee(ix::rv)), + make_arrinfo(th_post_cond(this->domain).reindex(this->zero)), + make_arrinfo(rv_post_cond(this->domain).reindex(this->zero)), std::map >() ); else if(params.backend == multi_CUDA) @@ -561,8 +570,8 @@ class slvr_lgrngn : public slvr_common &particles_t::step_cond, dynamic_cast*>(prtcls.get()), params.cloudph_opts, - make_arrinfo(this->mem->advectee(ix::th)), - make_arrinfo(this->mem->advectee(ix::rv)), + make_arrinfo(th_post_cond(this->domain).reindex(this->zero)), + make_arrinfo(rv_post_cond(this->domain).reindex(this->zero)), std::map >() ); assert(ftr.valid()); @@ -571,16 +580,9 @@ class slvr_lgrngn : public slvr_common { prtcls->step_cond( params.cloudph_opts, - make_arrinfo(this->mem->advectee(ix::th)), - make_arrinfo(this->mem->advectee(ix::rv)) + make_arrinfo(th_post_cond(this->domain).reindex(this->zero)), + make_arrinfo(rv_post_cond(this->domain).reindex(this->zero)) ); - // microphysics could have led to rv < 0 ? - negtozero(this->mem->advectee(ix::rv), "rv after step sync"); - - nancheck(this->mem->advectee(ix::th), "th after step sync"); - nancheck(this->mem->advectee(ix::rv), "rv after step sync"); - negcheck(this->mem->advectee(ix::th), "th after step sync"); - negcheck(this->mem->advectee(ix::rv), "rv after step sync"); } parent_t::tend = parent_t::clock::now(); From 25bdabfd9fc6e380944a76c2f7b9cbbead5fc25a Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 2 Aug 2018 11:28:12 +0200 Subject: [PATCH 62/74] lgrngn: better nancheck messages + barriers update --- src/slvr_lgrngn.hpp | 41 ++++++++++++++++++----------------------- 1 file changed, 18 insertions(+), 23 deletions(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index a8b0a2836a..67fe32301d 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -437,20 +437,20 @@ class slvr_lgrngn : public slvr_common } else assert(!ftr.valid()); #endif } + this->mem->barrier(); // add microphysics contribution to th and rv this->state(ix::rv)(this->ijk) += rv_post_cond(this->ijk) - rv_pre_cond(this->ijk); this->state(ix::th)(this->ijk) += th_post_cond(this->ijk) - th_pre_cond(this->ijk); + // microphysics could have led to rv < 0 ? + negtozero(this->mem->advectee(ix::rv)(this->ijk), "rv after condensation"); + nancheck(this->mem->advectee(ix::th)(this->ijk), "th after condensation"); + nancheck(this->mem->advectee(ix::rv)(this->ijk), "rv after condensation"); + negcheck(this->mem->advectee(ix::th)(this->ijk), "th after condensation"); + negcheck(this->mem->advectee(ix::rv)(this->ijk), "rv after condensation"); if (this->rank == 0) { - // microphysics could have led to rv < 0 ? - negtozero(this->mem->advectee(ix::rv), "rv after step sync"); - nancheck(this->mem->advectee(ix::th), "th after step sync"); - nancheck(this->mem->advectee(ix::rv), "rv after step sync"); - negcheck(this->mem->advectee(ix::th), "th after step sync"); - negcheck(this->mem->advectee(ix::rv), "rv after step sync"); - // running asynchronous stuff { using libcloudphxx::lgrngn::particles_t; @@ -483,7 +483,6 @@ class slvr_lgrngn : public slvr_common parent_t::tasync += std::chrono::duration_cast( parent_t::tend - parent_t::tbeg ); } } - this->mem->barrier(); } void hook_mixed_rhs_ante_step() @@ -494,11 +493,11 @@ class slvr_lgrngn : public slvr_common // rv might be negative due to large negative RHS from SD fluctuations + large-scale subsidence? // turn all negative rv into rv = 0... CHEATING - negtozero(this->mem->advectee(ix::rv), "rv before step sync"); - nancheck(this->mem->advectee(ix::th), "th before step sync"); - nancheck(this->mem->advectee(ix::rv), "rv before step sync"); - negcheck(this->mem->advectee(ix::th), "th before step sync"); - negcheck(this->mem->advectee(ix::rv), "rv before step sync"); + negtozero(this->mem->advectee(ix::rv)(this->ijk), "rv after mixed_rhs_ante_step apply rhs"); + nancheck(this->mem->advectee(ix::th)(this->ijk), "th after mixed_rhs_ante_step apply rhs"); + nancheck(this->mem->advectee(ix::rv)(this->ijk), "rv after mixed_rhs_ante_step apply rhs"); + negcheck(this->mem->advectee(ix::th)(this->ijk), "th after mixed_rhs_ante_step apply rhs"); + negcheck(this->mem->advectee(ix::rv)(this->ijk), "rv after mixed_rhs_ante_step apply rhs"); rv_pre_cond(this->ijk) = this->state(ix::rv)(this->ijk); th_pre_cond(this->ijk) = this->state(ix::th)(this->ijk); @@ -539,13 +538,9 @@ class slvr_lgrngn : public slvr_common this->n_dims == 2 ? libcloudphxx::lgrngn::arrinfo_t() : make_arrinfo(Cy), make_arrinfo(Cz) ); - } - this->mem->barrier(); - // start sync/async run of step_cond - // step_cond takes th and rv only for sync_out purposes - the values of th and rv before condensation come from sync_in, i.e. before apply_rhs - if (this->rank == 0) - { + // start sync/async run of step_cond + // step_cond takes th and rv only for sync_out purposes - the values of th and rv before condensation come from sync_in, i.e. before apply_rhs using libcloudphxx::lgrngn::particles_t; using libcloudphxx::lgrngn::CUDA; using libcloudphxx::lgrngn::multi_CUDA; @@ -599,10 +594,10 @@ class slvr_lgrngn : public slvr_common this->apply_rhs(this->dt); // no negtozero after apply, because only w changed here // TODO: add these nanchecks/negchecks to apply_rhs, since they are repeated twice now - nancheck(this->mem->advectee(ix::th), "th after step sync"); - nancheck(this->mem->advectee(ix::rv), "rv after step sync"); - negcheck(this->mem->advectee(ix::th), "th after step sync"); - negcheck(this->mem->advectee(ix::rv), "rv after step sync"); + nancheck(this->mem->advectee(ix::th)(this->ijk), "th after mixed_rhs_post_step apply rhs"); + nancheck(this->mem->advectee(ix::rv)(this->ijk), "rv after mixed_rhs_post_step apply rhs"); + negcheck(this->mem->advectee(ix::th)(this->ijk), "th after mixed_rhs_post_step apply rhs"); + negcheck(this->mem->advectee(ix::rv)(this->ijk), "rv after mixed_rhs_post_step apply rhs"); } void record_all() From f78ddb7d9d6fc31029c9f5791c60b89be8ccc2c6 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 2 Aug 2018 11:30:44 +0200 Subject: [PATCH 63/74] lgrngn: rhs at n done after running step_sync --- src/slvr_lgrngn.hpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 67fe32301d..f5f1aa2d7f 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -488,17 +488,6 @@ class slvr_lgrngn : public slvr_common void hook_mixed_rhs_ante_step() { - this->update_rhs(this->rhs, this->dt, 0); // TODO: update_rhs called twice per step causes halo filling twice (done by parent_t::update_rhs), probably not needed - we just need to set rhs to zero - this->apply_rhs(this->dt); - - // rv might be negative due to large negative RHS from SD fluctuations + large-scale subsidence? - // turn all negative rv into rv = 0... CHEATING - negtozero(this->mem->advectee(ix::rv)(this->ijk), "rv after mixed_rhs_ante_step apply rhs"); - nancheck(this->mem->advectee(ix::th)(this->ijk), "th after mixed_rhs_ante_step apply rhs"); - nancheck(this->mem->advectee(ix::rv)(this->ijk), "rv after mixed_rhs_ante_step apply rhs"); - negcheck(this->mem->advectee(ix::th)(this->ijk), "th after mixed_rhs_ante_step apply rhs"); - negcheck(this->mem->advectee(ix::rv)(this->ijk), "rv after mixed_rhs_ante_step apply rhs"); - rv_pre_cond(this->ijk) = this->state(ix::rv)(this->ijk); th_pre_cond(this->ijk) = this->state(ix::th)(this->ijk); @@ -586,6 +575,17 @@ class slvr_lgrngn : public slvr_common parent_t::tbeg = parent_t::clock::now(); } this->mem->barrier(); + + this->update_rhs(this->rhs, this->dt, 0); // TODO: update_rhs called twice per step causes halo filling twice (done by parent_t::update_rhs), probably not needed - we just need to set rhs to zero + this->apply_rhs(this->dt); + + // rv might be negative due to large negative RHS from SD fluctuations + large-scale subsidence? + // turn all negative rv into rv = 0... CHEATING + negtozero(this->mem->advectee(ix::rv)(this->ijk), "rv after mixed_rhs_ante_step apply rhs"); + nancheck(this->mem->advectee(ix::th)(this->ijk), "th after mixed_rhs_ante_step apply rhs"); + nancheck(this->mem->advectee(ix::rv)(this->ijk), "rv after mixed_rhs_ante_step apply rhs"); + negcheck(this->mem->advectee(ix::th)(this->ijk), "th after mixed_rhs_ante_step apply rhs"); + negcheck(this->mem->advectee(ix::rv)(this->ijk), "rv after mixed_rhs_ante_step apply rhs"); } void hook_mixed_rhs_post_step() From 55744d3125cc2cc9d8e719de804c165f86eba960 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 2 Aug 2018 12:37:55 +0200 Subject: [PATCH 64/74] lgrngn: diag_rl post cond, advect and subside it --- src/slvr_lgrngn.hpp | 51 ++++++++++++++++++++++++++++++++++++--------- 1 file changed, 41 insertions(+), 10 deletions(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index f5f1aa2d7f..4b3a671fb3 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -29,6 +29,18 @@ class slvr_lgrngn : public slvr_common &th_pre_cond, &th_post_cond; + void diag_rl() + { + if(this->rank != 0) return; + + prtcls->diag_all(); + prtcls->diag_wet_mom(3); + auto rl = parent_t::r_l(this->domain); // rl refrences subdomain of r_l + rl = typename parent_t::arr_t(prtcls->outbuf(), rl.shape(), blitz::duplicateData); // copy in data from outbuf; total liquid third moment of wet radius per kg of dry air [m^3 / kg] + nancheck(rl, "rl after copying from diag_wet_mom(3)"); + rl = rl * 4./3. * 1000. * 3.14159; // get mixing ratio [kg/kg] + } + // helper methods void diag() { @@ -375,7 +387,9 @@ class slvr_lgrngn : public slvr_common } void hook_mixed_rhs_ante_loop() - {} // empty, because update_rhs called before each step; defined, to avoid assert from default hook_mixed_rhs_ante_loop + { + diag_rl(); // init r_l + } #if defined(STD_FUTURE_WORKS) std::future ftr; @@ -399,15 +413,6 @@ class slvr_lgrngn : public slvr_common parent_t::tasync_wait += std::chrono::duration_cast( parent_t::tend - parent_t::tbeg ); } else assert(!ftr.valid()); #endif - - // store liquid water content to be used in update_rhs (if done in update_rhs, it fails on async runs) - prtcls->diag_all(); - prtcls->diag_wet_mom(3); - auto rl = parent_t::r_l(this->domain); // rl refrences subdomain of r_l - rl = typename parent_t::arr_t(prtcls->outbuf(), rl.shape(), blitz::duplicateData); // copy in data from outbuf; total liquid third moment of wet radius per kg of dry air [m^3 / kg] - nancheck(rl, "rl after copying from diag_wet_mom(3)"); - rl = rl * 4./3. * 1000. * 3.14159; // get mixing ratio [kg/kg] - } this->mem->barrier(); parent_t::hook_ante_step(); // includes RHS, which in turn launches sync_in and step_cond @@ -448,6 +453,9 @@ class slvr_lgrngn : public slvr_common nancheck(this->mem->advectee(ix::rv)(this->ijk), "rv after condensation"); negcheck(this->mem->advectee(ix::th)(this->ijk), "th after condensation"); negcheck(this->mem->advectee(ix::rv)(this->ijk), "rv after condensation"); + + // store liquid water content (post-cond, pre-adve and pre-subsidence) + diag_rl(); if (this->rank == 0) { @@ -483,6 +491,29 @@ class slvr_lgrngn : public slvr_common parent_t::tasync += std::chrono::duration_cast( parent_t::tend - parent_t::tbeg ); } } + this->mem->barrier(); + + // subsidence of rl + // TODO: very similar code to subsidence function in forcings.hppp + if(params.subsidence) + { + auto& tmp1(this->tmp1); + auto& r_l(this->r_l); + const auto& ijk(this->ijk); + const auto& params(this->params); + auto& F(this->F); + + tmp1(ijk) = r_l(ijk); + // fill halos for gradient calculation + // TODO: no need to xchng in horizontal, which potentially causes MPI communication + this->xchng_sclr(tmp1, this->ijk, this->halo); + this->vert_grad_cnt(tmp1, F, params.dz); + F(ijk).reindex(this->zero) *= - (*params.w_LS)(this->vert_idx); + r_l(ijk) += F(ijk) * this->dt; + } + + // advect r_l (1st-order) + this->self_advec_donorcell(this->r_l); } void hook_mixed_rhs_ante_step() From 7d61fc9b01b3fb7dc4a279ee4b91106652ce4853 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 29 Aug 2018 11:52:07 +0200 Subject: [PATCH 65/74] fix dry and moist thermals with pressure --- src/cases/DryThermalGMD2015.hpp | 22 +++++++++++++++++----- src/cases/MoistThermalGrabowskiClark99.hpp | 6 +++++- 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/src/cases/DryThermalGMD2015.hpp b/src/cases/DryThermalGMD2015.hpp index 83672b1426..a8b280adbb 100644 --- a/src/cases/DryThermalGMD2015.hpp +++ b/src/cases/DryThermalGMD2015.hpp @@ -87,11 +87,23 @@ namespace setup rhod = 1; th_e = 300; - //p_e = pow((rhod * si::kilograms / si::cubic_metres) * R_d() * (th_e * si::kelvins) * pow(p_1000(), R_d_over_c_pd()), 1. / (R_d_over_c_pd() + 1)) / si::pascals; -// const quantity p_theta = (quantity(rhod * si::kilograms / si::cubic_metres) * R_d() * quantity(th_e * si::kelvins))/si::pascals; - const quantity p_theta = (quantity(1. * si::kilograms / si::cubic_metres) * R_d() * quantity(300. * si::kelvins)); - /*const quantity*/ real_t p = pow( real_t(p_theta / si::pascals) * pow(real_t(p_1000() / si::pascals), R_d_over_c_pd()), 1. / (R_d_over_c_pd() + 1)) ;// * pow(p_1000(), R_d_over_c_pd()));//, 1. / (R_d_over_c_pd() + 1)) / si::pascals; - p_e = p; // total env pressure + + const quantity T( + libcloudphxx::common::theta_dry::T( + quantity(300 * si::kelvins), + quantity(1 * si::kilograms / si::cubic_metres) + ) + ); + + const quantity p( + libcloudphxx::common::theta_dry::p( + quantity(1 * si::kilograms / si::cubic_metres), + quantity(0.), + T + ) + ); + p_e = real_t(p / si::pascals); // total env pressure + th_ref = 300; } }; diff --git a/src/cases/MoistThermalGrabowskiClark99.hpp b/src/cases/MoistThermalGrabowskiClark99.hpp index 592fea8adf..dc310b284a 100644 --- a/src/cases/MoistThermalGrabowskiClark99.hpp +++ b/src/cases/MoistThermalGrabowskiClark99.hpp @@ -355,6 +355,10 @@ namespace setup { blitz::thirdIndex k; this->intcond_hlpr(solver, rhod, th_e, rv_e, rng_seed, k); + + arr_1D_t p_d_e(p_e - detail::calc_p_v()(p_e, rv_e)); + arr_1D_t T(th_e * pow(p_d_e / 1.e5, R_d_over_c_pd())); + using ix = typename concurr_t::solver_t::ix; int nz = solver.advectee().extent(2); real_t dz = (Z / si::metres) / (nz-1); @@ -362,7 +366,7 @@ namespace setup real_t dx = (X / si::metres) / (nx-1); int ny = solver.advectee().extent(1); real_t dy = (Y / si::metres) / (ny-1); - solver.advectee(ix::rv) = prtrb_rv(th_e, rhod, dz)( + solver.advectee(ix::rv) = prtrb_rv(T, p_e, dz)( sqrt( pow(blitz::tensor::i * dx - (X / si::metres / 2.), 2) + pow(blitz::tensor::j * dy - (Y / si::metres / 2.), 2) + From fd2c56b885efb64e01a42b9e05f6ef2ed0511404 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 29 Aug 2018 11:57:24 +0200 Subject: [PATCH 66/74] conform to codacy's stupid requests --- src/detail/checknan.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/detail/checknan.cpp b/src/detail/checknan.cpp index 24840a86ab..7c189de0cc 100644 --- a/src/detail/checknan.cpp +++ b/src/detail/checknan.cpp @@ -29,7 +29,7 @@ namespace nancheck_hlprs { template - void nancheck_hlpr(const arr_t &arr, std::string name) + void nancheck_hlpr(const arr_t &arr, const std::string &name) { if(!std::isfinite(sum(arr))) { @@ -43,7 +43,7 @@ namespace nancheck_hlprs } template - void nancheck2_hlpr(const arr_t &arrcheck, const arr_t &arrout, std::string name) + void nancheck2_hlpr(const arr_t &arrcheck, const arr_t &arrout, const std::string &name) { if(!std::isfinite(sum(arrcheck))) { @@ -58,7 +58,7 @@ namespace nancheck_hlprs } template - void negcheck_hlpr(const arr_t &arr, std::string name) + void negcheck_hlpr(const arr_t &arr, const std::string &name) { if(min(arr) < 0.) { @@ -72,7 +72,7 @@ namespace nancheck_hlprs } template - void negtozero_hlpr(arr_t arr, std::string name) + void negtozero_hlpr(arr_t arr, const std::string &name) { if(min(arr) < 0.) { From af94e8a9f9750b0b6b14b52e8f39d13b53fec426 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 29 Aug 2018 14:12:59 +0200 Subject: [PATCH 67/74] use rv_cc RH formula --- src/cases/MoistThermalGrabowskiClark99.hpp | 10 ++++++++-- src/opts_lgrngn.hpp | 2 ++ 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/cases/MoistThermalGrabowskiClark99.hpp b/src/cases/MoistThermalGrabowskiClark99.hpp index dc310b284a..f49d3597b4 100644 --- a/src/cases/MoistThermalGrabowskiClark99.hpp +++ b/src/cases/MoistThermalGrabowskiClark99.hpp @@ -32,13 +32,19 @@ namespace setup using libcloudphxx::common::theta_std::p_1000; - // RH T and p to rv + // RH T and p to rv assuming RH = r_v / r_vs quantity RH_T_p_to_rv(const real_t &RH, const quantity &T, const quantity &p) { - return moist_air::eps() * RH * const_cp::p_vs(T) / (p - RH * const_cp::p_vs(T)); + return RH * const_cp::r_vs(T, p); } /* + // RH T and p to rv assuming RH = p_v / p_vs + quantity RH_T_p_to_rv(const real_t &RH, const quantity &T, const quantity &p) + { + return moist_air::eps() * RH * const_cp::p_vs(T) / (p - RH * const_cp::p_vs(T)); + } + // rv(RH, th_dry, rhod) real_t RH_th_rhod_to_rv(const real_t &RH, const real_t &th_dry, const real_t &rhod) { diff --git a/src/opts_lgrngn.hpp b/src/opts_lgrngn.hpp index 9ae6f986ff..c39728ee60 100644 --- a/src/opts_lgrngn.hpp +++ b/src/opts_lgrngn.hpp @@ -163,6 +163,8 @@ void setopts_micro( // terminal velocity choice rt_params.cloudph_opts_init.terminal_velocity = libcloudphxx::lgrngn::vt_t::beard77fast; + rt_params.cloudph_opts_init.RH_formula = libcloudphxx::lgrngn::RH_formula_t::rv_cc; // use rv to be consistent with Lipps Hemler + // parsing --out_dry and --out_wet options values // the format is: "rmin:rmax|0,1,2;rmin:rmax|3;..." for (auto &opt : std::set({"out_dry", "out_wet"})) From 08def4912105ddbc2bfead826056541276b356c6 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Wed, 29 Aug 2018 14:35:08 +0200 Subject: [PATCH 68/74] moist thermal: following th_std --- src/cases/MoistThermalGrabowskiClark99.hpp | 10 ++++++---- src/opts_lgrngn.hpp | 2 -- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/cases/MoistThermalGrabowskiClark99.hpp b/src/cases/MoistThermalGrabowskiClark99.hpp index f49d3597b4..1fbab95733 100644 --- a/src/cases/MoistThermalGrabowskiClark99.hpp +++ b/src/cases/MoistThermalGrabowskiClark99.hpp @@ -287,12 +287,14 @@ namespace setup th_e = th_e / (1. + a * rv_e); // turn standard potential temp into dry potential temp +/* for(int k=0; k si_rv_e = rv_e(k); th_e(k) = libcloudphxx::common::theta_dry::std2dry(th_e(k) * si::kelvins, si_rv_e) / si::kelvins; real_t p_d = pre_ref(k) - libcloudphxx::common::moist_air::p_v(pre_ref(k) * si::pascals, rv_e(k)) / si::pascals; } +*/ th_ref = th_e;//th_std_fctr(th_std_0 / si::kelvins)(k * dz); } @@ -323,8 +325,8 @@ namespace setup blitz::secondIndex k; this->intcond_hlpr(solver, rhod, th_e, rv_e, rng_seed, k); - arr_1D_t p_d_e(p_e - detail::calc_p_v()(p_e, rv_e)); - arr_1D_t T(th_e * pow(p_d_e / 1.e5, R_d_over_c_pd())); +// arr_1D_t p_d_e(p_e - detail::calc_p_v()(p_e, rv_e)); + arr_1D_t T(th_e * pow(p_e / 1.e5, R_d_over_c_pd())); using ix = typename concurr_t::solver_t::ix; int nz = solver.advectee().extent(ix::w); @@ -362,8 +364,8 @@ namespace setup blitz::thirdIndex k; this->intcond_hlpr(solver, rhod, th_e, rv_e, rng_seed, k); - arr_1D_t p_d_e(p_e - detail::calc_p_v()(p_e, rv_e)); - arr_1D_t T(th_e * pow(p_d_e / 1.e5, R_d_over_c_pd())); +// arr_1D_t p_d_e(p_e - detail::calc_p_v()(p_e, rv_e)); + arr_1D_t T(th_e * pow(p_e / 1.e5, R_d_over_c_pd())); using ix = typename concurr_t::solver_t::ix; int nz = solver.advectee().extent(2); diff --git a/src/opts_lgrngn.hpp b/src/opts_lgrngn.hpp index 9818f936f0..e57cfe7ee8 100644 --- a/src/opts_lgrngn.hpp +++ b/src/opts_lgrngn.hpp @@ -165,8 +165,6 @@ void setopts_micro( rt_params.cloudph_opts_init.RH_formula = libcloudphxx::lgrngn::RH_formula_t::rv_cc; // use rv to be consistent with Lipps Hemler - rt_params.cloudph_opts_init.RH_formula = libcloudphxx::lgrngn::RH_formula_t::rv_cc; // use rv to be consistent with Lipps Hemler - // parsing --out_dry and --out_wet options values // the format is: "rmin:rmax|0,1,2;rmin:rmax|3;..." for (auto &opt : std::set({"out_dry", "out_wet"})) From 261273438f6a8812cccba45144cd173464d5ae3c Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 30 Aug 2018 12:00:32 +0200 Subject: [PATCH 69/74] set a valid Cy_domain in 2D to silence asserts --- src/slvr_dim.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/slvr_dim.hpp b/src/slvr_dim.hpp index b5300db218..9cd36be071 100644 --- a/src/slvr_dim.hpp +++ b/src/slvr_dim.hpp @@ -42,7 +42,7 @@ class slvr_dim< rng_t hrzntl_domain = this->mem->grid_size[0]; rng_t hrzntl_subdomain = this->i; idx_t<2> Cx_domain = idx_t<2>({this->mem->grid_size[0]^h, this->mem->grid_size[1]}); - idx_t<2> Cy_domain; + idx_t<2> Cy_domain = idx_t<2>({this->mem->grid_size[0], this->mem->grid_size[1]^h}); // just fill in with Cz_domain to avoid some asserts idx_t<2> Cz_domain = idx_t<2>({this->mem->grid_size[0], this->mem->grid_size[1]^h}); blitz::TinyVector zero = blitz::TinyVector({0,0}); From f102057beba8b6b3b083493ebe50be60b578aa9f Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 30 Aug 2018 12:47:38 +0200 Subject: [PATCH 70/74] MT test: tweak lgrngn results a bit to constant pressure changes --- tests/physics/moist_thermal.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/physics/moist_thermal.cpp b/tests/physics/moist_thermal.cpp index 4da4c40e71..e72e19bc43 100644 --- a/tests/physics/moist_thermal.cpp +++ b/tests/physics/moist_thermal.cpp @@ -75,8 +75,8 @@ int main(int ac, char** av) "rc_avg", // expected values map { - {"blk_1m", {{ 0, 0, 0.185215, 0.355712, 0.520675, 0.63674, 0.720535, 0.73, 0.68, 0.55, 0.44 }}}, - {"lgrngn", {{ 0, 0, 0.179877, 0.374551, 0.552171, 0.704754, 0.790675, 0.83, 0.83, 0.8, 0.73}}} + {"blk_1m", {{ 0, 0 , 0.185215, 0.355712, 0.520675, 0.63674, 0.720535, 0.73, 0.68, 0.55, 0.44 }}}, + {"lgrngn", {{ 0, 0.111, 0.19 , 0.374551, 0.552171, 0.704754, 0.790675, 0.83, 0.83, 0.8, 0.73}}} }, // epsilons map { @@ -142,7 +142,7 @@ int main(int ac, char** av) // expected values map { {"blk_1m", {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}}, - {"lgrngn", {{0, 0.476906, 0.428138, 0.374729, 0.333207, 0.277702, 0.224544, 0.213603, 0.251973, 0.298691, 0.351449}}} + {"lgrngn", {{0, 0.513454, 0.452787, 0.409429, 0.37853, 0.286295, 0.250789, 0.241939, 0.226568, 0.331526, 0.2}}} }, // epsilons map From a2d63187cabb9d8046a76651b8371fa921b975d9 Mon Sep 17 00:00:00 2001 From: pdziekan Date: Thu, 30 Aug 2018 15:56:42 +0200 Subject: [PATCH 71/74] MT test: blk_1m th_std exp results --- tests/physics/moist_thermal.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/physics/moist_thermal.cpp b/tests/physics/moist_thermal.cpp index e72e19bc43..226610892c 100644 --- a/tests/physics/moist_thermal.cpp +++ b/tests/physics/moist_thermal.cpp @@ -75,13 +75,13 @@ int main(int ac, char** av) "rc_avg", // expected values map { - {"blk_1m", {{ 0, 0 , 0.185215, 0.355712, 0.520675, 0.63674, 0.720535, 0.73, 0.68, 0.55, 0.44 }}}, - {"lgrngn", {{ 0, 0.111, 0.19 , 0.374551, 0.552171, 0.704754, 0.790675, 0.83, 0.83, 0.8, 0.73}}} + {"blk_1m", {{ 0, 0 , 0.19, 0.37 , 0.533 , 0.675 , 0.74 , 0.74, 0.66, 0.5, 0.39 }}}, + {"lgrngn", {{ 0, 0.111, 0.19, 0.374551, 0.552171, 0.704754, 0.790675, 0.83, 0.83, 0.8, 0.73}}} }, // epsilons map { {"blk_1m", {{1e-2, 2e-2, 2e-2, 2e-2, 2e-2, 2e-2, 2e-2, 2e-2, 3e-2, 5e-2, 2e-1}} }, - {"lgrngn", {{7e-2, 7e-2, 7e-2, 7e-2, 7e-2, 7e-2, 7e-2, 7e-2, 7e-2, 15e-2, 20e-2}} } + {"lgrngn", {{7e-2, 10e-2, 7e-2, 7e-2, 7e-2, 7e-2, 7e-2, 7e-2, 7e-2, 15e-2, 20e-2}} } } }); @@ -91,8 +91,8 @@ int main(int ac, char** av) "rc_std_dev", // expected values map { - {"blk_1m", {{0, 0, 0.0187748, 0.0711643, 0.136795, 0.211169, 0.264176, 0.311335, 0.329516, 0.290624, 0.237213}}}, - {"lgrngn", {{0, 0, 0.0551891, 0.12, 0.202856, 0.290047, 0.375072, 0.442025, 0.485523, 0.525853, 0.523559}}} + {"blk_1m", {{0, 0, 0.022 , 0.076, 0.154 , 0.206, 0.271, 0.332, 0.326, 0.262, 0.228}}}, + {"lgrngn", {{0, 0, 0.0551891, 0.12 , 0.202856, 0.290047, 0.375072, 0.442025, 0.485523, 0.525853, 0.523559}}} }, // epsilons map { @@ -142,7 +142,7 @@ int main(int ac, char** av) // expected values map { {"blk_1m", {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}}, - {"lgrngn", {{0, 0.513454, 0.452787, 0.409429, 0.37853, 0.286295, 0.250789, 0.241939, 0.226568, 0.331526, 0.2}}} + {"lgrngn", {{0, 0.513454, 0.452787, 0.409429, 0.37853, 0.286295, 0.250789, 0.241939, 0.226568, 0.331526, 0.}}} // 0 at the end because there are huge differences between runs }, // epsilons map @@ -241,7 +241,7 @@ int main(int ac, char** av) "clfrac", // expected values map { - {"blk_1m", {{0, 0, 0.0232409, 0.0249304, 0.0257979, 0.0268024, 0.0265285, 0.0248847, 0.0221908, 0.0195881, 0.0154787}}}, + {"blk_1m", {{0, 0, 0.0232409, 0.0246, 0.0257979, 0.0256, 0.0255, 0.0239, 0.0214, 0.0191, 0.0123}}}, {"lgrngn", {{0, 0, 0.0210493, 0.0241998, 0.024976, 0.0257522, 0.0259349, 0.0247477, 0.0239258, 0.0208666, 0.0168029}}} }, // epsilons map From 74286e021f82cbe564154bf8d2ef267f991e1c3e Mon Sep 17 00:00:00 2001 From: pdziekan Date: Tue, 4 Sep 2018 12:08:11 +0200 Subject: [PATCH 72/74] if usong large tail, set higher n_sd_max --- src/slvr_lgrngn.hpp | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index 4b3a671fb3..cc07a7bbc8 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -271,7 +271,12 @@ class slvr_lgrngn : public slvr_common params.cloudph_opts_init.z1 = (params.cloudph_opts_init.nz - .5) * this->dj; if(params.cloudph_opts_init.sd_conc) - params.cloudph_opts_init.n_sd_max = params.cloudph_opts_init.nx * params.cloudph_opts_init.nz * params.cloudph_opts_init.sd_conc; + { + if(params.cloudph_opts_init.sd_conc_large_tail) + params.cloudph_opts_init.n_sd_max = 1.2 * params.cloudph_opts_init.nx * params.cloudph_opts_init.nz * params.cloudph_opts_init.sd_conc; /// 1.2 to make space for large tail + else + params.cloudph_opts_init.n_sd_max = params.cloudph_opts_init.nx * params.cloudph_opts_init.nz * params.cloudph_opts_init.sd_conc; + } else params.cloudph_opts_init.n_sd_max = 1.1 * params.cloudph_opts_init.nx * params.cloudph_opts_init.nz * 1.e8 * params.cloudph_opts_init.dx * params.cloudph_opts_init.dz / params.cloudph_opts_init.sd_const_multi; // hardcoded N_a=100/cm^3 !! @@ -291,7 +296,12 @@ class slvr_lgrngn : public slvr_common params.cloudph_opts_init.z1 = (params.cloudph_opts_init.nz - .5) * this->dk; if(params.cloudph_opts_init.sd_conc) - params.cloudph_opts_init.n_sd_max = params.cloudph_opts_init.nx * params.cloudph_opts_init.ny *params.cloudph_opts_init.nz * params.cloudph_opts_init.sd_conc; + { + if(params.cloudph_opts_init.sd_conc_large_tail) + params.cloudph_opts_init.n_sd_max = 1.2 * params.cloudph_opts_init.nx * params.cloudph_opts_init.ny * params.cloudph_opts_init.nz * params.cloudph_opts_init.sd_conc; /// 1.2 to make space for large tail + else + params.cloudph_opts_init.n_sd_max = params.cloudph_opts_init.nx * params.cloudph_opts_init.ny * params.cloudph_opts_init.nz * params.cloudph_opts_init.sd_conc; + } else params.cloudph_opts_init.n_sd_max = 1.1 * params.cloudph_opts_init.nx * params.cloudph_opts_init.ny * params.cloudph_opts_init.nz * 1.e8 * params.cloudph_opts_init.dx * params.cloudph_opts_init.dy * params.cloudph_opts_init.dz / params.cloudph_opts_init.sd_const_multi; // hardcoded N_a=100/cm^3 !! From 87124ae8d5bd7fe098055f631f69a2a07e76216f Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Fri, 7 Sep 2018 15:23:40 +0200 Subject: [PATCH 73/74] fluxes: rhod and p_e from ground level --- src/calc_forces_common.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/calc_forces_common.hpp b/src/calc_forces_common.hpp index cdf8a44920..105bb72f95 100644 --- a/src/calc_forces_common.hpp +++ b/src/calc_forces_common.hpp @@ -39,7 +39,7 @@ void slvr_common::rv_src() // change of rv[1/s] = latent heating[W/m^3] / lat_heat_of_evap[J/kg] / density[kg/m^3] if(params.ForceParameters.surf_latent_flux_in_watts_per_square_meter) - alpha(ijk).reindex(this->zero) /= (libcloudphxx::common::const_cp::l_tri() * si::kilograms / si::joules) * (*params.rhod)(this->vert_idx); + alpha(ijk).reindex(this->zero) /= (libcloudphxx::common::const_cp::l_tri() * si::kilograms / si::joules) * (*params.rhod)(0);//this->vert_idx); // large-scale vertical wind subsidence(ix::rv); @@ -80,9 +80,9 @@ void slvr_common::th_src(typename parent_t::arr_t &rv) } // change of theta[K/s] = heating[W/m^3] / exner / c_p[J/K/kg] / this->rhod[kg/m^3] - alpha(ijk).reindex(this->zero) /= calc_exner()((*params.p_e)(this->vert_idx)) * + alpha(ijk).reindex(this->zero) /= calc_exner()((*params.p_e)(0)) * calc_c_p()(rv(ijk).reindex(this->zero)) * - (*params.rhod)(this->vert_idx); + (*params.rhod)(0); nancheck2(alpha(ijk), this->state(ix::th)(ijk), "change of theta"); From dc2be3472fe8a11146693fb299ad3bb3186458b2 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 10 Sep 2018 14:58:19 +0200 Subject: [PATCH 74/74] siurf flux: calc at ground --- src/calc_forces_common.hpp | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/calc_forces_common.hpp b/src/calc_forces_common.hpp index 105bb72f95..0f215faf3d 100644 --- a/src/calc_forces_common.hpp +++ b/src/calc_forces_common.hpp @@ -65,26 +65,36 @@ void slvr_common::th_src(typename parent_t::arr_t &rv) // radiation radiation(rv); - nancheck(beta(ijk), "radiation"); + nancheck(F(ijk), "radiation"); // sum of th flux, F(j) is upward flux through the bottom of the j-th cell this->vert_grad_fwd(F, alpha, params.dz); alpha(ijk) *= -1; // negative gradient means inflow nancheck(alpha(ijk), "sum of th flux"); + + // change of theta[K/s] = heating[W/m^3] / exner / c_p[J/K/kg] / this->rhod[kg/m^3] + alpha(ijk).reindex(this->zero) /= calc_exner()((*params.p_e)(this->vert_idx)) * + calc_c_p()(rv(ijk).reindex(this->zero)) * + (*params.rhod)(this->vert_idx); + + nancheck2(alpha(ijk), this->state(ix::th)(ijk), "change of theta"); // surface flux if(params.ForceParameters.surf_sensible_flux_in_watts_per_square_meter) { surf_sens(); nancheck(F(ijk), "sensible surf forcing"); - alpha(ijk) += F(ijk); - } + + constexpr int perm_no=ix::w; // 1 for 2D, 2 for 3D + auto ground = idxperm::pi(0, this->hrzntl_subdomain); - // change of theta[K/s] = heating[W/m^3] / exner / c_p[J/K/kg] / this->rhod[kg/m^3] - alpha(ijk).reindex(this->zero) /= calc_exner()((*params.p_e)(0)) * - calc_c_p()(rv(ijk).reindex(this->zero)) * - (*params.rhod)(0); + // change of theta[K/s] = heating[W/m^3] / exner / c_p[J/K/kg] / this->rhod[kg/m^3] + F(ijk).reindex(this->zero) /= calc_exner()((*params.p_e)(0)) * + calc_c_p()(rv(ijk).reindex(this->zero)) * // TODO: should be rv(ground) here! + (*params.rhod)(0); - nancheck2(alpha(ijk), this->state(ix::th)(ijk), "change of theta"); + nancheck(F(ijk), "sens surf force theta change"); + alpha(ijk) += F(ijk); + } // surf flux if it is specified as mean(theta*w) if(!params.ForceParameters.surf_sensible_flux_in_watts_per_square_meter)