Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
addncdata.cpp
Go to the documentation of this file.
1 /*
2  * addncdata.cpp
3  * this tool will take an existing h5m file and add data from an nc type file
4  * will support mainly showing the data associated with unstructured meshes (Homme, MPAS) with Visit
5  *
6  * example of usage:
7  * ./mbaddnc -i wholeFineATM.h5m -n surfdata_ne11np4_simyr1850_c160614.nc -o
8  * whole_LONGXY_surfdata.h5m -v LONGXY
9  *
10  * Basically, will output a new h5m file (whole_LONGXY_surfdata.h5m), which has an extra tag,
11  * corresponding to the variable LONGXY from the file surfdata_ne11np4_simyr1850_c160614.nc;
12  * matching is based on the global ids between what we think is the order on the original file
13  * (wholeFineATM.h5m) and the order of surfdata_ne11np4_simyr1850_c160614.nc
14  *
15  * file wholeFineATM.h5m is obtained from a coupled run in e3sm, with the ne 11, np 4,
16  *
17  * add an option to output the nc data to the original coarse ATM file, the one that also has the
18  * GLOBAL_DOFS tag with the global DOFs of the GLL points
19  */
20 
21 #include "moab/ProgOptions.hpp"
22 #include "moab/Core.hpp"
23 
24 #include "netcdf.h"
25 #include <cmath>
26 #include <sstream>
27 
28 using namespace moab;
29 
30 // copy from ReadNCDF.cpp some useful macros for reading from a netcdf file (exodus?)
31 // ncFile is an integer initialized when opening the nc file in read mode
32 
33 int ncFile;
34 
35 #define GET_DIM( ncdim, name, val ) \
36  { \
37  int gdfail = nc_inq_dimid( ncFile, name, &( ncdim ) ); \
38  if( NC_NOERR == gdfail ) \
39  { \
40  size_t tmp_val; \
41  gdfail = nc_inq_dimlen( ncFile, ncdim, &tmp_val ); \
42  if( NC_NOERR != gdfail ) \
43  { \
44  MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); \
45  } \
46  else \
47  ( val ) = tmp_val; \
48  } \
49  else \
50  ( val ) = 0; \
51  }
52 
53 #define GET_VAR( name, id, dims ) \
54  { \
55  ( id ) = -1; \
56  int gvfail = nc_inq_varid( ncFile, name, &( id ) ); \
57  if( NC_NOERR == gvfail ) \
58  { \
59  int ndims; \
60  gvfail = nc_inq_varndims( ncFile, id, &ndims ); \
61  if( NC_NOERR == gvfail ) \
62  { \
63  ( dims ).resize( ndims ); \
64  gvfail = nc_inq_vardimid( ncFile, id, &( dims )[0] ); \
65  if( NC_NOERR != gvfail ) \
66  { \
67  MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get variable dimension IDs" ); \
68  } \
69  } \
70  } \
71  }
72 
73 #define GET_1D_INT_VAR( name, id, vals ) \
74  { \
75  GET_VAR( name, id, vals ); \
76  if( -1 != ( id ) ) \
77  { \
78  size_t ntmp; \
79  int ivfail = nc_inq_dimlen( ncFile, ( vals )[0], &ntmp ); \
80  if( NC_NOERR != ivfail ) \
81  { \
82  MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); \
83  } \
84  ( vals ).resize( ntmp ); \
85  size_t ntmp1 = 0; \
86  ivfail = nc_get_vara_int( ncFile, id, &ntmp1, &ntmp, &( vals )[0] ); \
87  if( NC_NOERR != ivfail ) \
88  { \
89  MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << ( name ) ); \
90  } \
91  } \
92  }
93 
94 #define GET_1D_DBL_VAR( name, id, vals ) \
95  { \
96  std::vector< int > dum_dims; \
97  GET_VAR( name, id, dum_dims ); \
98  if( -1 != ( id ) ) \
99  { \
100  size_t ntmp; \
101  int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp ); \
102  if( NC_NOERR != dvfail ) \
103  { \
104  MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); \
105  } \
106  ( vals ).resize( ntmp ); \
107  size_t ntmp1 = 0; \
108  dvfail = nc_get_vara_double( ncFile, id, &ntmp1, &ntmp, &( vals )[0] ); \
109  if( NC_NOERR != dvfail ) \
110  { \
111  MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << ( name ) ); \
112  } \
113  } \
114  }
115 
116 #define GET_1D_FLT_VAR( name, id, vals ) \
117  { \
118  std::vector< int > dum_dims; \
119  GET_VAR( name, id, dum_dims ); \
120  if( -1 != ( id ) ) \
121  { \
122  size_t ntmp; \
123  int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp ); \
124  if( NC_NOERR != dvfail ) \
125  { \
126  MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); \
127  } \
128  ( vals ).resize( ntmp ); \
129  size_t ntmp1 = 0; \
130  dvfail = nc_get_vara_float( ncFile, id, &ntmp1, &ntmp, &( vals )[0] ); \
131  if( NC_NOERR != dvfail ) \
132  { \
133  MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << ( name ) ); \
134  } \
135  } \
136  }
137 
138 int main( int argc, char* argv[] )
139 {
140  ErrorCode rval;
141  ProgOptions opts;
142 
143  std::string inputfile, outfile( "out.h5m" ), netcdfFile, variable_name, sefile_name;
144 
145  opts.addOpt< std::string >( "input,i", "input mesh filename", &inputfile );
146  opts.addOpt< std::string >( "netcdfFile,n", "netcdf file aligned with the mesh input file", &netcdfFile );
147  opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );
148 
149  opts.addOpt< std::string >( "var,v", "variable to extract and add to output file", &variable_name );
150 
151  opts.addOpt< std::string >( "sefile,s", "spectral elements file (coarse SE mesh)", &sefile_name );
152  opts.parseCommandLine( argc, argv );
153 
154  Core* mb = new Core();
155 
156  MB_CHK_SET_ERR( mb->load_file( inputfile.c_str() ), "can't load input file" );
157 
158  std::cout << " opened " << inputfile << " with initial h5m data.\n";
159  // open the netcdf file, and see if it has that variable we are looking for
160 
161  Range nodes;
162  MB_CHK_SET_ERR( mb->get_entities_by_dimension( 0, 0, nodes ), "can't get nodes" );
163 
164  Range edges;
165  MB_CHK_SET_ERR( mb->get_entities_by_dimension( 0, 1, edges ), "can't get edges" );
166 
167  Range cells;
168  MB_CHK_SET_ERR( mb->get_entities_by_dimension( 0, 2, cells ), "can't get cells" );
169 
170  std::cout << " it has " << nodes.size() << " vertices " << edges.size() << " edges " << cells.size() << " cells\n";
171 
172  // construct maps between global id and handles
173  std::map< int, EntityHandle > vGidHandle;
174  std::map< int, EntityHandle > eGidHandle;
175  std::map< int, EntityHandle > cGidHandle;
176  std::vector< int > gids;
177  Tag gid;
178  MB_CHK_SET_ERR( mb->tag_get_handle( "GLOBAL_ID", gid ), "can't get global id tag" );
179  gids.resize( nodes.size() );
180  MB_CHK_SET_ERR( mb->tag_get_data( gid, nodes, &gids[0] ), "can't get global id on vertices" );
181  int i = 0;
182  for( Range::iterator vit = nodes.begin(); vit != nodes.end(); vit++ )
183  {
184  vGidHandle[gids[i++]] = *vit;
185  }
186 
187  gids.resize( edges.size() );
188  MB_CHK_SET_ERR( mb->tag_get_data( gid, edges, &gids[0] ), "can't get global id on edges" );
189  i = 0;
190  for( Range::iterator vit = edges.begin(); vit != edges.end(); vit++ )
191  {
192  eGidHandle[gids[i++]] = *vit;
193  }
194 
195  gids.resize( cells.size() );
196  MB_CHK_SET_ERR( mb->tag_get_data( gid, cells, &gids[0] ), "can't get global id on cells" );
197  i = 0;
198  for( Range::iterator vit = cells.begin(); vit != cells.end(); vit++ )
199  {
200  cGidHandle[gids[i++]] = *vit;
201  }
202 
203  // Open netcdf/exodus file
204  int fail = nc_open( netcdfFile.c_str(), 0, &ncFile );
205  if( NC_NOWRITE != fail )
206  {
207  MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf II file " << netcdfFile );
208  }
209 
210  std::cout << " opened " << netcdfFile << " with new data \n";
211  std::vector< int > dims;
212  int nc_var;
213 
214  size_t recs;
215  char recname[NC_MAX_NAME + 1];
216 
217  std::cout << " looking for variable " << variable_name << "\n";
218  GET_VAR( variable_name.c_str(), nc_var, dims );
219  std::cout << " it has " << dims.size() << " dimensions\n";
220 
221  int dimIndex = -1; // index of the dimension of interest
222  bool vertex_data = false;
223  bool edge_data = false;
224  bool cell_data = false;
225  for( size_t j = 0; j < dims.size(); j++ )
226  {
227  fail = nc_inq_dim( ncFile, dims[j], recname, &recs );
228  if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
229  std::string name_dim( recname );
230  std::cout << " dimension index " << j << " in file: " << dims[j] << " name: " << name_dim << " recs:" << recs
231  << "\n";
232  if( recs == nodes.size() )
233  {
234  dimIndex = j;
235  vertex_data = true;
236  }
237  if( recs == edges.size() )
238  {
239  dimIndex = j;
240  edge_data = true;
241  }
242  if( recs == cells.size() )
243  {
244  dimIndex = j;
245  cell_data = true;
246  }
247  }
248 
249  int otherDim = 1 - dimIndex; // used only if 2 dimensions ; could be 0 or 1;
250  size_t size_tag = 1; // good for one dimension
251  std::vector< int > evals; // size of size_tag
252  std::vector< double > dvals; // size of size_tag
253  nc_type dataType;
254 
255  // read the variable, and set it to the tag
256  fail = nc_inq_vartype( ncFile, nc_var, &dataType );
257  if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get variable type" );
258  DataType mbtype = MB_TYPE_DOUBLE;
259  bool float_var = false;
260  if( NC_INT == dataType )
261  mbtype = MB_TYPE_INTEGER;
262  else if( NC_DOUBLE != dataType && NC_FLOAT != dataType )
263  MB_CHK_SET_ERR( MB_FAILURE, "unknown type" );
264 
265  if( NC_FLOAT == dataType ) float_var = true;
266 
267  bool use_time = false;
268  int time_id = -1;
269  std::vector< float > times;
270  fail = nc_inq_varid( ncFile, "time", &time_id );
271  if( NC_NOERR == fail )
272  {
273  use_time = true;
274  int ii;
275  GET_1D_FLT_VAR( "time", ii, times );
276  }
277  Tag newTag;
278 
279  if( ( dims.size() >= 1 && dims.size() <= 2 ) && ( vertex_data || cell_data || edge_data ) )
280  {
281 
282  if( dims.size() == 2 )
283  {
284  fail = nc_inq_dim( ncFile, dims[otherDim], recname, &size_tag );
285  if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
286  }
287 
288  int def_val = 0;
289  MB_CHK_SET_ERR( mb->tag_get_handle( variable_name.c_str(), (int)size_tag, mbtype, newTag,
290  MB_TAG_CREAT | MB_TAG_DENSE, &def_val ),
291  "can't define new tag" );
292 
293  if( NC_INT == dataType )
294  {
295  std::vector< int > vals;
296  if( 1 == dims.size() )
297  {
298  GET_1D_INT_VAR( variable_name.c_str(), dims[0], vals );
299  if( vertex_data )
300  {
301  for( size_t k = 0; k < vals.size(); k++ )
302  {
303  EntityHandle vh = vGidHandle[k + 1];
304  MB_CHK_SET_ERR( mb->tag_set_data( newTag, &vh, 1, &vals[k] ), "can't set tag on vertex" );
305  }
306  }
307  else if( cell_data ) // cell_data
308  {
309  for( size_t k = 0; k < vals.size(); k++ )
310  {
311  EntityHandle ch = cGidHandle[k + 1]; // cell handle
312  MB_CHK_SET_ERR( mb->tag_set_data( newTag, &ch, 1, &vals[k] ), "can't set tag on cell" );
313  }
314  }
315  else if( edge_data ) // edge_data
316  {
317  for( size_t k = 0; k < vals.size(); k++ )
318  {
319  EntityHandle edgeh = eGidHandle[k + 1]; // edge handle
320  rval = mb->tag_set_data( newTag, &edgeh, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on edge" );
321  }
322  }
323  }
324  else // dims.size() == 2
325  {
326  // Single var for all coords
327  size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
328  count[dimIndex] = recs;
329  count[otherDim] = size_tag;
330  vals.resize( recs * size_tag );
331  fail = nc_get_vara_int( ncFile, nc_var, start, count, &vals[0] );
332  evals.resize( size_tag );
333 
334  if( vertex_data )
335  {
336  for( size_t k = 0; k < recs; k++ )
337  {
338  EntityHandle vh = vGidHandle[k + 1];
339  size_t start_in_vals = k * size_tag, stride = 1;
340  for( size_t j = 0; j < size_tag; j++ )
341  evals[j] = vals[start_in_vals + j * stride];
342  MB_CHK_SET_ERR( mb->tag_set_data( newTag, &vh, 1, &evals[0] ), "can't set tag on vertex" );
343  }
344  }
345  else // cell_data
346  {
347  for( size_t k = 0; k < recs; k++ )
348  {
349  EntityHandle ch = cGidHandle[k + 1]; // cell handle
350  size_t start_in_vals = k * size_tag, stride = 1;
351  for( size_t j = 0; j < size_tag; j++ )
352  evals[j] = vals[start_in_vals + j * stride];
353  MB_CHK_SET_ERR( mb->tag_set_data( newTag, &ch, 1, &evals[0] ), "can't set tag on cell" );
354  }
355  }
356  }
357  }
358  else
359  {
360  std::vector< double > vals;
361  if( 1 == dims.size() )
362  {
363  GET_1D_DBL_VAR( variable_name.c_str(), dims[0], vals );
364  if( vertex_data )
365  {
366  for( size_t k = 0; k < vals.size(); k++ )
367  {
368  EntityHandle vh = vGidHandle[k + 1];
369  MB_CHK_SET_ERR( mb->tag_set_data( newTag, &vh, 1, &vals[k] ), "can't set tag on vertex" );
370  }
371  }
372  else if( cell_data ) // cell_data
373  {
374  for( size_t k = 0; k < vals.size(); k++ )
375  {
376  EntityHandle ch = cGidHandle[k + 1]; // cell handle
377  MB_CHK_SET_ERR( mb->tag_set_data( newTag, &ch, 1, &vals[k] ), "can't set tag on cell" );
378  }
379  }
380  else // edge_data
381  {
382  for( size_t k = 0; k < vals.size(); k++ )
383  {
384  EntityHandle edgeh = eGidHandle[k + 1]; // edge handle
385  MB_CHK_SET_ERR( mb->tag_set_data( newTag, &edgeh, 1, &vals[k] ), "can't set tag on edge" );
386  }
387  }
388  }
389  else // dims.size() == 2
390  {
391  // Single var for all coords
392  size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
393  count[dimIndex] = recs;
394  count[otherDim] = size_tag;
395  vals.resize( recs * size_tag );
396  fail = nc_get_vara_double( ncFile, nc_var, start, count, &vals[0] );
397  dvals.resize( size_tag );
398 
399  if( vertex_data )
400  {
401  for( size_t k = 0; k < recs; k++ )
402  {
403  EntityHandle vh = vGidHandle[k + 1];
404  size_t start_in_vals = k * size_tag, stride = 1;
405  for( size_t j = 0; j < size_tag; j++ )
406  dvals[j] = vals[start_in_vals + j * stride];
407  MB_CHK_SET_ERR( mb->tag_set_data( newTag, &vh, 1, &dvals[0] ), "can't set tag on vertex" );
408  }
409  }
410  else // cell_data
411  {
412  for( size_t k = 0; k < recs; k++ )
413  {
414  EntityHandle ch = cGidHandle[k + 1]; // cell handle
415  size_t start_in_vals = k * size_tag, stride = 1;
416  for( size_t j = 0; j < size_tag; j++ )
417  dvals[j] = vals[start_in_vals + j * stride];
418  MB_CHK_SET_ERR( mb->tag_set_data( newTag, &ch, 1, &dvals[0] ), "can't set tag on cell" );
419  }
420  }
421  }
422  }
423  }
424  else if( ( dims.size() == 3 ) && vertex_data && dimIndex == 2 && mbtype == MB_TYPE_DOUBLE &&
425  use_time ) // the last one is the vertex
426  {
427  // the case when the last dim is ncol (for homme type mesh))
428  size_t dim0, dim1; // dim 2 is ncol..
429  char recname0[NC_MAX_NAME + 1], recname1[NC_MAX_NAME + 1];
430  fail = nc_inq_dim( ncFile, dims[0], recname0, &dim0 );
431  if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
432  fail = nc_inq_dim( ncFile, dims[1], recname1, &dim1 );
433  if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
434  std::string name0( recname0 );
435  std::string name1( recname1 );
436  std::string timestr( "time" );
437  size_t start[3] = { 0, 0, 0 };
438  size_t count[3] = { 1, 0, 0 };
439  count[1] = dim1;
440  count[2] = nodes.size();
441  // create a few tags, with name inserted
442  if( name0.compare( "time" ) == 0 )
443  {
444 
445  std::vector< double > dvalues( dim1 * count[2] );
446  std::vector< float > fvalues( dim1 * count[2] );
447  std::vector< double > transp( dim1 * count[2] );
448  // get the variable time
449  for( size_t k = 0; k < dim0; k++ )
450  {
451  // create a tag for each time, and
452  std::stringstream tag_name;
453  tag_name << variable_name << "_t" << times[k];
454  std::vector< double > defvals( dim1, 0. );
455  MB_CHK_SET_ERR( mb->tag_get_handle( tag_name.str().c_str(), (int)dim1, mbtype, newTag,
456  MB_TAG_CREAT | MB_TAG_DENSE, &defvals[0] ),
457  "can't define new tag" );
458  start[0] = k;
459 
460  if( float_var )
461  {
462  fail = nc_get_vara_float( ncFile, nc_var, start, count, &fvalues[0] );
463  // now arrange them in a tag, transpose data
464  for( size_t ii = 0; ii < dim1; ii++ )
465  for( size_t j = 0; j < count[2]; j++ )
466  {
467  transp[j * dim1 + ii] = fvalues[ii * count[2] + j];
468  }
469  }
470  else // double
471  {
472  fail = nc_get_vara_double( ncFile, nc_var, start, count, &dvalues[0] );
473  // now arrange them in a tag, transpose data
474  for( size_t ii = 0; ii < dim1; ii++ )
475  for( size_t j = 0; j < count[2]; j++ )
476  {
477  transp[j * dim1 + ii] = dvalues[ii * count[2] + j];
478  }
479  }
480  for( size_t ii = 0; ii < nodes.size(); ii++ )
481  {
482  EntityHandle vh = vGidHandle[ii + 1];
483  MB_CHK_SET_ERR( mb->tag_set_data( newTag, &vh, 1, &transp[ii * dim1] ), "can't set tag on nodes" );
484  }
485  }
486  }
487  }
488  MB_CHK_SET_ERR( mb->write_file( outfile.c_str() ), "can't write file" );
489  std::cout << " wrote file " << outfile << "\n";
490 
491  // now, if s option, load the coarse mesh and put data on each element, according to a matrix
492  if( !sefile_name.empty() )
493  {
494  // load the file, check for GLOBAL_DOFS tag, and create a new tag with the data associated
495  Core* mb2 = new Core();
496  MB_CHK_SET_ERR( mb2->load_file( sefile_name.c_str() ), "can't load spectral element file" );
497  std::cout << " loaded spectral file " << sefile_name << "\n";
498  // look for GLOBAL_DOFS tag
499  Tag gdofeTag;
500  MB_CHK_SET_ERR( mb2->tag_get_handle( "GLOBAL_DOFS", gdofeTag ), "file does not have GLOBAL_DOFS file" );
501  int sizeTag;
502  MB_CHK_SET_ERR( mb2->tag_get_length( gdofeTag, sizeTag ), "can't get size of tag" );
503  int np = (int)sqrt( 1.0 * sizeTag );
504  std::cout << " size of tag: " << sizeTag << " np = " << np << "\n";
505  std::vector< int > gdofs;
506  Range cells2;
507  MB_CHK_SET_ERR( mb2->get_entities_by_dimension( 0, 2, cells2 ), "can't get cells on spectral mesh" );
508  gdofs.resize( cells2.size() * sizeTag );
509  MB_CHK_SET_ERR( mb2->tag_get_data( gdofeTag, cells2, &gdofs[0] ), "can't get global dofs tag" );
510  // create a new tag for element data arranged as DFIELD
511 
512  std::vector< double > dfield;
513  dfield.resize( sizeTag, 0.0 );
514  Tag newTag2;
515  MB_CHK_SET_ERR( mb2->tag_get_handle( variable_name.c_str(), (int)sizeTag, mbtype, newTag2,
516  MB_TAG_CREAT | MB_TAG_DENSE, &dfield[0] ),
517  "can't define new tag" );
518 
519  int i1 = 0; // index in the gdofs array, per element
520 
521  // get the tag values from the other moab core, for newTag
522  int dataTagLen;
523  MB_CHK_SET_ERR( mb->tag_get_length( newTag, dataTagLen ), "can't get size of newTag" );
524  //
525  std::vector< double > oldData;
526  oldData.resize( dataTagLen * nodes.size() ); //
527 
528  // get the "old" values
529  MB_CHK_SET_ERR( mb->tag_get_data( newTag, nodes, &oldData[0] ), "can't get old values" );
530  for( Range::iterator it = cells2.begin(); it != cells2.end(); it++ )
531  {
532  EntityHandle cel = *it;
533  // gdofs per element are gdofs[i:i+np*np];
534  for( int k = 0; k < sizeTag; k++ )
535  {
536  int gdof = gdofs[i1 + k];
537  EntityHandle node = vGidHandle[gdof];
538  int index = nodes.index( node );
539  // get last layer of data
540  double val = oldData[index * dataTagLen + dataTagLen - 1]; //
541  dfield[k] = val;
542  }
543  i1 = i1 + sizeTag;
544  MB_CHK_SET_ERR( mb2->tag_set_data( newTag2, &cel, 1, &dfield[0] ), "can't set new tag" );
545  }
546 
547  // write the appended file with the new field:
548  MB_CHK_SET_ERR( mb2->write_file( "atm2.h5m" ), "can't write new spectral file" );
549  std::cout << " wrote file atm2.h5m \n";
550  }
551  return 0;
552 }