ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_NCPlotFile.cpp File Reference
#include <iomanip>
#include <iostream>
#include <string>
#include <ctime>
#include <AMReX_Utility.H>
#include <AMReX_MultiFab.H>
#include "ERF_NumericalConstants.H"
#include "ERF_DataStruct.H"
#include "ERF_NCPlotFile.H"
#include "ERF_NCInterface.H"
Include dependency graph for ERF_NCPlotFile.cpp:

Functions

void writeNCPlotFile (int lev, int which_subdomain, const std::string &dir, const Vector< const MultiFab * > &plotMF, const Vector< std::string > &plot_var_names, const Vector< int > &, Array< Real, AMREX_SPACEDIM > prob_lo, Array< Real, AMREX_SPACEDIM > prob_hi, Array< Real, AMREX_SPACEDIM > dx_in, const Box &subdomain, const double &time, const double &start_bdy_time, const SolverChoice &solverChoice, const Vector< Real > &zlevels_stag, const NCProvenance &provenance)
 

Function Documentation

◆ writeNCPlotFile()

void writeNCPlotFile ( int  lev,
int  which_subdomain,
const std::string &  dir,
const Vector< const MultiFab * > &  plotMF,
const Vector< std::string > &  plot_var_names,
const Vector< int > &  ,
Array< Real, AMREX_SPACEDIM >  prob_lo,
Array< Real, AMREX_SPACEDIM >  prob_hi,
Array< Real, AMREX_SPACEDIM >  dx_in,
const Box &  subdomain,
const double &  time,
const double &  start_bdy_time,
const SolverChoice &  solverChoice,
const Vector< Real > &  zlevels_stag,
const NCProvenance &  provenance 
)
33 {
34  //
35  // Set the full IO path for NetCDF output
36  //
37  std::string FullPath = dir;
38  if (lev == 0) {
39  const std::string& extension = amrex::Concatenate("_d",lev+1,2);
40  FullPath += extension + ".nc";
41  } else {
42  const std::string& extension = amrex::Concatenate("_d",lev+1+which_subdomain,2);
43  FullPath += extension + ".nc";
44  }
45 
46  Print() << "Writing level " << lev << " NetCDF plot file " << FullPath << std::endl;
47 
48  //
49  // Open netcdf file to write data
50  //
51  auto ncf = ncutils::NCFile::create_par(FullPath, NC_NETCDF4 | NC_MPIIO,
52  amrex::ParallelContext::CommunicatorSub(), MPI_INFO_NULL);
53 
54  auto ba = plotMF[lev]->boxArray();
55  auto dm = plotMF[lev]->DistributionMap();
56 
57  int nblocks = ba.size();
58 
59  int nx = subdomain.length(0);
60  int ny = subdomain.length(1);
61  int nz = subdomain.length(2);
62 
63  int num_pts = nx*ny*nz;
64 
65  int n_data_items = plotMF[lev]->nComp();
66 
67  //
68  // Start Define stuff
69  //
70  ncf.enter_def_mode();
71 
72  const std::string nt_name = "num_time_steps";
73  const std::string nb_name = "num_blocks";
74  const std::string np_name = "num_pts";
75  const std::string nx_name = "nx";
76  const std::string ny_name = "ny";
77  const std::string nz_name = "nz";
78 
79  const std::string ndim_name = "num_geo_dims";
80 
81  ncf.put_attr("title", "ERF NetCDF Plot data output");
82 
83  // Provenance of the run that produced the data. Each is written only when
84  // the caller supplied it, so a build without Noah-MP omits that attribute
85  // and the PlotfileToNetCDF converter writes none of the three.
86  if (!provenance.erf_git_hash.empty()) {
87  ncf.put_attr("erf_git_hash", provenance.erf_git_hash);
88  }
89  if (!provenance.amrex_git_hash.empty()) {
90  ncf.put_attr("amrex_git_hash", provenance.amrex_git_hash);
91  }
92  if (!provenance.noahmp_git_hash.empty()) {
93  ncf.put_attr("noahmp_git_hash", provenance.noahmp_git_hash);
94  }
95 
96  ncf.def_dim(ndim_name, AMREX_SPACEDIM);
97  ncf.def_dim(np_name , num_pts);
98  ncf.def_dim(nb_name , nblocks);
99 
100  ncf.def_dim(nt_name, NC_UNLIMITED);
101  ncf.def_dim(nx_name, nx);
102  ncf.def_dim(ny_name, ny);
103  ncf.def_dim(nz_name, nz);
104 
105  ncf.def_var("probLo" , NC_FLOAT, {ndim_name});
106  ncf.def_var("probHi" , NC_FLOAT, {ndim_name});
107 
108  ncf.def_var("Geom.smallend", NC_INT, {ndim_name});
109  ncf.def_var("Geom.bigend" , NC_INT, {ndim_name});
110  ncf.def_var("CellSize" , NC_FLOAT, {ndim_name});
111 
112  ncf.def_var("x_grid", NC_DOUBLE, {np_name});
113  ncf.def_var("y_grid", NC_DOUBLE, {np_name});
114  ncf.def_var("z_grid", NC_DOUBLE, {np_name});
115 
116  for (int i = 0; i < plot_var_names.size(); i++) {
117  ncf.def_var(plot_var_names[i], NC_DOUBLE, {nz_name, ny_name, nx_name});
118  }
119 
120  ncf.exit_def_mode();
121  //
122  // End Define stuff
123  //
124 
125  // We are doing single-level writes but it doesn't have to be level 0
126  //
127  // Write out the netcdf plotfile head information.
128  //
129  if (n_data_items == 0) {
130  amrex::Error("Must specify at least one valid data item to plot");
131  }
132 
133  ncf.put_attr("number_variables", std::vector<int>{n_data_items});
134  ncf.put_attr("space_dimension", std::vector<int>{AMREX_SPACEDIM});
135  ncf.put_attr("current_time", std::vector<double>{time});
136  ncf.put_attr("start_time", std::vector<double>{start_bdy_time});
137  ncf.put_attr("CurrentLevel", std::vector<int>{lev});
138 
139  Real dx[AMREX_SPACEDIM];
140  for (int i = 0; i < AMREX_SPACEDIM; i++) {
141  dx[i] = dx_in[i];
142  }
143 
144  amrex::Vector<Real> probLo;
145  amrex::Vector<Real> probHi;
146  for (int i = 0; i < AMREX_SPACEDIM; i++) {
147  probLo.push_back(prob_lo[i]);
148  probHi.push_back(prob_hi[i]);
149  }
150 
151  auto nc_probLo = ncf.var("probLo");
152  nc_probLo.par_access(NC_COLLECTIVE);
153  nc_probLo.put(probLo.data(), {0}, {AMREX_SPACEDIM});
154 
155  auto nc_probHi = ncf.var("probHi");
156  nc_probHi.par_access(NC_COLLECTIVE);
157  nc_probHi.put(probHi.data(), {0}, {AMREX_SPACEDIM});
158 
159  amrex::Vector<int> smallend;
160  amrex::Vector<int> bigend;
161  smallend.clear(); bigend.clear();
162  for (int j = 0; j < AMREX_SPACEDIM; j++) {
163  smallend.push_back(subdomain.smallEnd(j));
164  bigend.push_back(subdomain.bigEnd(j));
165  }
166 
167  auto nc_Geom_smallend = ncf.var("Geom.smallend");
168  nc_Geom_smallend.par_access(NC_COLLECTIVE);
169  nc_Geom_smallend.put(smallend.data(), {0}, {AMREX_SPACEDIM});
170 
171  auto nc_Geom_bigend = ncf.var("Geom.bigend");
172  nc_Geom_bigend.par_access(NC_COLLECTIVE);
173  nc_Geom_bigend.put(bigend.data(), {0}, {AMREX_SPACEDIM});
174 
175  amrex::Vector<Real> CellSize;
176  CellSize.clear();
177  for (Real& j : dx) {
178  CellSize.push_back(j);
179  }
180  auto nc_CellSize = ncf.var("CellSize");
181  nc_CellSize.par_access(NC_COLLECTIVE);
182  nc_CellSize.put(CellSize.data(), {0}, {AMREX_SPACEDIM});
183 
184  ncf.put_attr("DefaultGeometry", std::vector<int>{amrex::DefaultGeometry().Coord()});
185 
186  std::vector<Real> x_grid;
187  std::vector<Real> y_grid;
188  std::vector<Real> z_grid;
189  long unsigned goffset = 0;
190  long unsigned glen = 0;
191 
192  // *******************************************************************************
193  // NOTE: the (x,y,z) output here are for a mesh withOUT terrain-fitted coordinates
194  // *******************************************************************************
195  if (solverChoice.mesh_type == MeshType::ConstantDz) {
196  for (int i = 0; i < ba.size(); ++i) {
197  auto bx = ba[i];
198  if (subdomain.contains(bx)) {
199  // The loop indices below are box-local, so they must be offset by the
200  // low corner of the box to give the global location of each cell
201  const auto bx_lo = bx.smallEnd();
202  x_grid.clear(); y_grid.clear(); z_grid.clear();
203  for (auto k3 = 0; k3 < bx.length(2); ++k3) {
204  for (auto k2 = 0; k2 < bx.length(1); ++k2) {
205  for (auto k1 = 0; k1 < bx.length(0); ++k1) {
206  x_grid.push_back(prob_lo[0]+dx[0]*(static_cast<Real>(bx_lo[0]+k1)+myhalf));
207  y_grid.push_back(prob_lo[1]+dx[1]*(static_cast<Real>(bx_lo[1]+k2)+myhalf));
208  z_grid.push_back(prob_lo[2]+dx[2]*(static_cast<Real>(bx_lo[2]+k3)+myhalf));
209  }
210  }
211  }
212 
213  goffset += glen;
214  glen = bx.numPts();
215 
216  auto nc_x_grid = ncf.var("x_grid");
217  auto nc_y_grid = ncf.var("y_grid");
218  auto nc_z_grid = ncf.var("z_grid");
219 
220  nc_x_grid.par_access(NC_COLLECTIVE);
221  nc_y_grid.par_access(NC_COLLECTIVE);
222  nc_z_grid.par_access(NC_COLLECTIVE);
223 
224  nc_x_grid.put(x_grid.data(), {goffset}, {glen});
225  nc_y_grid.put(y_grid.data(), {goffset}, {glen});
226  nc_z_grid.put(z_grid.data(), {goffset}, {glen});
227  }
228  }
229  }
230  else if (solverChoice.mesh_type == MeshType::StretchedDz)
231  {
232  for (int i = 0; i < ba.size(); ++i) {
233  auto bx = ba[i];
234  if (subdomain.contains(bx)) {
235  // The loop indices below are box-local, so they must be offset by the
236  // low corner of the box to give the global location of each cell;
237  // zlevels_stag is indexed globally as well
238  const auto bx_lo = bx.smallEnd();
239  x_grid.clear(); y_grid.clear(); z_grid.clear();
240  for (auto k3 = 0; k3 < bx.length(2); ++k3) {
241  for (auto k2 = 0; k2 < bx.length(1); ++k2) {
242  for (auto k1 = 0; k1 < bx.length(0); ++k1) {
243  x_grid.push_back(prob_lo[0]+dx[0]*(static_cast<Real>(bx_lo[0]+k1)+myhalf));
244  y_grid.push_back(prob_lo[1]+dx[1]*(static_cast<Real>(bx_lo[1]+k2)+myhalf));
245  z_grid.push_back(myhalf * (zlevels_stag[bx_lo[2]+k3] + zlevels_stag[bx_lo[2]+k3+1]));
246  }
247  }
248  }
249 
250  goffset += glen;
251  glen = bx.numPts();
252 
253  auto nc_x_grid = ncf.var("x_grid");
254  auto nc_y_grid = ncf.var("y_grid");
255  auto nc_z_grid = ncf.var("z_grid");
256 
257  nc_x_grid.par_access(NC_COLLECTIVE);
258  nc_y_grid.par_access(NC_COLLECTIVE);
259  nc_z_grid.par_access(NC_COLLECTIVE);
260 
261  nc_x_grid.put(x_grid.data(), {goffset}, {glen});
262  nc_y_grid.put(y_grid.data(), {goffset}, {glen});
263  nc_z_grid.put(z_grid.data(), {goffset}, {glen});
264  }
265  }
266  }
267 
268  const int ncomp = plotMF[lev]->nComp();
269 
270  // *******************************************************************************
271  // The data below is written from inside an MFIter loop, so each rank makes one
272  // put call per box that it owns. Ranks can own different numbers of boxes when
273  // the domain does not decompose evenly, so the sequence of put calls does not
274  // match across ranks -- we must therefore use independent rather than collective
275  // access, otherwise the write hangs waiting for ranks that have no box left.
276  //
277  // Note also that the access mode is a property of the variable, so it only needs
278  // to be set once per variable rather than once per box per component.
279  // *******************************************************************************
280  for (int k(0); k < ncomp; ++k) {
281  ncf.var(plot_var_names[k]).par_access(NC_INDEPENDENT);
282  }
283 
284  for (MFIter mfi(*plotMF[lev]); mfi.isValid(); ++mfi)
285  {
286  auto bx = mfi.validbox();
287 
288  if (subdomain.contains(bx))
289 
290  {
291  //
292  // These are the dimensions of the data we write for only this box
293  //
294  long unsigned local_nx = bx.length()[0];
295  long unsigned local_ny = bx.length()[1];
296  long unsigned local_nz = bx.length()[2];
297 
298  long unsigned local_start_x = static_cast<long unsigned>(bx.smallEnd()[0]-subdomain.smallEnd()[0]);
299  long unsigned local_start_y = static_cast<long unsigned>(bx.smallEnd()[1]-subdomain.smallEnd()[1]);
300  long unsigned local_start_z = static_cast<long unsigned>(bx.smallEnd()[2]-subdomain.smallEnd()[2]);
301 
302  for (int k(0); k < ncomp; ++k) {
303  FArrayBox tmp;
304  tmp.resize(bx, 1, amrex::The_Pinned_Arena());
305  tmp.template copy<RunOn::Device>((*plotMF[lev])[mfi.index()], k, 0, 1);
306  Gpu::streamSynchronize();
307 
308  auto nc_plot_var = ncf.var(plot_var_names[k]);
309  nc_plot_var.put(tmp.dataPtr(), {local_start_z,local_start_y,local_start_x},
310  {local_nz, local_ny, local_nx});
311  }
312  }
313  }
314  ncf.close();
315 }
nx
const int nx
Definition: ERF_InitCustomPertVels_CloudChamber.H:14
ny
const int ny
Definition: ERF_InitCustomPertVels_CloudChamber.H:15
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
prob_lo
const GpuArray< Real, AMREX_SPACEDIM > prob_lo
Definition: ERF_InitCustomPert_CloudChamber.H:33
prob_hi
const amrex::Real * prob_hi
Definition: ERF_InitCustomPert_DataAssimilation_ISV.H:17
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
WSM6SedCellScratch::tmp
@ tmp
Definition: ERF_AdvanceWSM6.cpp:116
module_mp_morr_two_moment::k1
real(c_double), private k1
Definition: ERF_module_mp_morr_two_moment.F90:213
SolverChoice::mesh_type
static MeshType mesh_type
Vertical mesh representation.
Definition: ERF_DataStruct.H:2236

Referenced by main_main(), ERF::Write2DPlotFile(), and ERF::Write3DPlotFile().

Here is the caller graph for this function: