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/drawbicyc.cpp b/drawbicyc/drawbicyc.cpp index b1f8812a25..3a0cc235bd 100644 --- a/drawbicyc/drawbicyc.cpp +++ b/drawbicyc/drawbicyc.cpp @@ -56,6 +56,7 @@ int main(int argc, char** argv) if(NDims == 2) { + if(flag_series) plot_series(PlotterMicro_t<2>(h5, micro), Plots(type)); if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), Plots(type)); if(flag_fields) plot_fields(PlotterMicro_t<2>(h5, micro), Plots(type)); @@ -68,8 +69,8 @@ int main(int argc, char** argv) if(flag_fields) plot_fields(PlotterMicro_t<3>(h5, micro), Plots(type)); if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro), Plots(type)); if(flag_lgrngn_spec) { - plot_lgrngn_spec_positions(PlotterMicro_t<3>(h5, "lgrngn"), Plots(type), 12); - plot_lgrngn_spec(PlotterMicro_t<3>(h5, "lgrngn"), Plots(type), 12); + plot_lgrngn_spec_positions(PlotterMicro_t<3>(h5, "lgrngn"), Plots(type)); + plot_lgrngn_spec(PlotterMicro_t<3>(h5, "lgrngn"), Plots(type)); } } else diff --git a/drawbicyc/plot_fields.hpp b/drawbicyc/plot_fields.hpp index 8f6100df34..1a4de8305b 100644 --- a/drawbicyc/plot_fields.hpp +++ b/drawbicyc/plot_fields.hpp @@ -31,7 +31,7 @@ void plot_fields(Plotter_t plotter, Plots plots) { std::cout << at * n["outfreq"] << " : " << plt << std::endl; Gnuplot gp; - init(gp, plotter.file + ".plot/" + plt + "/" + zeropad(at * n["outfreq"]) + ".svg", 1, 1, n); + init(gp, plotter.file + ".plot/" + plt + "/" + zeropad(at * n["outfreq"]) + ".svg", 1, 1, n, 1, 0.666666); if (plt == "rl") { @@ -205,6 +205,26 @@ 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 == "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/plot_lgrngn_spec.hpp b/drawbicyc/plot_lgrngn_spec.hpp index eb1e2b8a03..911bf19f20 100644 --- a/drawbicyc/plot_lgrngn_spec.hpp +++ b/drawbicyc/plot_lgrngn_spec.hpp @@ -4,21 +4,33 @@ #include "bins.hpp" #include +namespace +{ + po::variables_map vm; +}; //plot spectra positions template -void plot_lgrngn_spec_positions(Plotter_t plotter, Plots plots, int at) +void plot_lgrngn_spec_positions(Plotter_t plotter, Plots plots) { + // read opts + po::options_description opts("spectra plotting options"); + opts.add_options() + ("spectra_step", po::value()->required() , "time step number at which spectra are plotted") + ; + handle_opts(opts, vm); + const int spectra_step = vm["spectra_step"].as(); + auto& n = plotter.map; Gnuplot gp; try{ // cloud water content - auto tmp = plotter.h5load_ract_timestep(at * n["outfreq"]) * 1e3; + auto tmp = plotter.h5load_ract_timestep(spectra_step * n["outfreq"]) * 1e3; std::string title = "cloud water mixing ratio [g/kg]"; - gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n"; - init(gp, plotter.file + "_spectra_positions_at" + zeropad(at) + ".svg", 1, 1, n, 2.); + gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(spectra_step) * n["outfreq"] * n["dt"] / 60.) << "min'\n"; + init(gp, plotter.file + "_spectra_positions_at" + zeropad(spectra_step) + ".svg", 1, 1, n, 2.); { for (auto &fcs : focus_3d) @@ -60,8 +72,10 @@ void plot_lgrngn_spec_positions(Plotter_t plotter, Plots plots, int at) //plot spectra template -void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at) +void plot_lgrngn_spec(Plotter_t plotter, Plots plots) { + const int spectra_step = vm["spectra_step"].as(); + const double r_c_adiab = 7.8; // adiabatic water content [g/kg] coming from Wojtek auto& n = plotter.map; for(auto elem : n) @@ -71,10 +85,10 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at) Gnuplot gp; // int off = 2; // TODO!!! int off = 0; // TODO!!! - string file = plotter.file + "_spectra_at" + zeropad(at) + ".svg"; + string file = plotter.file + "_spectra_at" + zeropad(spectra_step) + ".svg"; - int hor = min(focus_3d.size(), 2); - int ver = double(focus_3d.size()) / 2. + 0.99999; + int hor = min(focus_3d.size(), 5); + int ver = double(focus_3d.size()) / 5. + 0.99999; init_prof(gp, file, ver, hor); @@ -82,7 +96,6 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at) gp << "set xrange [0:30]\n"; gp << "set yrange [0:50]\n"; gp << "set xlabel 'r [μm]'\n"; - gp << "set ylabel 'n [cm^{-3} μm^{-1}]'\n"; // read in density auto tmp = plotter.h5load(plotter.file + "/const.h5", "G"); @@ -90,19 +103,27 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at) // focus to the gridbox from where the size distribution is plotted char lbl = 'a'; + bool setylabel = 1; for (auto &fcs : focus_3d) { + if(setylabel) + { + gp << "set ylabel 'n [cm^{-3} μm^{-1}]'\n"; + setylabel = 0; + } + else + gp << "unset ylabel\n"; const int &x = fcs[0], &y = fcs[1], &z = fcs[2]; const blitz::RectDomain<3> focusBox({x-box_size, y-box_size, z-box_size}, {x+box_size, y+box_size, z+box_size}); - gp << "set label 1 '(" << lbl << ")' at graph .1, .93 font ',20'\n"; + gp << "set label 1 '(" << lbl << ")' at graph .1, .93 font \",15\"\n"; lbl += 1; // calc ratio of water content to adiabatic water content { - auto tmp = plotter.h5load_ract_timestep(at * n["outfreq"]) * 1e3; + auto tmp = plotter.h5load_ract_timestep(spectra_step * n["outfreq"]) * 1e3; double ratio = mean(tmp(focusBox)) / r_c_adiab; - gp << "set label 4 'r_c / r_c^{adiab} = " << std::setprecision(2) << ratio << "' at graph .2, .93 font ',20'\n"; + gp << "set label 4 'AF = " << std::fixed << std::setprecision(2) << ratio << "' at graph .2, .63 font \",15\"\n"; } // calc mean and std dev of radius of acivated droplets in the box @@ -111,7 +132,7 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at) double act_rw_std_dev = 0.; try { - auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]); + auto tmp = plotter.h5load_timestep("actrw_rw_mom0", spectra_step * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); act_conc = mean(snap(focusBox)); } @@ -119,7 +140,7 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at) // mean droplet radius in the box try { - auto tmp = plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]); + auto tmp = plotter.h5load_timestep("actrw_rw_mom1", spectra_step * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); // 1st raw moment / mass [m / kg] if(act_conc > 0) act_rw_mean = mean(snap(focusBox)) / act_conc; @@ -130,13 +151,13 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at) // std deviation of distribution of radius in the box try { - auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]); + auto tmp = plotter.h5load_timestep("actrw_rw_mom0", spectra_step * n["outfreq"]); typename Plotter_t::arr_t zeroth_raw_mom(tmp); // 0th raw moment / mass [1 / kg] - tmp = plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]); + tmp = plotter.h5load_timestep("actrw_rw_mom1", spectra_step * n["outfreq"]); typename Plotter_t::arr_t first_raw_mom(tmp); // 1st raw moment / mass [m / kg] - tmp = plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]); + tmp = plotter.h5load_timestep("actrw_rw_mom2", spectra_step * n["outfreq"]); typename Plotter_t::arr_t second_raw_mom(tmp); // 2nd raw moment / mass [m^2 / kg] - tmp = plotter.h5load_timestep("sd_conc", at * n["outfreq"]); + tmp = plotter.h5load_timestep("sd_conc", spectra_step * n["outfreq"]); typename Plotter_t::arr_t sd_conc(tmp); // number of SDs if(act_conc > 0) { @@ -165,8 +186,8 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at) catch(...) {;} - gp << "set label 2 '_{act} = " << std::setprecision(2) << act_rw_mean * 1e6 <<"µm' at graph .1, .83 font '15'\n"; - gp << "set label 3 'σ(r)_{act} = " << std::setprecision(2) << act_rw_std_dev * 1e6 <<"µm' at graph .1, .73 font ',15'\n"; + gp << "set label 2 ' = " << std::fixed << std::setprecision(2) << act_rw_mean * 1e6 <<" µm' at graph .2, .83 font \",15\"\n"; + gp << "set label 3 'σ = " << std::fixed << std::setprecision(2) << act_rw_std_dev * 1e6 <<" µm' at graph .2, .73 font \",15\"\n"; std::map focus_w; @@ -178,7 +199,7 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at) for (int i = 0; i < nsw; ++i) { const string name = "rw_rng" + zeropad(i + off) + "_mom0"; - auto tmp_w = plotter.h5load_timestep(name, at * n["outfreq"]); + auto tmp_w = plotter.h5load_timestep(name, spectra_step * n["outfreq"]); focus_w[(left_edges_rw[i] + left_edges_rw[i+1]) / 2. / 1e-6 / si::metres] = sum((tmp_w*rhod)(focusBox)) @@ -187,7 +208,7 @@ void plot_lgrngn_spec(Plotter_t plotter, Plots plots, int at) / ((left_edges_rw[i+1] - left_edges_rw[i]) / 1e-6 / si::metres); // per micrometre } const string name = "rw_rng" + zeropad(nsw + off) + "_mom0"; - auto tmp_w = plotter.h5load_timestep(name, at * n["outfreq"]); + auto tmp_w = plotter.h5load_timestep(name, spectra_step * n["outfreq"]); notice_macro("setting-up plot parameters"); std::cout << "larger drops conc: " << (tmp_w * rhod)(x,y,z) * 1e-6 << endl; diff --git a/drawbicyc/plot_prof.hpp b/drawbicyc/plot_prof.hpp index 17fe40a6bc..a709eaee79 100644 --- a/drawbicyc/plot_prof.hpp +++ b/drawbicyc/plot_prof.hpp @@ -41,7 +41,7 @@ void plot_profiles(Plotter_t plotter, Plots plots) int last_timestep = vm["prof_end"].as() / n["dt"] / n["outfreq"]; // some ugly constants - const double p_1000 = 1000.; + const double p_1000 = 100000.; const double L = 2.5e6; const double R_d = 287.; const double c_p = 1004; @@ -49,6 +49,8 @@ void plot_profiles(Plotter_t plotter, Plots plots) double z_i; + bool res_pos_out_done = false; + for (auto &plt : plots.profs) { blitz::firstIndex i; @@ -69,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") @@ -330,49 +332,39 @@ 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; + res_tmp = (snap -1) * 100; } +/* + // for mean over updraught cells { auto tmp = plotter.h5load_timestep("w", at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); res_tmp2 = isupdraught(snap); res_tmp *= res_tmp2; } +*/ // mean only over updraught cells - res_pos = plotter.horizontal_sum(res_tmp2); // number of downdraft cells on a given level - res_prof += where(res_pos > 0 , plotter.horizontal_sum(res_tmp) / res_pos, 0); +// res_pos = plotter.horizontal_sum(res_tmp2); // number of downdraft cells on a given level +// res_prof += where(res_pos > 0 , plotter.horizontal_sum(res_tmp) / res_pos, 0); + +// mean over all cells +// res_prof += plotter.horizontal_sum(res_tmp); res += res_tmp; - gp << "set title 'supersaturation RH-based in updrafts only'\n"; - gp << "set yrange [0.45:1.]\n"; - gp << "set xrange [0.000:*]\n"; + //gp << "set title 'supersaturation RH-based in updrafts only'\n"; + gp << "set title 'supersaturation RH-based'\n"; +// gp << "set yrange [0.45:1.]\n"; + // gp << "set xrange [0.000:*]\n"; } else if (plt == "00rtot") { - // total water content (vapor + cloud + rain, missing droplets with r<0.5um!) - /* - { - 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 @@ -383,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 @@ -391,27 +383,50 @@ void plot_profiles(Plotter_t plotter, Plots plots) } gp << "set title 'cloud droplets ( 0.5um < r < 25um) concentration [1/cm^3]'\n"; } + else if (plt == "cl_nc") + { + // cloud droplet (0.5um < r < 25 um) concentration in cloudy grid cells + try + { + // cloud fraction (cloudy if N_c > 20/cm^3) + 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 + typename Plotter_t::arr_t snap2; + snap2.resize(snap.shape()); + snap2=snap; + snap = iscloudy(snap); // cloudiness mask + snap2 *= snap; + + // mean only over cloudy cells + res_pos = plotter.horizontal_sum(snap); // number of cloudy cells on a given level + res_prof += where(res_pos > 0 , plotter.horizontal_sum(snap2) / res_pos, 0); + } + catch(...){;} + } else if (plt == "thl") { // liquid potential temp [K] { - { - auto tmp = plotter.h5load_timestep("cloud_rw_mom3", at * n["outfreq"]) * 4./3 * 3.14 * 1e3; - typename Plotter_t::arr_t snap(tmp); - res_tmp2 = snap; - } - { - auto tmp = plotter.h5load_timestep("rain_rw_mom3", at * n["outfreq"]) * 4./3 * 3.14 * 1e3; - typename Plotter_t::arr_t snap(tmp); - res_tmp2 += snap; - } - // res_tmp2 is now q_l (liq water content) - auto tmp = plotter.h5load_timestep("th", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); // snap is theta_dry - res_tmp = pow(snap * pow(rhod * R_d / (p_1000 * 100), R_d / c_pd), c_pd / (c_pd - R_d)); // res_tmp is now temperature; 1 bar = 100 000Pa - snap *= (res_tmp - res_tmp2 * L / c_p) / res_tmp; - res += snap; -// res += res_tmp2; + 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"])); + ///auto tmp = plotter.h5load_timestep("rv", at * n["outfreq"]); + typename Plotter_t::arr_t rv(plotter.h5load_timestep("rv", 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)); + // 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); +// res += ql; } gp << "set title 'liquid potential temp [K]'\n"; } @@ -419,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 @@ -432,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) /* @@ -491,17 +501,23 @@ void plot_profiles(Plotter_t plotter, Plots plots) z_i = (double(k_i)-0.5) / (last_timestep - first_timestep + 1) * n["dz"]; std::cout << "average inversion height " << z_i; res_pos = i * n["dz"] / z_i; + if(!res_pos_out_done) + { + oprof_file << res_pos; + res_pos_out_done = true; + } + if (plt != "act_rd" && plt != "act_conc") { - if (plt == "ugccn_rw_down" || plt == "sat_RH" || plt=="gccn_rw_down" || plt=="non_gccn_rw_down" || plt=="gccn_rw_up" || plt=="non_gccn_rw_up") + if (plt == "ugccn_rw_down" || /*plt == "sat_RH" ||*/ plt=="gccn_rw_down" || plt=="non_gccn_rw_down" || plt=="gccn_rw_up" || plt=="non_gccn_rw_up" || plt == "cl_nc") // these are plots that are done only in up/downdrafts/cloudy cells (sat_RH now calculated over all cells) res_prof /= last_timestep - first_timestep + 1; else res_prof = plotter.horizontal_mean(res); // average in x // res_prof(0) is ground level - we dont know what is there? surf fluxes shouldnt be added to it?! anyway, set res_prof(0)=res_prof(1) for plotting purposes - res_prof(0) = res_prof(1); +// res_prof(0) = res_prof(1); gp << "plot '-' with line\n"; gp.send1d(boost::make_tuple(res_prof, res_pos)); @@ -520,8 +536,8 @@ void plot_profiles(Plotter_t plotter, Plots plots) res_prof /= last_timestep - first_timestep + 1; res_prof2 /= last_timestep - first_timestep + 1; // res_prof(0) is ground level - we dont know what is there? surf fluxes shouldnt be added to it?! anyway, set res_prof(0)=res_prof(1) for plotting purposes - res_prof(0) = res_prof(1); - res_prof2(0) = res_prof2(1); +// res_prof(0) = res_prof(1); + // res_prof2(0) = res_prof2(1); gp.send1d(boost::make_tuple(res_prof, res_pos)); gp.send1d(boost::make_tuple(res_prof2, res_pos)); oprof_file << res_prof ; diff --git a/drawbicyc/plot_series.hpp b/drawbicyc/plot_series.hpp index 7ff1431be1..a77a6898d5 100644 --- a/drawbicyc/plot_series.hpp +++ b/drawbicyc/plot_series.hpp @@ -83,6 +83,7 @@ void plot_series(Plotter_t plotter, Plots plots) if (plt == "clfrac") { +/* try { // cloud fraction (cloudy if q_c > 0.1 g/kg) @@ -92,6 +93,16 @@ void plot_series(Plotter_t plotter, Plots plots) res_tmp = iscloudy_rc(snap); // find cells with rc>1e-5 res_prof(at) = blitz::mean(res_tmp); } +*/ + try + { + auto tmp = plotter.h5load_ract_timestep(at * n["outfreq"]) * 1e3; //g/kg + typename Plotter_t::arr_t snap(tmp); + snap *= rhod; // water per cubic metre (should be wet density...) + plotter.k_i = blitz::sum(snap, plotter.LastIndex) * n["dz"]; // LWP [g/m2] in the column + plotter.k_i = where(plotter.k_i > 20 , 1 , 0); // cloudiness as in Ackermann et al. + res_prof(at) = blitz::mean(plotter.k_i); + } catch(...){;} } // max RH in the domain diff --git a/drawbicyc/plots.hpp b/drawbicyc/plots.hpp index 2482ee8b26..8d53484a80 100644 --- a/drawbicyc/plots.hpp +++ b/drawbicyc/plots.hpp @@ -37,7 +37,8 @@ std::vector profs_dycoms({ //, "act_conc" ,"clfrac" //, "N_c", -//"sat_RH" +,"cl_nc" +,"sat_RH" }); // rtot has to be first std::vector profs_moist_thermal({ @@ -51,7 +52,8 @@ std::vector fields_dycoms({ "th", "rv", "u", "w", "sd_conc",//, "r_dry", -"RH", "supersat" +"RH", "supersat", +"lib_pres", "lib_temp" }); std::vector fields_moist_thermal({ diff --git a/src/bicycles.cpp b/src/bicycles.cpp index 195cd7614f..efe8670e20 100644 --- a/src/bicycles.cpp +++ b/src/bicycles.cpp @@ -17,22 +17,41 @@ #include "opts_lgrngn.hpp" #include "opts_blk_1m.hpp" + #include "panic.hpp" #include // 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 &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, const int halo) { - 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.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); - p.rhod = new setup::arr_1D_t(rhod.dataFirst(), rhod.shape(), blitz::neverDeleteData); - p.w_LS = new setup::arr_1D_t(w_LS.dataFirst(), w_LS.shape(), blitz::neverDeleteData); + const auto nz = p.grid_size.back(); + + std::vector, std::reference_wrapper>> tobecopied = { + {p.hgt_fctr_sclr, hgt_s }, + {p.hgt_fctr_vctr, hgt_v }, + {p.th_e , th_e }, + {p.p_e , p_e }, + {p.rv_e , rv_e }, + {p.th_ref , th_ref }, + {p.pre_ref , pre_ref}, + {p.rhod , rhod }, + {p.w_LS , w_LS } + }; + + for (auto dst_src : tobecopied) + { + auto a = new setup::arr_1D_t(blitz::Range(-halo, nz - 1 + halo)); + dst_src.first.get() = new setup::arr_1D_t(a->dataFirst(), a->shape(), blitz::neverDeleteData); + dst_src.first.get()->reindexSelf(a->base()); + + // copy profile data + for (int k = 0; k < nz; ++k) + { + (*dst_src.first.get())(k) = dst_src.second.get()(k); + } + } } // 2D model run logic - the same for any microphysics @@ -102,13 +121,14 @@ 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), p_e(nz), rv_e(nz), th_ref(nz), pre_ref(nz), rhod(nz), 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, 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, solver_t::halo); + // set outvars p.outvars.insert({solver_t::ix::rv, {"rv", "[kg kg-1]"}}); @@ -122,17 +142,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 @@ -213,13 +233,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, 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); + copy_profiles(th_e, p_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_vctr, hgt_fctr_sclr, p, solver_t::halo); // set outvars p.outvars.insert({solver_t::ix::rv, {"rv", "[kg kg-1]"}}); @@ -235,17 +255,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 @@ -338,6 +358,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 @@ -348,6 +369,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 @@ -424,13 +446,15 @@ int main(int argc, char** argv) ("nt", po::value()->default_value(3600) , "timestep count") ("rng_seed", po::value()->default_value(-1) , "rng seed, negative for random") ("dt", po::value()->required() , "timestep length") - ("z_rlx_sclr", po::value()->default_value(10) , "scalars surface flux charasteristic heihjt") + ("z_rlx_sclr", po::value()->default_value(25) , "scalars surface flux charasteristic heihjt") ("outdir", po::value(), "output file name (netCDF-compatible HDF5)") ("outfreq", po::value(), "output rate (timestep interval)") ("spinup", po::value()->default_value(2400) , "number of initial timesteps during which rain formation is to be turned off") ("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") @@ -489,6 +513,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(); @@ -513,7 +539,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/calc_forces_blk_1m.hpp b/src/calc_forces_blk_1m.hpp new file mode 100644 index 0000000000..6bd248e3bd --- /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) += (*this->mem->vab_coeff)(ijk) * this->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) += (*this->mem->vab_coeff)(ijk) * this->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.hpp b/src/calc_forces_common.hpp similarity index 68% rename from src/calc_forces.hpp rename to src/calc_forces_common.hpp index c68d8e8b97..a68d8bb401 100644 --- a/src/calc_forces.hpp +++ b/src/calc_forces_common.hpp @@ -17,7 +17,14 @@ struct calc_T BZ_DECLARE_FUNCTOR2(calc_T) }; -// forcing functions +struct calc_exner +{ + setup::real_t operator()(setup::real_t p) const + {return libcloudphxx::common::theta_std::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() @@ -28,11 +35,14 @@ void slvr_common::rv_src() // surface flux surf_latent(); // sum of rv flux - this->vert_grad_fwd(F, alpha, params.dz); + //this->vert_grad_fwd(F, alpha, params.dz); + alpha(ijk) = F(ijk); // 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) /= (libcloudphxx::common::const_cp::l_tri() * si::kilograms / si::joules) * rhod(this->vert_idx); +// else + // alpha(ijk) *= -1; // negative gradient means inflow // 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); @@ -45,8 +55,8 @@ void slvr_common::rv_src() 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); + //alpha(ijk) += (*this->mem->vab_coeff)(ijk) * rv_e(this->vert_idx); // TODO: its a constant, cache it + //beta(ijk) = - (*this->mem->vab_coeff)(ijk); } template @@ -57,6 +67,7 @@ void slvr_common::th_src(typename parent_t::arr_t &rv) { // -- heating -- // surface flux + /* if(params.ForceParameters.surf_sensible_flux_in_watts_per_square_meter) { surf_sens(); @@ -66,33 +77,41 @@ void slvr_common::th_src(typename parent_t::arr_t &rv) } 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"); +// 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"); + alpha(ijk) *= -1; // negative gradient means inflow + nancheck(alpha(ijk), "sum of th flux"); + + if(params.ForceParameters.surf_sensible_flux_in_watts_per_square_meter) + { + surf_sens(); + nancheck(F(ijk), "sensible surf forcing"); + // beta as tmp storage + alpha(ijk) += F(ijk); + } - // 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) / - calc_c_p()(rv(ijk).reindex(this->zero)) / - calc_T()(this->state(ix::th)(ijk).reindex(this->zero), (*params.rhod)(this->vert_idx)) / - (*params.rhod)(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) /= calc_exner()(p_e(this->vert_idx)) * calc_c_p()(rv(ijk)) * rhod(this->vert_idx); - nancheck2(alpha(ijk), this->state(ix::th)(ijk), "change of theta"); + 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); +// this->vert_grad_fwd(F, beta, params.dz); + // nancheck(beta(ijk), "grad of sensible surf forcing"); + // alpha(ijk) += - beta(ijk); // negative gradient of upward flux means inflow + alpha(ijk) += F(ijk); } // large-scale vertical wind @@ -108,7 +127,7 @@ void slvr_common::th_src(typename parent_t::arr_t &rv) 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); + //alpha(ijk) += (*this->mem->vab_coeff)(ijk) * th_e(this->vert_idx); //beta(ijk) = - (*this->mem->vab_coeff)(ijk); } @@ -126,5 +145,3 @@ void slvr_common::w_src(typename parent_t::arr_t &th, typename pare alpha(ijk) += F(ijk); } - - diff --git a/src/cases/CasesCommon.hpp b/src/cases/CasesCommon.hpp index 0eb057f1b1..b518a5ca8a 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 @@ -63,8 +63,8 @@ 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 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 214d9f05ec..3f1f250411 100644 --- a/src/cases/DYCOMS98.hpp +++ b/src/cases/DYCOMS98.hpp @@ -21,7 +21,7 @@ namespace setup Y = 6400 * si::metres; // DYCOMS: 6400 const real_t z_abs = 1250; const real_t z_i = 795; //initial inversion height - const quantity z_rlx_vctr = 1 * si::metres; + const quantity z_rlx_vctr = 25 * si::metres; // liquid water potential temperature at height z quantity th_l(const real_t &z) @@ -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; @@ -167,9 +169,10 @@ 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) + // alse set w_LS and hgt_fctrs + // TODO: move hgt_fctrs from cases to main code + 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,16 +180,14 @@ namespace setup using libcloudphxx::common::const_cp::l_tri; using libcloudphxx::common::theta_std::p_1000; - // pressure profile - arr_1D_t pre(nz); - // temperature profile + // 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(0) = th_l(0.) / si::kelvins * pow(p_0 / p_1000(), R_d_over_c_pd()); - pre(0) = p_0 / si::pascals; + p_e(0) = p_0 / si::pascals; th_e(0) = th_l(0.) / si::kelvins; rv_e(0) = rt(0.); @@ -198,26 +199,35 @@ 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::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 = 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 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(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); + + 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; @@ -229,6 +239,8 @@ 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) // 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; @@ -236,14 +248,20 @@ 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 - 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? + 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 + + // 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; @@ -253,6 +271,7 @@ namespace setup // surface sources relaxation factors // for vectors + /* real_t z_0 = 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; @@ -260,6 +279,14 @@ namespace setup 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; } void update_surf_flux_sens(typename concurr_t::solver_t::arr_sub_t &surf_flux_sens, int timestep, real_t dt) @@ -299,7 +326,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); @@ -320,7 +347,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/DryThermalGMD2015.hpp b/src/cases/DryThermalGMD2015.hpp index 9a314b3823..83672b1426 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; @@ -76,10 +78,20 @@ 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 &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; + //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 th_ref = 300; } }; @@ -98,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); @@ -121,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 439d22ead4..9adbf8adfb 100644 --- a/src/cases/LasherTrapp2001.hpp +++ b/src/cases/LasherTrapp2001.hpp @@ -29,11 +29,12 @@ namespace setup X = 10000 * si::metres, // DYCOMS: 6400 Y = 10000 * si::metres; // DYCOMS: 6400 const real_t z_abs = 7000; - const quantity z_rlx_vctr = 1 * si::metres; + const quantity z_rlx_vctr = 25 * si::metres; // env profiles of th and rv from the sounding arr_1D_t th_dry_env; arr_1D_t th_std_env; + arr_1D_t p_env; arr_1D_t rv_env; @@ -59,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; @@ -116,7 +119,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 &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; @@ -179,13 +182,15 @@ 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); + 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(); + 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; th_e = th_std_env; // temp to calc rhod @@ -217,6 +222,7 @@ namespace setup // surface sources relaxation factors // for vectors + /* real_t z_0 = 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; @@ -224,6 +230,15 @@ namespace setup 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; } // functions that set surface fluxes per timestep @@ -289,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); @@ -310,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 a8517c00bb..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 @@ -166,6 +145,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; @@ -184,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; @@ -199,14 +178,13 @@ namespace setup // initial potential temperature solver.advectee(ix::th) = th_e(index); -//libcloudphxx::common::theta_dry::std2dry(th_std_0, rv_0) / si::kelvins; } 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 &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 @@ -246,10 +224,11 @@ 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; + p_e(0) = pre_ref(0); T(0.) = T_0 / si::kelvins; for(int k=1; 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); } @@ -315,39 +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; - } - */ - }; // 2d/3d children @@ -365,24 +312,27 @@ 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) ), blitz::tensor::j * dz ); - + } }; @@ -401,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); @@ -424,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/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/forcings.hpp b/src/forcings.hpp index 75f350c339..cbd832f374 100644 --- a/src/forcings.hpp +++ b/src/forcings.hpp @@ -12,19 +12,20 @@ void slvr_common::buoyancy(typename parent_t::arr_t &th, typename p namespace moist_air = libcloudphxx::common::moist_air; const real_t eps = moist_air::R_v() / moist_air::R_d() - 1.; if(params.buoyancy_wet) - tmp1(ijk).reindex(this->zero) = + tmp1(ijk) = (libcloudphxx::common::earth::g() / si::metres_per_second_squared) * ( - (th(ijk).reindex(this->zero) - (*params.th_e)(this->vert_idx)) / (*params.th_ref)(this->vert_idx) - + eps * (rv(ijk).reindex(this->zero) - (*params.rv_e)(this->vert_idx)) - - r_l(ijk).reindex(this->zero) + (th(ijk) - th_e(this->vert_idx)) / th_ref(this->vert_idx) + + eps * (rv(ijk) - rv_e(this->vert_idx)) + - r_l(ijk) ); else - tmp1(ijk).reindex(this->zero) = + tmp1(ijk) = (libcloudphxx::common::earth::g() / si::metres_per_second_squared) * ( - (th(ijk).reindex(this->zero) - (*params.th_e)(this->vert_idx)) / (*params.th_ref)(this->vert_idx) + (th(ijk) - th_e(this->vert_idx)) / th_ref(this->vert_idx) ); - this->smooth(tmp1, F); +// this->smooth(tmp1, F); + F(ijk) = tmp1(ijk); } template @@ -43,13 +44,13 @@ void slvr_common::radiation(typename parent_t::arr_t &rv) // index of first cell above inversion tmp1(ijk) = rv(ijk) + r_l(ijk); - k_i(this->hrzntl_subdomain) = blitz::first( tmp1(ijk).reindex(this->zero) < params.ForceParameters.q_i, this->vert_idx); // vertical index of first cell above inversion (inversion is at the lower edge of this cell) + k_i(this->hrzntl_subdomain) = blitz::first( tmp1(ijk) < params.ForceParameters.q_i, this->vert_idx); // vertical index of first cell above inversion (inversion is at the lower edge of this cell) // calc Eqs. 5 and 6 from Ackerman et al 2009 // calc sum of r_l above certain level and store it in tmp1 tmp1(ijk) = r_l(ijk); - tmp1(ijk).reindex(this->zero) *= params.ForceParameters.heating_kappa * (*params.rhod)(this->vert_idx); + tmp1(ijk) *= params.ForceParameters.heating_kappa * rhod(this->vert_idx); for(int z = nz-2 ; z >= 0; --z) tmp1(idxperm::pi(z, this->hrzntl_subdomain)) += tmp1(idxperm::pi(z+1, this->hrzntl_subdomain)); @@ -60,14 +61,14 @@ void slvr_common::radiation(typename parent_t::arr_t &rv) // multiply by distance from the bottom of the cell to the top of the domain tmp1(ground) *= - (nz - 1) * params.dz * tmp1(ground); - tmp1(noground).reindex(this->zero) *= - (nz - this->vert_idx - 1.5) * params.dz; // vert_idx starts from 0, but its 2nd cell from ground; do not merge this line with F_0 * exp(...) since it gives some strange values! + tmp1(noground) *= - (nz - this->vert_idx - 1.5) * params.dz; // vert_idx starts from 0, but its 2nd cell from ground; do not merge this line with F_0 * exp(...) since it gives some strange values! F(ijk) = params.ForceParameters.F_0 * exp(tmp1(ijk)); // calc sum of r_l below certain level and store it in tmp1 tmp1(ijk) = r_l(ijk); - tmp1(ijk).reindex(this->zero) *= params.ForceParameters.heating_kappa * (*params.rhod)(this->vert_idx); + tmp1(ijk) *= params.ForceParameters.heating_kappa * rhod(this->vert_idx); // copy one cell upwards for(int z = nz-1 ; z >= 1; --z) @@ -78,12 +79,12 @@ void slvr_common::radiation(typename parent_t::arr_t &rv) tmp1(idxperm::pi(z, this->hrzntl_subdomain)) += tmp1(idxperm::pi(z-1, this->hrzntl_subdomain)); // multiply by distance from the bottom of the cell to the bottom of the domain - tmp1(noground).reindex(this->zero) *= - (this->vert_idx + 0.5) * params.dz; + tmp1(noground) *= - (this->vert_idx + 0.5) * params.dz; tmp1(ground) = 0.; F(ijk) += params.ForceParameters.F_1 * exp(tmp1(ijk)); // free atmosphere part - F(ijk).reindex(this->zero) += where(this->vert_idx > k_i(this->hrzntl_subdomain)(blitz::tensor::i, blitz::tensor::j), // works even in 2D ?!?! + F(ijk) += where(this->vert_idx > k_i(this->hrzntl_subdomain)(blitz::tensor::i, blitz::tensor::j), // works even in 2D ?!?! (libcloudphxx::common::moist_air::c_pd() / si::joules * si::kilograms * si::kelvins) * params.ForceParameters.rho_i * params.ForceParameters.D * (0.25 * pow((this->vert_idx - 0.5) * params.dz - (k_i(this->hrzntl_subdomain)(blitz::tensor::i, blitz::tensor::j) - .5) * params.dz, 4./3) + (k_i(this->hrzntl_subdomain)(blitz::tensor::i, blitz::tensor::j) - .5) * params.dz * pow((this->vert_idx - 0.5) * params.dz - (k_i(this->hrzntl_subdomain)(blitz::tensor::i, blitz::tensor::j) - .5) * params.dz, 1./3)) @@ -102,8 +103,8 @@ void slvr_common::surf_sens() //TODO: each thread has surf_flux_sens of the size of the domain of all threads and each updates all of it // either make it shared among threads and updated by one all make it of the size of hrzntl_subdomain params.update_surf_flux_sens(surf_flux_sens, this->timestep, this->dt); - F(ijk).reindex(this->zero) = - surf_flux_sens(this->hrzntl_subdomain)(blitz::tensor::i, blitz::tensor::j) // "-" because negative gradient means inflow - * (*params.hgt_fctr_sclr)(this->vert_idx); + F(ijk) = surf_flux_sens(this->hrzntl_subdomain)(blitz::tensor::i, blitz::tensor::j) + * hgt_fctr_sclr(this->vert_idx); // tmp1(ijk)=F(ijk); //TODO: unnecessary copy //this->smooth(tmp1, F); @@ -116,8 +117,8 @@ void slvr_common::surf_latent() //TODO: each thread has surf_flux_sens of the size of the domain of all threads and each updates all of it // either make it shared among threads and updated by one all make it of the size of hrzntl_subdomain params.update_surf_flux_lat(surf_flux_lat, this->timestep, this->dt); - F(ijk).reindex(this->zero) = - surf_flux_lat(this->hrzntl_subdomain)(blitz::tensor::i, blitz::tensor::j) - * (*params.hgt_fctr_sclr)(this->vert_idx); + F(ijk) = surf_flux_lat(this->hrzntl_subdomain)(blitz::tensor::i, blitz::tensor::j) + * hgt_fctr_sclr(this->vert_idx); // tmp1(ijk)=F(ijk); //TODO: unnecessary copy //this->smooth(tmp1, F); @@ -131,10 +132,10 @@ void slvr_common::subsidence(const int &type) // large-scale vertic { tmp1(ijk) = this->state(type)(ijk); this->vert_grad_cnt(tmp1, F, params.dz); - F(ijk).reindex(this->zero) *= - (*params.w_LS)(this->vert_idx); + F(ijk) *= - 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.; 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..9ae6f986ff 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 @@ -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; diff --git a/src/slvr_blk_1m.hpp b/src/slvr_blk_1m.hpp index c2b272c8d0..2911f63c4d 100644 --- a/src/slvr_blk_1m.hpp +++ b/src/slvr_blk_1m.hpp @@ -6,6 +6,7 @@ #include #include + template class slvr_blk_1m_common : public slvr_common { @@ -15,8 +16,12 @@ 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 array with copy of the environmental total pressure of dry air + typename parent_t::arr_t &p_e_nd; + void condevap() { auto @@ -25,11 +30,21 @@ 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_nd(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, th, rv, rc, rr, this->dt + ); +*/ + + libcloudphxx::blk_1m::adj_cellwise_nwtrph( + opts, p_e_arg, th, rv, rc, this->dt + ); this->mem->barrier(); } @@ -40,16 +55,52 @@ 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) { // if uninitialised fill with zeros zero_if_uninitialised(ix::rc); zero_if_uninitialised(ix::rr); + + // init the p_e array + p_e_nd(this->ijk) = this->p_e(this->vert_idx); // deal with initial supersaturation condevap(); @@ -69,14 +120,24 @@ 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); } } void hook_ante_step() { parent_t::hook_ante_step(); + + 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 - 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( @@ -84,24 +145,78 @@ 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? + this->mem->barrier(); + if(this->rank == 0) + this->tbeg = clock::now(); + // cell-wise + // 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_nd(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 + 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 ); } } // 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 } @@ -114,13 +229,30 @@ class slvr_blk_1m_common : public slvr_common libcloudphxx::blk_1m::opts_t cloudph_opts; }; + protected: + + // per-thread copy of params + // TODO: but slvr_common also has a copy of it's params.... + rt_params_t params; + + 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_nd + } + // ctor slvr_blk_1m_common( typename parent_t::ctor_args_t args, const rt_params_t &p ) : parent_t(args, p), - opts(p.cloudph_opts) + params(p), + opts(p.cloudph_opts), + puddle(0), + p_e_nd(args.mem->tmp[__FILE__][0][0]) {} }; @@ -141,6 +273,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( @@ -158,6 +291,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) { @@ -166,7 +303,15 @@ 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(); + if(this->rank == 0) + { + 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 ); } } }; @@ -181,6 +326,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( @@ -198,6 +344,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) @@ -207,7 +357,15 @@ 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(); + if(this->rank == 0) + { + 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 ); + } } }; diff --git a/src/slvr_common.hpp b/src/slvr_common.hpp index e9dc4d9f4a..c53beb7775 100644 --- a/src/slvr_common.hpp +++ b/src/slvr_common.hpp @@ -6,7 +6,6 @@ #include #include "../git_revision.h" - template class slvr_common : public slvr_dim { @@ -20,8 +19,11 @@ 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 - clock::time_point tbeg, tend, tbeg1, tbeg_loop; - std::chrono::milliseconds tdiag, tupdate, tsync, tasync, tasync_wait, tloop, tvip_rhs; + // 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, tbeg_loop; + std::chrono::milliseconds tdiag, tupdate, tsync, tsync_wait, tasync, tasync_wait, tloop, tvip_rhs, tnondelayed_step; int spinup; // number of timesteps @@ -40,6 +42,8 @@ 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 + setup::arr_1D_t th_e, p_e, rv_e, th_ref, pre_ref, rhod, w_LS, hgt_fctr_sclr, hgt_fctr_vctr; + // surface precip stuff std::ofstream f_puddle; // output precipitation file @@ -257,19 +261,43 @@ class slvr_common : public slvr_dim // loop over horizontal dimensions for(int it = 0; it < parent_t::n_dims-1; ++it) { - F(this->ijk).reindex(this->zero) = + this->vip_rhs[it](this->ijk) += + where(U_ground(blitz::tensor::i, blitz::tensor::j) == 0., 0., + -2 * pow(params.ForceParameters.u_fric,2) * // const, cache it + this->vip_ground[it](blitz::tensor::i, blitz::tensor::j) / // u_i at z=0 + U_ground(blitz::tensor::i, blitz::tensor::j) * // |U| at z=0 + hgt_fctr_vctr(this->vert_idx) // hgt_fctr + ); +/* + * + this->vip_rhs[it](this->ijk).reindex(this->zero) += where(U_ground(blitz::tensor::i, blitz::tensor::j) == 0., 0., - -pow(params.ForceParameters.u_fric,2) * // const, cache it + -2 * pow(params.ForceParameters.u_fric,2) * // const, cache it this->vip_ground[it](blitz::tensor::i, blitz::tensor::j) / // u_i at z=0 U_ground(blitz::tensor::i, blitz::tensor::j) * // |U| at z=0 (*params.hgt_fctr_vctr)(this->vert_idx) // hgt_fctr ); +*/ // du/dt = sum of kinematic momentum fluxes * dt - this->vert_grad_fwd(F, this->vip_rhs[it], params.dz); +// this->vert_grad_fwd(F, this->vip_rhs[it], params.dz); // multiplied by 2 here because it is later multiplied by 0.5 * dt - this->vip_rhs[it](this->ijk) *= -2; + // this->vip_rhs[it](this->ijk) *= -2; + // } + +/* + for (int j = this->j.first(); j <= this->j.last(); ++j) + { + this->vip_rhs[0](this->i, j).reindex({0}) += where(U_ground(blitz::tensor::i) == 0., 0., + -2 * pow(params.ForceParameters.u_fric,2) * // const, cache it + this->vip_ground[0](blitz::tensor::i) / // u_i at z=0 + U_ground(blitz::tensor::i) * // |U| at z=0 + (*params.hgt_fctr_vctr)(this->vert_idx) // hgt_fctr + ); + } +*/ + this->mem->barrier(); if(this->rank == 0) { @@ -298,8 +326,10 @@ 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; + << "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; } } } @@ -311,7 +341,8 @@ 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; + 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; 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 @@ -337,7 +368,16 @@ 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]), + th_e((*p.th_e)(this->vert_rng)), + p_e((*p.p_e)(this->vert_rng)), + rv_e((*p.rv_e)(this->vert_rng)), + th_ref((*p.th_ref)(this->vert_rng)), + pre_ref((*p.pre_ref)(this->vert_rng)), + rhod((*p.rhod)(this->vert_rng)), + w_LS((*p.w_LS)(this->vert_rng)), + hgt_fctr_sclr((*p.hgt_fctr_sclr)(this->vert_rng)), + hgt_fctr_vctr((*p.hgt_fctr_vctr)(this->vert_rng)) { 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_dim.hpp b/src/slvr_dim.hpp index e3c9d136d8..5b81def4ff 100644 --- a/src/slvr_dim.hpp +++ b/src/slvr_dim.hpp @@ -47,16 +47,34 @@ class slvr_dim< blitz::TinyVector zero = blitz::TinyVector({0,0}); blitz::secondIndex vert_idx; + const rng_t &vert_rng = this->j; libmpdataxx::arrvec_t vip_ground; std::set hori_vel = std::set{ix::u}; + void GCtoC(arrvec_t &C) + { + using namespace arakawa_c; + const auto &i(this->i), &j(this->j); + const auto &im(this->im), &jm(this->jm); + C[0](im + h, j) = this->mem->GC[0](im + h, j) / (0.5 * ((*this->mem->G)(im, j) + (*this->mem->G)(im + 1, j))); + C[1](i, jm + h) = this->mem->GC[1](i, jm + h) / (0.5 * ((*this->mem->G)(i, jm) + (*this->mem->G)(i, jm + 1))); + } + void vert_grad_fwd(typename parent_t::arr_t &in, typename parent_t::arr_t &out, setup::real_t dz) { - in(this->i, this->j.last() + 1) = in(this->i, this->j.last()); + // extrapolate upward, top cell is two times lower + in(this->i, this->j.last() + 1) = 1.5*in(this->i, this->j.last()) - .5 * in(this->i, this->j.last()-1); out(this->i, this->j) = ( in(this->i, this->j+1) - in(this->i, this->j)) / dz; - // top and bottom cells are two times lower - out(this->i, 0) *= 2; + // top nad bottom cells are two times lower out(this->i, this->j.last()) *= 2; + out(this->i, 0) *= 2; + + // we don't want surf fluxes to change values at the ground ? (j=0) + /* + out(this->i, 0) = 0; + // instead, the flux from 0-th level goes to the first level above the gfround + out(this->i, 1) = ( in(this->i, 2) - in(this->i, 0)) / (1.5*dz); + */ } void vert_grad_cnt(typename parent_t::arr_t &in, typename parent_t::arr_t &out, setup::real_t dz) @@ -66,7 +84,9 @@ class slvr_dim< out(this->i, this->j) = ( in(this->i, this->j+1) - in(this->i, this->j-1)) / 2./ dz; // top and bottom cells are two times lower out(this->i, 0) *= 2; - out(this->i, this->j.last()) *= 2; + //out(this->i, this->j.last()) *= 2; + // set to 0 at top level to have no subsidence there - TODO: it messes with other possible uses of this function + out(this->i, this->j.last()) = 0; } void smooth(typename parent_t::arr_t &in, typename parent_t::arr_t &out) @@ -82,7 +102,7 @@ class slvr_dim< auto calc_U_ground() return_macro(, - abs(this->state(ix::vip_i)(this->i, 0).reindex({0})) + abs(this->state(ix::vip_i)(this->i, 0)) ) // ctor @@ -92,7 +112,7 @@ class slvr_dim< ) : parent_t(args, p) { - vip_ground.push_back(new arr_sub_t(this->state(ix::vip_i)(this->i, 0).reindex({0}))); + vip_ground.push_back(new arr_sub_t(this->state(ix::vip_i)(this->i, 0))); } }; @@ -119,16 +139,35 @@ class slvr_dim< blitz::TinyVector zero = blitz::TinyVector({0,0,0}); blitz::thirdIndex vert_idx; + const rng_t &vert_rng = this->k; libmpdataxx::arrvec_t vip_ground; std::set hori_vel = std::set{ix::u, ix::v}; + void GCtoC(arrvec_t &C) + { + using namespace arakawa_c; + const auto &i(this->i), &j(this->j), &k(this->k); + const auto &im(this->im), &jm(this->jm), &km(this->km); + C[0](im + h, j, k) = this->mem->GC[0](im + h, j, k) / (0.5 * ((*this->mem->G)(im, j, k) + (*this->mem->G)(im + 1, j, k))); + C[1](i, jm + h, k) = this->mem->GC[1](i, jm + h, k) / (0.5 * ((*this->mem->G)(i, jm, k) + (*this->mem->G)(i, jm + 1, k))); + C[2](i, j, km + h) = this->mem->GC[2](i, j, km + h) / (0.5 * ((*this->mem->G)(i, j, km) + (*this->mem->G)(i, j, km + 1))); + } + void vert_grad_fwd(typename parent_t::arr_t &in, typename parent_t::arr_t &out, setup::real_t dz) { - in(this->i, this->j, this->k.last() + 1) = in(this->i, this->j, this->k.last()); + // extrapolate upward + in(this->i, this->j, this->k.last() + 1) = 1.5*in(this->i, this->j, this->k.last()) - 0.5*in(this->i, this->j, this->k.last()-1); out(this->i, this->j, this->k) = ( in(this->i, this->j, this->k+1) - in(this->i, this->j, this->k)) / dz; // top and bottom cells are two times lower - out(this->i, this->j, 0) *= 2; out(this->i, this->j, this->k.last()) *= 2; + out(this->i, this->j, 0) *= 2; + +/* + // we don't want surf fluxes to change values at the ground + out(this->i, this->j, 0) = 0; + // instead, the flux from 0-th level goes to the first level above the gfround + out(this->i, this->j, 1) = ( in(this->i, this->j, 2) - in(this->i, this->j, 0)) / (1.5*dz); +*/ } void vert_grad_cnt(typename parent_t::arr_t &in, typename parent_t::arr_t &out, setup::real_t dz) @@ -138,7 +177,8 @@ class slvr_dim< out(this->i, this->j, this->k) = ( in(this->i, this->j, this->k+1) - in(this->i, this->j, this->k-1)) / 2./ dz; // top and bottom cells are two times lower out(this->i, this->j, 0) *= 2; - out(this->i, this->j, this->k.last()) *= 2; + //out(this->i, this->j, this->k.last()) *= 2; + out(this->i, this->j, this->k.last()) = 0; } void smooth(typename parent_t::arr_t &in, typename parent_t::arr_t &out) @@ -154,7 +194,7 @@ class slvr_dim< auto calc_U_ground() return_macro(, - sqrt(pow2(this->state(ix::vip_i)(this->i, this->j, 0).reindex({0,0})) + pow2(this->state(ix::vip_j)(this->i, this->j, 0).reindex({0,0}))) + sqrt(pow2(this->state(ix::vip_i)(this->i, this->j, 0)) + pow2(this->state(ix::vip_j)(this->i, this->j, 0))) ) // ctor @@ -164,8 +204,8 @@ class slvr_dim< ) : parent_t(args, p) { - vip_ground.push_back(new arr_sub_t(this->state(ix::vip_i)(this->i, this->j, 0).reindex({0,0}))); - vip_ground.push_back(new arr_sub_t(this->state(ix::vip_j)(this->i, this->j, 0).reindex({0,0}))); + vip_ground.push_back(new arr_sub_t(this->state(ix::vip_i)(this->i, this->j, 0))); + vip_ground.push_back(new arr_sub_t(this->state(ix::vip_j)(this->i, this->j, 0))); } }; diff --git a/src/slvr_lgrngn.hpp b/src/slvr_lgrngn.hpp index baa1bd1a66..c6b2606e24 100644 --- a/src/slvr_lgrngn.hpp +++ b/src/slvr_lgrngn.hpp @@ -42,6 +42,14 @@ 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 temperature + prtcls->diag_temperature(); + this->record_aux("libcloud_temperature", prtcls->outbuf()); + // recording precipitation rate per grid cel prtcls->diag_all(); prtcls->diag_precip_rate(); @@ -219,6 +227,12 @@ class slvr_lgrngn : public slvr_common { params.flag_coal = params.cloudph_opts.coal; + // use stat_field array to temporarily store 1d pressure profile in a 3d field + // beacuse we need 3d field to init particles + auto& p_e_nd = this->stat_field; + + p_e_nd(this->ijk) = this->p_e(this->vert_idx); + // TODO: barrier? this->mem->barrier(); if (this->rank == 0) @@ -277,15 +291,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); - - prtcls->init( - make_arrinfo(this->mem->advectee(ix::th)), - make_arrinfo(this->mem->advectee(ix::rv)), - make_arrinfo(rhod) - ); + prtcls->init( + make_arrinfo(this->mem->advectee(ix::th)), + make_arrinfo(this->mem->advectee(ix::rv)), + make_arrinfo((*this->mem->G)(this->domain).reindex(this->zero)), + make_arrinfo(p_e_nd(this->domain).reindex(this->zero)) + ); // writing diagnostic data for the initial condition parent_t::tbeg_loop = parent_t::clock::now(); @@ -329,6 +340,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); @@ -361,9 +373,16 @@ class slvr_lgrngn : public slvr_common { parent_t::hook_ante_step(); // includes output this->mem->barrier(); + + // using fluxes array to temporarily store Courant number fields + auto& C = this->flux; + + this->GCtoC(C); + + 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 ( @@ -388,22 +407,14 @@ class slvr_lgrngn : public slvr_common 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(); + auto Cx = C[0](this->Cx_domain).reindex(this->zero); + auto Cy = C[1](this->Cy_domain).reindex(this->zero); + auto Cz = C[this->n_dims - 1](this->Cz_domain).reindex(this->zero); + 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(); @@ -415,27 +426,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; @@ -503,6 +593,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: