ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Classes | Typedefs | Enumerations | Functions
ERF_NCWpsFile.H File Reference
#include <sstream>
#include <string>
#include <atomic>
#include <cfenv>
#include "AMReX_FArrayBox.H"
#include "AMReX_IArrayBox.H"
#include "ERF_EpochTime.H"
#include "ERF_NCInterface.H"
Include dependency graph for ERF_NCWpsFile.H:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes

struct  NDArray< DataType >
 

Typedefs

using PlaneVector = amrex::Vector< amrex::FArrayBox >
 

Enumerations

enum class  NC_Data_Dims_Type {
  Time_SL_SN_WE , Time_BT_SN_WE , Time_SN_WE , Time_BT ,
  Time_SL , Time , Time_BdyWidth_BT_SN , Time_BdyWidth_BT_WE ,
  Time_BdyWidth_SN , Time_BdyWidth_WE
}
 

Functions

int BuildFABsFromWRFBdyFile (const std::string &fname, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_xlo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_xhi, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_ylo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_yhi)
 
template<typename DType >
void ReadTimeSliceFromNetCDFFile (const std::string &fname, const int tidx, amrex::Vector< std::string > names, amrex::Vector< NDArray< DType > > &arrays, amrex::Vector< int > &success)
 
template<typename DType >
void ReadNetCDFFile (const std::string &fname, amrex::Vector< std::string > names, amrex::Vector< NDArray< DType > > &arrays, amrex::Vector< int > &success)
 
template<class FAB , typename DType >
void fill_fab_from_arrays (int iv, const amrex::Box &domain, amrex::Vector< NDArray< float >> &nc_arrays, const std::string &var_name, NC_Data_Dims_Type &NC_dim_type, FAB &temp)
 
template<class FAB , typename DType >
void BuildFABsFromNetCDFFile (const amrex::Box &domain, const std::string &fname, amrex::Vector< std::string > nc_var_names, amrex::Vector< enum NC_Data_Dims_Type > NC_dim_types, amrex::Vector< FAB * > fab_vars, amrex::Vector< int > &success)
 

Typedef Documentation

◆ PlaneVector

using PlaneVector = amrex::Vector<amrex::FArrayBox>

Enumeration Type Documentation

◆ NC_Data_Dims_Type

enum NC_Data_Dims_Type
strong
Enumerator
Time_SL_SN_WE 
Time_BT_SN_WE 
Time_SN_WE 
Time_BT 
Time_SL 
Time 
Time_BdyWidth_BT_SN 
Time_BdyWidth_BT_WE 
Time_BdyWidth_SN 
Time_BdyWidth_WE 
33  {
34  Time_SL_SN_WE, // Time, Soil Layers, South-North, West-East
35  Time_BT_SN_WE,
36  Time_SN_WE,
37  Time_BT,
38  Time_SL, // Time, Soil layers
39  Time,
40  Time_BdyWidth_BT_SN,
41  Time_BdyWidth_BT_WE,
42  Time_BdyWidth_SN,
43  Time_BdyWidth_WE
44 };
NC_Data_Dims_Type::Time_BdyWidth_WE
@ Time_BdyWidth_WE
NC_Data_Dims_Type::Time_BdyWidth_BT_SN
@ Time_BdyWidth_BT_SN
NC_Data_Dims_Type::Time_BT
@ Time_BT
NC_Data_Dims_Type::Time_SN_WE
@ Time_SN_WE
NC_Data_Dims_Type::Time_BdyWidth_SN
@ Time_BdyWidth_SN
NC_Data_Dims_Type::Time_BT_SN_WE
@ Time_BT_SN_WE
NC_Data_Dims_Type::Time
@ Time
NC_Data_Dims_Type::Time_BdyWidth_BT_WE
@ Time_BdyWidth_BT_WE
NC_Data_Dims_Type::Time_SL_SN_WE
@ Time_SL_SN_WE
NC_Data_Dims_Type::Time_SL
@ Time_SL

Function Documentation

◆ BuildFABsFromNetCDFFile()

template<class FAB , typename DType >
void BuildFABsFromNetCDFFile ( const amrex::Box &  domain,
const std::string &  fname,
amrex::Vector< std::string >  nc_var_names,
amrex::Vector< enum NC_Data_Dims_Type >  NC_dim_types,
amrex::Vector< FAB * >  fab_vars,
amrex::Vector< int > &  success 
)

Function to read NetCDF variables and fill the corresponding Array4's

Parameters
domainBox specifying the domain covered by the NetCDF fields
fnameName of the NetCDF file to be read
nc_var_namesVariable names in the NetCDF file
NC_dim_typesNetCDF data dimension types
fab_varsFab data we are to fill
successSuccess flags for each requested variable
363 {
364  int ioproc = amrex::ParallelDescriptor::IOProcessorNumber(); // I/O rank
365 
366  amrex::Vector<NDArray<float>> nc_arrays(nc_var_names.size());
367 
368  if (amrex::ParallelDescriptor::IOProcessor())
369  {
370  ReadNetCDFFile(fname, nc_var_names, nc_arrays, success);
371  }
372 
373  amrex::ParallelDescriptor::Bcast(success.dataPtr(), success.size(), ioproc);
374 
375  for (int iv = 0; iv < nc_var_names.size(); iv++)
376  {
377  if (success[iv] == 1) {
378  FAB tmp;
379  if (amrex::ParallelDescriptor::IOProcessor()) {
380  fill_fab_from_arrays<FAB,DType>(iv, domain, nc_arrays, nc_var_names[iv], NC_dim_types[iv], tmp);
381  }
382 
383  int ncomp = tmp.nComp();
384  amrex::Box box = tmp.box();
385 
386  amrex::ParallelDescriptor::Bcast(&box, 1, ioproc);
387  amrex::ParallelDescriptor::Bcast(&ncomp, 1, ioproc);
388 
389  if (!amrex::ParallelDescriptor::IOProcessor()) {
390 #ifdef AMREX_USE_GPU
391  tmp.resize(box,ncomp,amrex::The_Pinned_Arena());
392 #else
393  tmp.resize(box,ncomp);
394 #endif
395  }
396 
397  amrex::ParallelDescriptor::Bcast(tmp.dataPtr(), tmp.size(), ioproc);
398 
399  // Shift box by the domain lower corner
400  amrex::Box fab_bx = tmp.box();
401  amrex::Dim3 dom_lb = lbound(domain);
402  fab_bx += amrex::IntVect(dom_lb.x,dom_lb.y,dom_lb.z);
403  // fab_vars points to data on device
404  fab_vars[iv]->resize(fab_bx,1);
405 #ifdef AMREX_USE_GPU
406  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
407  tmp.dataPtr(), tmp.dataPtr() + tmp.size(),
408  fab_vars[iv]->dataPtr());
409 #else
410  // Provided by BaseFab inheritance through FArrayBox
411  fab_vars[iv]->copy(tmp,tmp.box(),0,fab_bx,0,1);
412 #endif
413  } // success
414  } // iv
415 }
ReadNetCDFFile
void ReadNetCDFFile(const std::string &fname, amrex::Vector< std::string > names, amrex::Vector< NDArray< DType > > &arrays, amrex::Vector< int > &success)
Definition: ERF_NCWpsFile.H:156
WSM6SedCellScratch::tmp
@ tmp
Definition: ERF_AdvanceWSM6.cpp:116
Here is the call graph for this function:

◆ BuildFABsFromWRFBdyFile()

int BuildFABsFromWRFBdyFile ( const std::string &  fname,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  bdy_data_xlo,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  bdy_data_xhi,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  bdy_data_ylo,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  bdy_data_yhi 
)

◆ fill_fab_from_arrays()

template<class FAB , typename DType >
void fill_fab_from_arrays ( int  iv,
const amrex::Box &  domain,
amrex::Vector< NDArray< float >> &  nc_arrays,
const std::string &  var_name,
NC_Data_Dims_Type &  NC_dim_type,
FAB &  temp 
)

Helper function for reading data from NetCDF file into a provided FAB.

Parameters
ivIndex for which variable we are going to fill
domainBox specifying the horizontal and vertical extent to read
nc_arraysArrays of data from NetCDF file
var_nameVariable name
NC_dim_typeDimension type for the variable as stored in the NetCDF file
tempFAB where we store the variable data from the NetCDF Arrays
254 {
255  // bottom_top, south_north, west_east (these can be staggered or unstaggered).
256  // Initialised because amrex::Abort is not [[noreturn]], so the compiler
257  // cannot tell that the final else never falls through.
258  int ns1 = 0, ns2 = 0, ns3 = 0;
259  if (NC_dim_type == NC_Data_Dims_Type::Time_BT) {
260  ns1 = nc_arrays[iv].get_vshape()[1];
261  ns2 = 1;
262  ns3 = 1;
263  // amrex::Print() << "TYPE BT " << ns1 << std::endl;
264  } else if (NC_dim_type == NC_Data_Dims_Type::Time_SN_WE) {
265  ns1 = 1;
266  ns2 = nc_arrays[iv].get_vshape()[1];
267  ns3 = nc_arrays[iv].get_vshape()[2];
268  // amrex::Print() << "TYPE SN WE " << ns2 << " " << ns3 << std::endl;
269  } else if (NC_dim_type == NC_Data_Dims_Type::Time_BT_SN_WE) {
270  ns1 = nc_arrays[iv].get_vshape()[1];
271  ns2 = nc_arrays[iv].get_vshape()[2];
272  ns3 = nc_arrays[iv].get_vshape()[3];
273  // amrex::Print() << "TYPE BT SN WE " << ns1 << " " << ns2 << " " << ns3 << std::endl;
274  } else if (NC_dim_type == NC_Data_Dims_Type::Time_SL_SN_WE) {
275  ns1 = nc_arrays[iv].get_vshape()[1];
276  ns2 = nc_arrays[iv].get_vshape()[2];
277  ns3 = nc_arrays[iv].get_vshape()[3];
278  // amrex::Print() << "TYPE SL SN WE " << ns1 << " " << ns2 << " " << ns3 << std::endl;
279  } else if (NC_dim_type == NC_Data_Dims_Type::Time_SL) {
280  ns1 = nc_arrays[iv].get_vshape()[1];
281  ns2 = 1;
282  ns3 = 1;
283  // amrex::Print() << "TYPE SL " << ns1 << std::endl;
284  }
285  else {
286  amrex::Abort("fill_fab_from_arrays: variable " + var_name +
287  " has NC_Data_Dims_Type " + std::to_string(static_cast<int>(NC_dim_type)) +
288  ", which this reader does not handle (only Time_BT, Time_SN_WE,"
289  " Time_BT_SN_WE, Time_SL_SN_WE and Time_SL)");
290  }
291 
292  // allow fewer z levels to be read in
293  // - ns1 may be the bottom_top or bottom_top_stag dim
294  // - domain.bigEnd(2) is the input number of cells in z
295  // - khi may be nz==bottom_top, nz < bottom_top, or 0 for 2D fields
296  int khi = std::min(domain.bigEnd(2), ns1-1);
297 
298  // TODO: The box will only start at (0,0,0) at level 0 -- we need to generalize this
299  amrex::Box in_box(amrex::IntVect(0,0,0), amrex::IntVect(ns3-1,ns2-1,ns1-1)); // ns1 may or may not be staggered
300  amrex::Box my_box(amrex::IntVect(0,0,0), amrex::IntVect(ns3-1,ns2-1,khi));
301 
302  if (var_name == "PH" || var_name == "PHB") {
303  my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
304  my_box.setBig(2, khi+1);
305  }
306  else if (var_name == "U" || var_name == "UU" || var_name == "MAPFAC_U") {
307  my_box.setType(amrex::IndexType(amrex::IntVect(1,0,0)));
308  }
309  else if (var_name == "V" || var_name == "VV" || var_name == "MAPFAC_V") {
310  my_box.setType(amrex::IndexType(amrex::IntVect(0,1,0)));
311  }
312  else if (var_name == "W" || var_name == "WW") {
313  my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
314  my_box.setBig(2, khi+1);
315  }
316  else if (var_name == "TSLB" || var_name == "SMOIS" || var_name == "SH2O" || var_name == "ZS" || var_name == "DZS") {
317  // soil variables are staggered in z
318  my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
319  }
320 
321  amrex::Arena* Arena_Used = amrex::The_Arena();
322 #ifdef AMREX_USE_GPU
323  // Make sure temp lives on CPU since nc_arrays lives on CPU only
324  Arena_Used = amrex::The_Pinned_Arena();
325 #endif
326  temp.resize(my_box,1, Arena_Used);
327  amrex::Array4<DType> fab_arr = temp.array();
328 
329  int ioff = temp.box().smallEnd()[0];
330  int joff = temp.box().smallEnd()[1];
331  int kmax = temp.box().bigEnd()[2]; // highest z level to be read in
332  //amrex::Print() << var_name << " domain khi, kmax, ns1 = " << khi << " " << kmax << " " << ns1 << std::endl;
333 
334  auto num_pts = in_box.numPts();
335 
336  for (int n(0); n < num_pts; ++n) {
337  int k = n / (ns2*ns3);
338  if (k > kmax) continue;
339  int j = (n - k*(ns2*ns3)) / ns3 + joff;
340  int i = n - k*(ns2*ns3) - (j-joff) * ns3 + ioff;
341  fab_arr(i,j,k,0) = static_cast<DType>(*(nc_arrays[iv].get_data()+n));
342  }
343 }
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
Arena_Used
Arena * Arena_Used
Definition: ERF_Morrison_Advance_F.H:23

◆ ReadNetCDFFile()

template<typename DType >
void ReadNetCDFFile ( const std::string &  fname,
amrex::Vector< std::string >  names,
amrex::Vector< NDArray< DType > > &  arrays,
amrex::Vector< int > &  success 
)
158 {
159  AMREX_ASSERT(arrays.size() == names.size());
160 
161  if (amrex::ParallelDescriptor::IOProcessor())
162  {
163  auto prev_fpe_except = amrex::disableFPExcept(amrex::FPExcept::all);
164 
165  auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4);
166 
167  /*
168  // get the dimension information
169  int Time = static_cast<int>(ncf.dim("Time").len());
170  int DateStrLen = static_cast<int>(ncf.dim("DateStrLen").len());
171  int west_east = static_cast<int>(ncf.dim("west_east").len());
172  int south_north = static_cast<int>(ncf.dim("south_north").len());
173  int bottom_top = static_cast<int>(ncf.dim("bottom_top").len());
174  int bottom_top_stag = static_cast<int>(ncf.dim("bottom_top_stag").len());
175  int west_east_stag = static_cast<int>(ncf.dim("west_east_stag").len());
176  int south_north_stag = static_cast<int>(ncf.dim("south_north_stag").len());
177  int bdy_width = static_cast<int>(ncf.dim("bdy_width").len());
178  */
179 
180  // amrex::Print() << "Reading the dimensions from the netcdf file " << "\n";
181 
182  for (auto n=0; n<arrays.size(); ++n)
183  {
184  std::string vname_to_write = names[n];
185  std::string vname_to_read = names[n];
186 
187  if (vname_to_read.substr(0,2) == "R_") {
188  vname_to_read = names[n+4]; // This allows us to read "T" instead -- we will over-write this later
189  }
190 
191  success[n] = ncf.has_var(vname_to_read);
192  if (success[n] == 0) {
193  amrex::Print() << " Skipping " << vname_to_read << std::endl;
194  } else {
195  amrex::Print() << " Reading " << vname_to_read << std::endl;
196  }
197 
198  if (success[n] == 1) {
199 
200  std::vector<std::string> dimnames = ncf.var(vname_to_read).dimnames();
201  AMREX_ALWAYS_ASSERT(dimnames[0] == "Time" || dimnames[0] == "time");
202 
203  std::vector<size_t> shape = ncf.var(vname_to_read).shape();
204  arrays[n] = NDArray<DType>(vname_to_read,shape);
205  DType* dataPtr = arrays[n].get_data();
206 
207  std::vector<size_t> start(shape.size(), 0);
208 
209 #if 0
210  auto numPts = arrays[n].ndim();
211  amrex::Print() << "NetCDF Variable name = " << vname_to_read << std::endl;
212  amrex::Print() << "numPts read from NetCDF file/var = " << numPts << std::endl;
213  amrex::Print() << "Dims in var = " << ncf.var(vname_to_read).ndim() << std::endl;
214  amrex::Print() << "Dim names = (";
215  for (auto &dim:dimnames)
216  amrex::Print() << dim << ", " ;
217  amrex::Print() << ")" << std::endl;
218  amrex::Print() << "Dims of the variable = (";
219  for (auto &dim:shape)
220  amrex::Print() << dim << ", " ;
221  amrex::Print() << ")" << std::endl;
222 #endif
223 
224  ncf.var(vname_to_read).get(dataPtr, start, shape);
225 
226  } // has_var
227  }
228  ncf.close();
229 
230  std::feclearexcept(FE_ALL_EXCEPT);
231  amrex::setFPExcept(prev_fpe_except);
232  }
233 }
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
NDArray
Definition: ERF_NCWpsFile.H:56

Referenced by BuildFABsFromNetCDFFile().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ReadTimeSliceFromNetCDFFile()

template<typename DType >
void ReadTimeSliceFromNetCDFFile ( const std::string &  fname,
const int  tidx,
amrex::Vector< std::string >  names,
amrex::Vector< NDArray< DType > > &  arrays,
amrex::Vector< int > &  success 
)
114 {
115  amrex::Print() << "Reading time slice " << tidx << " from " << fname << std::endl;
116  AMREX_ASSERT(arrays.size() == names.size());
117 
118  if (amrex::ParallelDescriptor::IOProcessor())
119  {
120  auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4);
121 
122  int ntimes = ncf.dim("Time").len();
123  AMREX_ALWAYS_ASSERT((tidx >= 0) && (tidx < ntimes));
124 
125  for (auto n=0; n<arrays.size(); ++n)
126  {
127  std::string vname_to_write = names[n];
128  std::string vname_to_read = names[n];
129  if (vname_to_read.substr(0,2) == "R_") {
130  vname_to_read = names[n+4]; // This allows us to read "T" instead -- we will over-write this later
131  }
132 
133  success[n] = ncf.has_var(vname_to_read);
134 
135  if (success[n] == 1)
136  {
137  std::vector<std::string> dimnames = ncf.var(vname_to_read).dimnames();
138  AMREX_ALWAYS_ASSERT(dimnames[0] == "Time");
139 
140  std::vector<size_t> count = ncf.var(vname_to_read).shape();
141  std::vector<size_t> start(count.size(), 0);
142  start[0] = tidx;
143  count[0] = 1;
144 
145  arrays[n] = NDArray<DType>(vname_to_read, count);
146  DType* dataPtr = arrays[n].get_data();
147 
148  ncf.var(vname_to_read).get(dataPtr, start, count);
149  } // has_var
150  }
151  ncf.close();
152  }
153 }
Here is the call graph for this function: