From 9ef2d1a237de3204c7e2ce6fb84a84fc75854edd Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Fri, 9 Dec 2022 15:27:05 +0100 Subject: [PATCH 1/8] multiply by rhod for refined arrays --- UWLCM_plotters/include/Plotter3d.hpp | 25 ++++++-- UWLCM_plotters/include/PlotterCommon.hpp | 7 +++ UWLCM_plotters/include/PlotterMicro.hpp | 20 +++++++ drawbicyc/include/plot_series.hpp | 72 ++++++++++++------------ 4 files changed, 83 insertions(+), 41 deletions(-) diff --git a/UWLCM_plotters/include/Plotter3d.hpp b/UWLCM_plotters/include/Plotter3d.hpp index 9d80fed..d481604 100644 --- a/UWLCM_plotters/include/Plotter3d.hpp +++ b/UWLCM_plotters/include/Plotter3d.hpp @@ -18,7 +18,7 @@ class Plotter_t<3> : public PlotterCommon using parent_t = PlotterCommon; hsize_t n[3]; enum {x, y, z}; - arr_t tmp, tmp_srfc; + arr_t tmp, tmp_srfc, tmp_ref; blitz::Range yrange; public: @@ -36,15 +36,17 @@ class Plotter_t<3> : public PlotterCommon cnt[3] = { n[x], n[y], srfc ? 1 : n[z] }, off[3] = { 0, 0, 0 }; this->h5s.selectHyperslab( H5S_SELECT_SET, cnt, off); + + arr_t &arr(tmp.extent(0) == n[x] ? tmp : tmp_ref); // crude check if it is refined or normal data; assume surface data is never refined hsize_t ext[3] = { - hsize_t(tmp.extent(0)), - hsize_t(tmp.extent(1)), - hsize_t(srfc ? tmp_srfc.extent(2) : tmp.extent(2)) + hsize_t(arr.extent(0)), + hsize_t(arr.extent(1)), + hsize_t(srfc ? tmp_srfc.extent(2) : arr.extent(2)) }; - this->h5d.read(srfc ? tmp_srfc.data() : tmp.data(), H5::PredType::NATIVE_FLOAT, H5::DataSpace(3, ext), h5s); + this->h5d.read(srfc ? tmp_srfc.data() : arr.data(), H5::PredType::NATIVE_FLOAT, H5::DataSpace(3, ext), h5s); - return blitz::safeToReturn((srfc ? tmp_srfc : tmp) + 0); + return blitz::safeToReturn((srfc ? tmp_srfc : arr) + 0); } auto h5load_timestep( @@ -235,6 +237,17 @@ class Plotter_t<3> : public PlotterCommon // other dataset are of the size x*z, resize tmp tmp.resize(n[0]-1, n[1]-1, n[2]-1); tmp_srfc.resize(n[0]-1, n[1]-1, 1); + + // init refined data + h5load(file + "/const.h5", "X refined"); + this->map["dx refined"] = tmp(1,0,0) - tmp(0,0,0); + h5load(file + "/const.h5", "Y refined"); + this->map["dy refined"] = tmp(0,1,0) - tmp(0,0,0); + h5load(file + "/const.h5", "Z refined"); + this->map["dz refined"] = tmp(0,0,1) - tmp(0,0,0); + this->CellVol_ref = this->map["dx refined"] * this->map["dy refined"] * this->map["dz refined"]; + this->h5f.openDataSet("X refined").getSpace().getSimpleExtentDims(n, NULL); + tmp_ref.resize(n[0]-1, n[1]-1, n[2]-1); } }; diff --git a/UWLCM_plotters/include/PlotterCommon.hpp b/UWLCM_plotters/include/PlotterCommon.hpp index 1c75c83..29dfd2c 100644 --- a/UWLCM_plotters/include/PlotterCommon.hpp +++ b/UWLCM_plotters/include/PlotterCommon.hpp @@ -93,6 +93,7 @@ class PlotterCommon map["dt"] = h5load_attr(file + "/const.h5", "dt", "advection"); map["outfreq"] = h5load_attr(file + "/const.h5", "outfreq", "user_params"); map["MPI_compiler"] = h5load_attr(file + "/const.h5", "MPI compiler (true/false)", "MPI details"); + map["n_fra_iter"] = h5load_attr(file + "/const.h5", "n_fra_iter", "fractal"); // read number of timesteps hsize_t n; @@ -120,6 +121,12 @@ class PlotterCommon h5s.getSimpleExtentDims(&n, NULL); map_prof.emplace("rhod", arr_prof_t(n)); h5d.read(map_prof["rhod"].data(), H5::PredType::NATIVE_FLOAT); + + // read dry air density profile on refined grid + h5load(file + "/const.h5", "refined rhod"); + h5s.getSimpleExtentDims(&n, NULL); + map_prof.emplace("refined rhod", arr_prof_t(n)); + h5d.read(map_prof["refined rhod"].data(), H5::PredType::NATIVE_FLOAT); } } }; diff --git a/UWLCM_plotters/include/PlotterMicro.hpp b/UWLCM_plotters/include/PlotterMicro.hpp index 1945a10..e920d81 100644 --- a/UWLCM_plotters/include/PlotterMicro.hpp +++ b/UWLCM_plotters/include/PlotterMicro.hpp @@ -20,6 +20,26 @@ class PlotterMicro_t : public Plotter_t const double L_evap = 2264.76e3; // latent heat of evaporation [J/kg] public: + void multiply_by_rhod(arr_t &arr) + { + if(arr.extent(NDims-1) == this->map_prof["rhod"].extent(0)) + arr *= this->map_prof["rhod"](this->LastIndex); + else if(arr.extent(NDims-1) == this->map_prof["rhod refined"].extent(0)) + arr *= this->map_prof["rhod refined"](this->LastIndex); + else + throw std::runtime_error("multiply_by_rhod: input array is neither normal grid size nor refined grid size"); + } + + void multiply_by_CellVol(arr_t &arr) + { + if(arr.extent(NDims-1) == this->map_prof["rhod"].extent(0)) + arr *= this->CellVol; + else if(arr.extent(NDims-1) == this->map_prof["rhod refined"].extent(0)) + arr *= this->CellVol_ref; + else + throw std::runtime_error("multiply_by_CellVol: input array is neither normal grid size nor refined grid size"); + } + // functions for diagnosing fields // // aerosol droplets mixing ratio diff --git a/drawbicyc/include/plot_series.hpp b/drawbicyc/include/plot_series.hpp index 24d3f31..2577dc3 100644 --- a/drawbicyc/include/plot_series.hpp +++ b/drawbicyc/include/plot_series.hpp @@ -129,7 +129,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; //g/kg typename Plotter_t::arr_t snap(tmp); snap += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; //g/kg - snap *= rhod; // water per cubic metre (should be wet density...) + plotter.multiply_by_rhod(snap); 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_series[plt](at) = blitz::mean(plotter.k_i); @@ -224,19 +224,19 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { auto tmp = plotter.h5load_timestep("cloud_rw_mom3", at * n["outfreq"]) * 4./3. * 3.1416 * 1e3; typename Plotter_t::arr_t snap(tmp); - snap *= rhod; + plotter.multiply_by_rhod(snap); res_series[plt](at) += blitz::mean(snap); } { auto tmp = plotter.h5load_timestep("rain_rw_mom3", at * n["outfreq"]) * 4./3. * 3.1416 * 1e3; typename Plotter_t::arr_t snap(tmp); - snap *= rhod; + plotter.multiply_by_rhod(snap); res_series[plt](at) += blitz::mean(snap); } { auto tmp = plotter.h5load_timestep("rv", at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); - snap *= rhod; + plotter.multiply_by_rhod(snap); res_series[plt](at) += blitz::mean(snap); } } @@ -321,7 +321,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) auto tmp2 = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]); typename Plotter_t::arr_t snap_mom(tmp2); com_N_c(at) = snap_mom(com_x_idx(at), com_z_idx(at)); // 0th raw moment / mass [1/kg] - snap_mom *= rhod; // now per m^3 + plotter.multiply_by_rhod(snap); res_series[plt](at) = snap_mom(com_x_idx(at), com_z_idx(at)); } else @@ -436,7 +436,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); snap /= 1e6; // per cm^3 - snap *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap); res_series[plt](at) = blitz::mean(snap); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -449,7 +449,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) auto tmp = plotter.h5load_timestep("all_rw_mom0", at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); snap /= 1e6; // per cm^3 - snap *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap); res_series[plt](at) = blitz::mean(snap); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -462,7 +462,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) auto tmp = plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); snap /= 1e6; // per cm^3 - snap *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap); res_series[plt](at) = blitz::mean(snap); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -475,7 +475,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // cloud fraction (cloudy if N_c > 20/cm^3) auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 typename Plotter_t::arr_t snap2; snap2.resize(snap.shape()); @@ -501,7 +501,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux typename Plotter_t::arr_t snap2(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - snap2 *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap2); snap2 /= 1e6; // per cm^3 snap2 *= snap; @@ -520,12 +520,12 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // cloud fraction (cloudy if N_c > 20/cm^3) auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask auto tmp2 = plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"]); typename Plotter_t::arr_t snap2(tmp2); - snap2 *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap2); snap2 /= 1e6; // per cm^3 snap2 *= snap; if(blitz::sum(snap) > 0) @@ -548,7 +548,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) auto tmp2 = plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"]); typename Plotter_t::arr_t snap2(tmp2); - snap2 *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap2); snap2 /= 1e6; // per cm^3 snap2 *= snap; @@ -567,7 +567,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // cloud fraction (cloudy if N_c > 20/cm^3) auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux @@ -592,7 +592,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // cloud fraction (cloudy if N_c > 20/cm^3) auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux @@ -691,7 +691,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { auto tmp = plotter.h5load_timestep("rd_rng000_mom3", at * n["outfreq"]) * 4./3. * 3.14 * rho_dry * 1e3; typename Plotter_t::arr_t snap(tmp); - snap *= rhod * plotter.CellVol; // turn mixing ratio in g/kg to total mass in g + plotter.multiply_by_rhod(snap); + plotter.multiply_by_CellVol(snap); res_series[plt](at) = blitz::sum(snap); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -724,7 +725,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { // cloud fraction (cloudy if N_c > 20/cm^3) typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask @@ -750,7 +751,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { // cloud fraction (cloudy if N_c > 20/cm^3) typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask @@ -823,7 +824,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; //g/kg typename Plotter_t::arr_t snap(tmp); snap += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; //g/kg - snap *= rhod; // water per cubic metre (should be wet density...) + plotter.multiply_by_rhod(snap); res_series[plt](at) = blitz::mean(snap); } } @@ -836,7 +837,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { { typename Plotter_t::arr_t snap(plotter.h5load_rr_timestep(at * n["outfreq"])); - snap *= rhod * 1e3; // water per cubic metre (should be wet density...) & g/kg + snap *= 1e3; + plotter.multiply_by_rhod(snap); res_series[plt](at) = blitz::mean(snap); } } @@ -1069,11 +1071,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { // cloud fraction (cloudy if N_c > 20/cm^3) typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask typename Plotter_t::arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"])); - snap2 *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap2); snap2 /= 1e6; // per cm^3 snap2 *= snap; if(blitz::sum(snap) > 0) @@ -1090,11 +1092,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { // cloud fraction (cloudy if N_c > 20/cm^3) typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask typename Plotter_t::arr_t snap2(plotter.h5load_timestep("non_gccn_rw_mom0", at * n["outfreq"])); - snap2 *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap2); snap2 /= 1e6; // per cm^3 snap2 *= snap; if(blitz::sum(snap) > 0) @@ -1110,7 +1112,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { // cloud fraction (cloudy if N_c > 20/cm^3) typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask typename Plotter_t::arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"])); @@ -1133,7 +1135,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { // cloud fraction (cloudy if N_c > 20/cm^3) typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask typename Plotter_t::arr_t snap_m0(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"])); @@ -1155,7 +1157,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { typename Plotter_t::arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"])); snap2 /= 1e6; // per cm^3 - snap2 *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap2); res_series[plt](at) = blitz::mean(snap2); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -1167,7 +1169,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { typename Plotter_t::arr_t snap2(plotter.h5load_timestep("non_gccn_rw_mom0", at * n["outfreq"])); snap2 /= 1e6; // per cm^3 - snap2 *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap2); res_series[plt](at) = blitz::mean(snap2); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -1181,11 +1183,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { // cloud fraction (cloudy if N_c > 20/cm^3) typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_geq_0.8um_rw_mom0", at * n["outfreq"])); - snap2 *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap2); snap2 /= 1e6; // per cm^3 snap2 *= snap; if(blitz::sum(snap) > 0) @@ -1202,7 +1204,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_geq_0.8um_rw_mom0", at * n["outfreq"])); snap2 /= 1e6; // per cm^3 - snap2 *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap2); res_series[plt](at) = blitz::mean(snap2); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -1215,11 +1217,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { // cloud fraction (cloudy if N_c > 20/cm^3) typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_lt_0.8um_rw_mom0", at * n["outfreq"])); - snap2 *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap2); snap2 /= 1e6; // per cm^3 snap2 *= snap; if(blitz::sum(snap) > 0) @@ -1235,7 +1237,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_lt_0.8um_rw_mom0", at * n["outfreq"])); snap2 /= 1e6; // per cm^3 - snap2 *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap2); res_series[plt](at) = blitz::mean(snap2); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -1248,7 +1250,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { // cloud fraction (cloudy if N_c > 20/cm^3) typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - snap *= rhod; // b4 it was specific moment + plotter.multiply_by_rhod(snap); snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask typename Plotter_t::arr_t snap_m0(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); From 15ed230456bd6dd5b2e36532db279b5db5843663 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Fri, 9 Dec 2022 16:03:58 +0100 Subject: [PATCH 2/8] rico series fixes --- drawbicyc/include/cases/RICO11/plots.hpp | 4 ++-- drawbicyc/include/gnuplot_series_set_labels.hpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/drawbicyc/include/cases/RICO11/plots.hpp b/drawbicyc/include/cases/RICO11/plots.hpp index e841b4f..72cd397 100644 --- a/drawbicyc/include/cases/RICO11/plots.hpp +++ b/drawbicyc/include/cases/RICO11/plots.hpp @@ -18,8 +18,8 @@ const std::vector series_rico({ "rwp", "surf_precip", "acc_precip" -,"nc" -,"nr" +//,"nc" +//,"nr" ,"cl_nc_rico" ,"cl_nr_rico" ,"cloud_avg_supersat" diff --git a/drawbicyc/include/gnuplot_series_set_labels.hpp b/drawbicyc/include/gnuplot_series_set_labels.hpp index b258c9d..2720752 100644 --- a/drawbicyc/include/gnuplot_series_set_labels.hpp +++ b/drawbicyc/include/gnuplot_series_set_labels.hpp @@ -151,7 +151,7 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt) gp << "set title 'average cloud drop conc [1/cm^3]'\n"; else if (plt == "ntot") gp << "set title 'average particle conc [1/cm^3]'\n"; - else if (plt == "cl_nc") + else if (plt == "cl_nc" || plt == "cl_nc_rico") { gp << "set title 'average cloud drop conc [1/cm^3] in cloudy cells'\n"; gp << "set xlabel ''\n"; @@ -163,7 +163,7 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt) gp << "set xlabel ''\n"; gp << "set ylabel ''\n"; } - else if (plt == "cl_nr") + else if (plt == "cl_nr" || plt == "cl_nr_rico") { gp << "set title 'average rain drop conc [1/cm^3] in cloudy cells'\n"; gp << "set xlabel ''\n"; From 45ed9737450ba8eb225d06b30e61feeb3e2cc913 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Fri, 9 Dec 2022 18:37:08 +0100 Subject: [PATCH 3/8] refined plotting wip' --- UWLCM_plotters/include/Plotter3d.hpp | 13 ++++--- UWLCM_plotters/include/PlotterCommon.hpp | 2 +- UWLCM_plotters/include/PlotterMicro.hpp | 45 +++++++++++++----------- drawbicyc/include/cases/RICO11/plots.hpp | 26 +++++++------- 4 files changed, 46 insertions(+), 40 deletions(-) diff --git a/UWLCM_plotters/include/Plotter3d.hpp b/UWLCM_plotters/include/Plotter3d.hpp index d481604..e2738e5 100644 --- a/UWLCM_plotters/include/Plotter3d.hpp +++ b/UWLCM_plotters/include/Plotter3d.hpp @@ -239,15 +239,14 @@ class Plotter_t<3> : public PlotterCommon tmp_srfc.resize(n[0]-1, n[1]-1, 1); // init refined data - h5load(file + "/const.h5", "X refined"); - this->map["dx refined"] = tmp(1,0,0) - tmp(0,0,0); - h5load(file + "/const.h5", "Y refined"); - this->map["dy refined"] = tmp(0,1,0) - tmp(0,0,0); - h5load(file + "/const.h5", "Z refined"); - this->map["dz refined"] = tmp(0,0,1) - tmp(0,0,0); - this->CellVol_ref = this->map["dx refined"] * this->map["dy refined"] * this->map["dz refined"]; this->h5f.openDataSet("X refined").getSpace().getSimpleExtentDims(n, NULL); tmp_ref.resize(n[0]-1, n[1]-1, n[2]-1); + this->map["refined dx"] = tmp(1,0,0) - tmp(0,0,0); + this->h5f.openDataSet("Y refined").getSpace().getSimpleExtentDims(n, NULL); + this->map["refined dy"] = tmp(0,1,0) - tmp(0,0,0); + this->h5f.openDataSet("Z refined").getSpace().getSimpleExtentDims(n, NULL); + this->map["refined dz"] = tmp(0,0,1) - tmp(0,0,0); + this->CellVol_ref = this->map["refined dx"] * this->map["refined dy"] * this->map["refined dz"]; } }; diff --git a/UWLCM_plotters/include/PlotterCommon.hpp b/UWLCM_plotters/include/PlotterCommon.hpp index 29dfd2c..3f763bb 100644 --- a/UWLCM_plotters/include/PlotterCommon.hpp +++ b/UWLCM_plotters/include/PlotterCommon.hpp @@ -14,7 +14,7 @@ class PlotterCommon std::map map; std::map map_prof; blitz::Array timesteps; - double CellVol, DomainSurf, DomainVol; + double CellVol, DomainSurf, DomainVol, CellVol_ref; protected: H5::H5File h5f; diff --git a/UWLCM_plotters/include/PlotterMicro.hpp b/UWLCM_plotters/include/PlotterMicro.hpp index e920d81..5e60852 100644 --- a/UWLCM_plotters/include/PlotterMicro.hpp +++ b/UWLCM_plotters/include/PlotterMicro.hpp @@ -24,8 +24,8 @@ class PlotterMicro_t : public Plotter_t { if(arr.extent(NDims-1) == this->map_prof["rhod"].extent(0)) arr *= this->map_prof["rhod"](this->LastIndex); - else if(arr.extent(NDims-1) == this->map_prof["rhod refined"].extent(0)) - arr *= this->map_prof["rhod refined"](this->LastIndex); + else if(arr.extent(NDims-1) == this->map_prof["refined rhod"].extent(0)) + arr *= this->map_prof["refined rhod"](this->LastIndex); else throw std::runtime_error("multiply_by_rhod: input array is neither normal grid size nor refined grid size"); } @@ -34,7 +34,7 @@ class PlotterMicro_t : public Plotter_t { if(arr.extent(NDims-1) == this->map_prof["rhod"].extent(0)) arr *= this->CellVol; - else if(arr.extent(NDims-1) == this->map_prof["rhod refined"].extent(0)) + else if(arr.extent(NDims-1) == this->map_prof["refined rhod"].extent(0)) arr *= this->CellVol_ref; else throw std::runtime_error("multiply_by_CellVol: input array is neither normal grid size nor refined grid size"); @@ -45,41 +45,43 @@ class PlotterMicro_t : public Plotter_t // aerosol droplets mixing ratio auto h5load_ra_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) //-> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") - res = this->h5load_timestep("aerosol_rw_mom3", at) * 4./3. * 3.1416 * 1e3; + return arr_t(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); + return res; + // return blitz::safeToReturn(res + 0); + } } // cloud droplets mixing ratio auto h5load_rc_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) //-> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") - res = this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3; + return arr_t(this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3); else if(this->micro == "blk_1m") - res = this->h5load_timestep("rc", at); + return arr_t(this->h5load_timestep("rc", at)); else if(this->micro == "blk_2m") - res = this->h5load_timestep("rc", at); - return blitz::safeToReturn(res + 0); + return arr_t(this->h5load_timestep("rc", at)); } // rain droplets mixing ratio auto h5load_rr_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) //-> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") - res = this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3; + return arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3); else if(this->micro == "blk_1m") - res = this->h5load_timestep("rr", at); + return arr_t(this->h5load_timestep("rr", at)); else if(this->micro == "blk_2m") - res = this->h5load_timestep("rr", at); - return blitz::safeToReturn(res + 0); + return arr_t(this->h5load_timestep("rr", at)); } // activated drops mixing ratio @@ -108,15 +110,18 @@ class PlotterMicro_t : public Plotter_t // cloud droplets concentration [1/kg] auto h5load_nc_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) //-> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") - res = this->h5load_timestep("cloud_rw_mom0", at); + return arr_t(this->h5load_timestep("cloud_rw_mom0", at)); else if(this->micro == "blk_1m") + { res = 0; + return res; + // return blitz::safeToReturn(res + 0); + } else if(this->micro == "blk_2m") - res = this->h5load_timestep("nc", at); - return blitz::safeToReturn(res + 0); + return arr_t(this->h5load_timestep("nc", at); } // precipitation flux [W/m2] diff --git a/drawbicyc/include/cases/RICO11/plots.hpp b/drawbicyc/include/cases/RICO11/plots.hpp index 72cd397..e0828a4 100644 --- a/drawbicyc/include/cases/RICO11/plots.hpp +++ b/drawbicyc/include/cases/RICO11/plots.hpp @@ -1,30 +1,33 @@ #pragma once const std::vector series_rico({ -// "clb_bigrain_mean_rd", -// "clb_bigrain_mean_kappa", -// "clb_bigrain_mean_conc", -// "clb_bigrain_mean_inclt", -//, "clb_bigrain_mean_gccn_fraction" - -//"cl_acnv25_rico", -//"cl_accr25_rico", "RH_max", "cloud_cover_rico", -"min_cloud_base_rico", +"min_cloud_base_rico" "inversion_height_rico", "tot_water", "lwp", "rwp", "surf_precip", "acc_precip" -//,"nc" -//,"nr" ,"cl_nc_rico" ,"cl_nr_rico" ,"cloud_avg_supersat" ,"wvarmax" ,"cl_meanr" //TODO: zmienic maske na rico +,"sd_conc" + + +// "clb_bigrain_mean_rd", +// "clb_bigrain_mean_kappa", +// "clb_bigrain_mean_conc", +// "clb_bigrain_mean_inclt", +//, "clb_bigrain_mean_gccn_fraction" + +//"cl_acnv25_rico", +//"cl_accr25_rico", +//,"nc" +//,"nr" //TODO (po usprawnieniu cloud mask i ujednoliceniu tego: /* ,"cl_avg_cloud_rad" @@ -34,7 +37,6 @@ const std::vector series_rico({ //,"rd_geq_0.8um_conc" //,"cl_rd_lt_0.8um_conc" //,"cl_rd_geq_0.8um_conc" -,"sd_conc" /* "cloud_base", From 274bdca897efb75c8099f593e1e73dd9da687cc2 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 12 Dec 2022 12:07:21 +0100 Subject: [PATCH 4/8] compilation fix --- UWLCM_plotters/include/PlotterMicro.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/UWLCM_plotters/include/PlotterMicro.hpp b/UWLCM_plotters/include/PlotterMicro.hpp index 5e60852..a3a4f4d 100644 --- a/UWLCM_plotters/include/PlotterMicro.hpp +++ b/UWLCM_plotters/include/PlotterMicro.hpp @@ -121,7 +121,7 @@ class PlotterMicro_t : public Plotter_t // return blitz::safeToReturn(res + 0); } else if(this->micro == "blk_2m") - return arr_t(this->h5load_timestep("nc", at); + return arr_t(this->h5load_timestep("nc", at)); } // precipitation flux [W/m2] From bfc92a43ab9fd73ed8a76121012990b36f6d0daf Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 12 Dec 2022 12:20:05 +0100 Subject: [PATCH 5/8] plotter micro return without res cd --- UWLCM_plotters/include/PlotterMicro.hpp | 27 +++++++++++++++---------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/UWLCM_plotters/include/PlotterMicro.hpp b/UWLCM_plotters/include/PlotterMicro.hpp index a3a4f4d..0ec057a 100644 --- a/UWLCM_plotters/include/PlotterMicro.hpp +++ b/UWLCM_plotters/include/PlotterMicro.hpp @@ -87,12 +87,14 @@ class PlotterMicro_t : public Plotter_t // activated drops mixing ratio auto h5load_ract_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) { if(this->micro == "lgrngn") { - res = this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3; - res += arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3); + return arr_t( + arr_t(this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3) + + arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3) + ); } else if(this->micro == "blk_1m") { @@ -104,7 +106,8 @@ class PlotterMicro_t : public Plotter_t res = this->h5load_timestep("rc", at); res += arr_t(this->h5load_timestep("rr", at)); } - return blitz::safeToReturn(res + 0); + // return blitz::safeToReturn(res + 0); + return res; } // cloud droplets concentration [1/kg] @@ -127,14 +130,14 @@ class PlotterMicro_t : public Plotter_t // precipitation flux [W/m2] auto h5load_prflux_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + )// -> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") { - res = this->h5load_timestep("precip_rate", at) + return arr_t(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 - * L_evap; + * L_evap); } else if(this->micro == "blk_1m") try @@ -150,19 +153,21 @@ class PlotterMicro_t : public Plotter_t { res = 0; } - return blitz::safeToReturn(res + 0); + // return blitz::safeToReturn(res + 0); + return res; } // RH auto h5load_RH_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) //-> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") - res = this->h5load_timestep("RH", at); + return arr_t(this->h5load_timestep("RH", at)); else if(this->micro == "blk_1m") res = 0; - return blitz::safeToReturn(res + 0); + // return blitz::safeToReturn(res + 0); + return res; } From a62b28fc14198379ddf75cc990279ffcbcfc4419 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 12 Dec 2022 13:24:01 +0100 Subject: [PATCH 6/8] fixes for refined plotting --- UWLCM_plotters/include/Plotter3d.hpp | 26 ++++++++++++++++++------ UWLCM_plotters/include/PlotterCommon.hpp | 2 +- drawbicyc/include/cases/RICO11/plots.hpp | 16 +++++++-------- drawbicyc/include/plot_series.hpp | 10 ++++----- 4 files changed, 34 insertions(+), 20 deletions(-) diff --git a/UWLCM_plotters/include/Plotter3d.hpp b/UWLCM_plotters/include/Plotter3d.hpp index e2738e5..86692f7 100644 --- a/UWLCM_plotters/include/Plotter3d.hpp +++ b/UWLCM_plotters/include/Plotter3d.hpp @@ -240,13 +240,27 @@ class Plotter_t<3> : public PlotterCommon // init refined data this->h5f.openDataSet("X refined").getSpace().getSimpleExtentDims(n, NULL); - tmp_ref.resize(n[0]-1, n[1]-1, n[2]-1); - this->map["refined dx"] = tmp(1,0,0) - tmp(0,0,0); - this->h5f.openDataSet("Y refined").getSpace().getSimpleExtentDims(n, NULL); - this->map["refined dy"] = tmp(0,1,0) - tmp(0,0,0); - this->h5f.openDataSet("Z refined").getSpace().getSimpleExtentDims(n, NULL); - this->map["refined dz"] = tmp(0,0,1) - tmp(0,0,0); + this->map["refined x"] = n[0]-1; + this->map["refined y"] = n[1]-1; + this->map["refined z"] = n[2]-1; + tmp_ref.resize(n[0], n[1], n[2]); + h5load(file + "/const.h5", "X refined"); + this->map["refined dx"] = tmp_ref(1,0,0) - tmp_ref(0,0,0); + h5load(file + "/const.h5", "Y refined"); + this->map["refined dy"] = tmp_ref(0,1,0) - tmp_ref(0,0,0); + h5load(file + "/const.h5", "Z refined"); + this->map["refined dz"] = tmp_ref(0,0,1) - tmp_ref(0,0,0); this->CellVol_ref = this->map["refined dx"] * this->map["refined dy"] * this->map["refined dz"]; + tmp_ref.resize(n[0]-1, n[1]-1, n[2]-1); + + for (auto const& x : this->map) +{ + std::cout << x.first // string (key) + << ':' + << x.second // string's value + << std::endl; +} + } }; diff --git a/UWLCM_plotters/include/PlotterCommon.hpp b/UWLCM_plotters/include/PlotterCommon.hpp index 3f763bb..0110c5d 100644 --- a/UWLCM_plotters/include/PlotterCommon.hpp +++ b/UWLCM_plotters/include/PlotterCommon.hpp @@ -61,6 +61,7 @@ class PlotterCommon auto attr = h5g.openAttribute(attr_name); notice_macro(std::string("about to read attribute value")) attr.read(attr.getDataType(), &ret); + notice_macro(std::string("attribute value read: ") + std::to_string(ret)) return ret; } @@ -93,7 +94,6 @@ class PlotterCommon map["dt"] = h5load_attr(file + "/const.h5", "dt", "advection"); map["outfreq"] = h5load_attr(file + "/const.h5", "outfreq", "user_params"); map["MPI_compiler"] = h5load_attr(file + "/const.h5", "MPI compiler (true/false)", "MPI details"); - map["n_fra_iter"] = h5load_attr(file + "/const.h5", "n_fra_iter", "fractal"); // read number of timesteps hsize_t n; diff --git a/drawbicyc/include/cases/RICO11/plots.hpp b/drawbicyc/include/cases/RICO11/plots.hpp index e0828a4..0a4ccd3 100644 --- a/drawbicyc/include/cases/RICO11/plots.hpp +++ b/drawbicyc/include/cases/RICO11/plots.hpp @@ -3,19 +3,19 @@ const std::vector series_rico({ "RH_max", "cloud_cover_rico", -"min_cloud_base_rico" +"min_cloud_base_rico", "inversion_height_rico", "tot_water", "lwp", "rwp", "surf_precip", -"acc_precip" -,"cl_nc_rico" -,"cl_nr_rico" -,"cloud_avg_supersat" -,"wvarmax" -,"cl_meanr" //TODO: zmienic maske na rico -,"sd_conc" +"acc_precip", +"cl_nc_rico", +"cl_nr_rico", +"cloud_avg_supersat", +"wvarmax", +"cl_meanr", //TODO: zmienic maske na rico +"sd_conc" // "clb_bigrain_mean_rd", diff --git a/drawbicyc/include/plot_series.hpp b/drawbicyc/include/plot_series.hpp index 2577dc3..054b4f4 100644 --- a/drawbicyc/include/plot_series.hpp +++ b/drawbicyc/include/plot_series.hpp @@ -130,7 +130,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) typename Plotter_t::arr_t snap(tmp); snap += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; //g/kg plotter.multiply_by_rhod(snap); - plotter.k_i = blitz::sum(snap, plotter.LastIndex) * n["dz"]; // LWP [g/m2] in the column + plotter.k_i = blitz::sum(snap, plotter.LastIndex) * n["refined dz"]; // LWP [g/m2] in the column plotter.k_i = where(plotter.k_i > 20 , 1 , 0); // cloudiness as in Ackermann et al. res_series[plt](at) = blitz::mean(plotter.k_i); } @@ -577,7 +577,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) cloudy_column = where(cloudy_column > 0, 1, 0); plotter.k_i = where(cloudy_column == 0, 0, plotter.k_i); if(blitz::sum(cloudy_column) > 0) - res_series[plt](at) = double(blitz::sum(plotter.k_i)) / blitz::sum(cloudy_column) * n["dz"]; + res_series[plt](at) = double(blitz::sum(plotter.k_i)) / blitz::sum(cloudy_column) * n["refined dz"]; else res_series[plt](at) = 0; } @@ -633,7 +633,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) cloudy_column = where(cloudy_column > 0, 1, 0); plotter.k_i = where(cloudy_column == 0, 1e6, plotter.k_i); // 1e6 denotes no clouds in the column if(blitz::sum(cloudy_column) > 0) - res_series[plt](at) = blitz::min(plotter.k_i) * n["dz"]; + res_series[plt](at) = blitz::min(plotter.k_i) * n["refined dz"]; else res_series[plt](at) = 0; } @@ -1828,11 +1828,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) } else if (plt == "lwp") { - res_series[plt] *= (n["z"] - 1) * n["dz"]; // top and bottom cells are smaller + res_series[plt] *= (n["refined z"] - 1) * n["refined dz"]; // top and bottom cells are smaller } else if (plt == "rwp") { - res_series[plt] *= (n["z"] - 1) * n["dz"]; // top and bottom cells are smaller + res_series[plt] *= (n["refined z"] - 1) * n["refined dz"]; // top and bottom cells are smaller } else if (plt == "er") { From 428c88cbddfc9371bb515e5b66adeb779f9b31d6 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Mon, 12 Dec 2022 13:33:58 +0100 Subject: [PATCH 7/8] lwp plot fix --- drawbicyc/include/plot_series.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/drawbicyc/include/plot_series.hpp b/drawbicyc/include/plot_series.hpp index 054b4f4..2e48c77 100644 --- a/drawbicyc/include/plot_series.hpp +++ b/drawbicyc/include/plot_series.hpp @@ -821,11 +821,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { { - auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; //g/kg - typename Plotter_t::arr_t snap(tmp); - snap += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; //g/kg - plotter.multiply_by_rhod(snap); - res_series[plt](at) = blitz::mean(snap); + typename Plotter_t::arr_t rl(plotter.h5load_rc_timestep(at * n["outfreq"])); + typename Plotter_t::arr_t rr(plotter.h5load_rr_timestep(at * n["outfreq"])); + rl = (rl + rr) * 1e3; // g/kg + plotter.multiply_by_rhod(rl); + res_series[plt](at) = blitz::mean(rl); } } catch(...) {if(at==first_timestep) data_found[plt]=0;} From 84c5a6fd7f00ee74492d415b4687bc444f5f05cc Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 14 Dec 2022 15:41:10 +0100 Subject: [PATCH 8/8] handle output from non-refined UWLCM --- UWLCM_plotters/include/Plotter3d.hpp | 40 ++++++++++++++++-------- UWLCM_plotters/include/PlotterCommon.hpp | 16 +++++++--- 2 files changed, 39 insertions(+), 17 deletions(-) diff --git a/UWLCM_plotters/include/Plotter3d.hpp b/UWLCM_plotters/include/Plotter3d.hpp index 86692f7..4f99ed8 100644 --- a/UWLCM_plotters/include/Plotter3d.hpp +++ b/UWLCM_plotters/include/Plotter3d.hpp @@ -239,19 +239,33 @@ class Plotter_t<3> : public PlotterCommon tmp_srfc.resize(n[0]-1, n[1]-1, 1); // init refined data - this->h5f.openDataSet("X refined").getSpace().getSimpleExtentDims(n, NULL); - this->map["refined x"] = n[0]-1; - this->map["refined y"] = n[1]-1; - this->map["refined z"] = n[2]-1; - tmp_ref.resize(n[0], n[1], n[2]); - h5load(file + "/const.h5", "X refined"); - this->map["refined dx"] = tmp_ref(1,0,0) - tmp_ref(0,0,0); - h5load(file + "/const.h5", "Y refined"); - this->map["refined dy"] = tmp_ref(0,1,0) - tmp_ref(0,0,0); - h5load(file + "/const.h5", "Z refined"); - this->map["refined dz"] = tmp_ref(0,0,1) - tmp_ref(0,0,0); - this->CellVol_ref = this->map["refined dx"] * this->map["refined dy"] * this->map["refined dz"]; - tmp_ref.resize(n[0]-1, n[1]-1, n[2]-1); + try + { + this->h5f.openDataSet("X refined").getSpace().getSimpleExtentDims(n, NULL); + this->map["refined x"] = n[0]-1; + this->map["refined y"] = n[1]-1; + this->map["refined z"] = n[2]-1; + tmp_ref.resize(n[0], n[1], n[2]); + h5load(file + "/const.h5", "X refined"); + this->map["refined dx"] = tmp_ref(1,0,0) - tmp_ref(0,0,0); + h5load(file + "/const.h5", "Y refined"); + this->map["refined dy"] = tmp_ref(0,1,0) - tmp_ref(0,0,0); + h5load(file + "/const.h5", "Z refined"); + this->map["refined dz"] = tmp_ref(0,0,1) - tmp_ref(0,0,0); + this->CellVol_ref = this->map["refined dx"] * this->map["refined dy"] * this->map["refined dz"]; + tmp_ref.resize(n[0]-1, n[1]-1, n[2]-1); + } + catch(...) // for pre-refinement simulations that didnt store refined stuff + { + this->map["refined x"] = this->map["x"]; + this->map["refined y"] = this->map["y"]; + this->map["refined z"] = this->map["z"]; + this->map["refined dx"] = this->map["dx"]; + this->map["refined dy"] = this->map["dy"]; + this->map["refined dz"] = this->map["dz"]; + this->CellVol_ref = this->CellVol; + tmp_ref.resize(tmp.shape()); + } for (auto const& x : this->map) { diff --git a/UWLCM_plotters/include/PlotterCommon.hpp b/UWLCM_plotters/include/PlotterCommon.hpp index 0110c5d..5501425 100644 --- a/UWLCM_plotters/include/PlotterCommon.hpp +++ b/UWLCM_plotters/include/PlotterCommon.hpp @@ -123,10 +123,18 @@ class PlotterCommon h5d.read(map_prof["rhod"].data(), H5::PredType::NATIVE_FLOAT); // read dry air density profile on refined grid - h5load(file + "/const.h5", "refined rhod"); - h5s.getSimpleExtentDims(&n, NULL); - map_prof.emplace("refined rhod", arr_prof_t(n)); - h5d.read(map_prof["refined rhod"].data(), H5::PredType::NATIVE_FLOAT); + try + { + h5load(file + "/const.h5", "refined rhod"); + h5s.getSimpleExtentDims(&n, NULL); + map_prof.emplace("refined rhod", arr_prof_t(n)); + h5d.read(map_prof["refined rhod"].data(), H5::PredType::NATIVE_FLOAT); + } + catch(...) // for pre-refinement simulations, use rhod as refined rhod + { + map_prof.emplace("refined rhod", map_prof["rhod"].copy()); + std::cerr << "refined rhod as copy of rhod: " << map_prof["refined rhod"]; + } } } };