Mesh Oriented datABase  (version 5.5.1)
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  rval = mb->load_file( inputfile.c_str() );MB_CHK_SET_ERR( rval, "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  rval = mb->get_entities_by_dimension( 0, 0, nodes );MB_CHK_SET_ERR( rval, "can't get nodes" );
163 
164  Range edges;
165  rval = mb->get_entities_by_dimension( 0, 1, edges );MB_CHK_SET_ERR( rval, "can't get edges" );
166 
167  Range cells;
168  rval = mb->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "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  rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_SET_ERR( rval, "can't get global id tag" );
179  gids.resize( nodes.size() );
180  rval = mb->tag_get_data( gid, nodes, &gids[0] );MB_CHK_SET_ERR( rval, "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  rval = mb->tag_get_data( gid, edges, &gids[0] );MB_CHK_SET_ERR( rval, "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  rval = mb->tag_get_data( gid, cells, &gids[0] );MB_CHK_SET_ERR( rval, "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 cell_data = false;
224  for( size_t j = 0; j < dims.size(); j++ )
225  {
226  fail = nc_inq_dim( ncFile, dims[j], recname, &recs );
227  if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
228  std::string name_dim( recname );
229  std::cout << " dimension index " << j << " in file: " << dims[j] << " name: " << name_dim << " recs:" << recs
230  << "\n";
231  if( recs == nodes.size() )
232  {
233  dimIndex = j;
234  vertex_data = true;
235  }
236  if( recs == cells.size() )
237  {
238  dimIndex = j;
239  cell_data = true;
240  }
241  }
242 
243  int otherDim = 1 - dimIndex; // used only if 2 dimensions ; could be 0 or 1;
244  size_t size_tag = 1; // good for one dimension
245  std::vector< int > evals; // size of size_tag
246  std::vector< double > dvals; // size of size_tag
247  nc_type dataType;
248 
249  // read the variable, and set it to the tag
250  fail = nc_inq_vartype( ncFile, nc_var, &dataType );
251  if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get variable type" );
252  DataType mbtype = MB_TYPE_DOUBLE;
253  bool float_var = false;
254  if( NC_INT == dataType )
255  mbtype = MB_TYPE_INTEGER;
256  else if( NC_DOUBLE != dataType && NC_FLOAT != dataType )
257  MB_CHK_SET_ERR( MB_FAILURE, "unknown type" );
258 
259  if( NC_FLOAT == dataType ) float_var = true;
260 
261  bool use_time = false;
262  int time_id = -1;
263  std::vector< float > times;
264  fail = nc_inq_varid( ncFile, "time", &time_id );
265  if( NC_NOERR == fail )
266  {
267  use_time = true;
268  int ii;
269  GET_1D_FLT_VAR( "time", ii, times );
270  }
271  Tag newTag;
272 
273  if( ( dims.size() >= 1 && dims.size() <= 2 ) && ( vertex_data || cell_data ) )
274  {
275 
276  if( dims.size() == 2 )
277  {
278  fail = nc_inq_dim( ncFile, dims[otherDim], recname, &size_tag );
279  if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
280  }
281 
282  int def_val = 0;
283  rval = mb->tag_get_handle( variable_name.c_str(), (int)size_tag, mbtype, newTag, MB_TAG_CREAT | MB_TAG_DENSE,
284  &def_val );MB_CHK_SET_ERR( rval, "can't define new tag" );
285 
286  if( NC_INT == dataType )
287  {
288  std::vector< int > vals;
289  if( 1 == dims.size() )
290  {
291  GET_1D_INT_VAR( variable_name.c_str(), dims[0], vals );
292  if( vertex_data )
293  {
294  for( size_t k = 0; k < vals.size(); k++ )
295  {
296  EntityHandle vh = vGidHandle[k + 1];
297  rval = mb->tag_set_data( newTag, &vh, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
298  }
299  }
300  else // cell_data
301  {
302  for( size_t k = 0; k < vals.size(); k++ )
303  {
304  EntityHandle ch = cGidHandle[k + 1]; // cell handle
305  rval = mb->tag_set_data( newTag, &ch, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on cell" );
306  }
307  }
308  }
309  else // dims.size() == 2
310  {
311  // Single var for all coords
312  size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
313  count[dimIndex] = recs;
314  count[otherDim] = size_tag;
315  vals.resize( recs * size_tag );
316  fail = nc_get_vara_int( ncFile, nc_var, start, count, &vals[0] );
317  evals.resize( size_tag );
318 
319  if( vertex_data )
320  {
321  for( size_t k = 0; k < recs; k++ )
322  {
323  EntityHandle vh = vGidHandle[k + 1];
324  size_t start_in_vals = k * size_tag, stride = 1;
325  for( size_t j = 0; j < size_tag; j++ )
326  evals[j] = vals[start_in_vals + j * stride];
327  rval = mb->tag_set_data( newTag, &vh, 1, &evals[0] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
328  }
329  }
330  else // cell_data
331  {
332  for( size_t k = 0; k < recs; k++ )
333  {
334  EntityHandle ch = cGidHandle[k + 1]; // cell handle
335  size_t start_in_vals = k * size_tag, stride = 1;
336  for( size_t j = 0; j < size_tag; j++ )
337  evals[j] = vals[start_in_vals + j * stride];
338  rval = mb->tag_set_data( newTag, &ch, 1, &evals[0] );MB_CHK_SET_ERR( rval, "can't set tag on cell" );
339  }
340  }
341  }
342  }
343  else
344  {
345  std::vector< double > vals;
346  if( 1 == dims.size() )
347  {
348  GET_1D_DBL_VAR( variable_name.c_str(), dims[0], vals );
349  if( vertex_data )
350  {
351  for( size_t k = 0; k < vals.size(); k++ )
352  {
353  EntityHandle vh = vGidHandle[k + 1];
354  rval = mb->tag_set_data( newTag, &vh, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
355  }
356  }
357  else // cell_data
358  {
359  for( size_t k = 0; k < vals.size(); k++ )
360  {
361  EntityHandle ch = cGidHandle[k + 1]; // cell handle
362  rval = mb->tag_set_data( newTag, &ch, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
363  }
364  }
365  }
366  else // dims.size() == 2
367  {
368  // Single var for all coords
369  size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
370  count[dimIndex] = recs;
371  count[otherDim] = size_tag;
372  vals.resize( recs * size_tag );
373  fail = nc_get_vara_double( ncFile, nc_var, start, count, &vals[0] );
374  dvals.resize( size_tag );
375 
376  if( vertex_data )
377  {
378  for( size_t k = 0; k < recs; k++ )
379  {
380  EntityHandle vh = vGidHandle[k + 1];
381  size_t start_in_vals = k * size_tag, stride = 1;
382  for( size_t j = 0; j < size_tag; j++ )
383  dvals[j] = vals[start_in_vals + j * stride];
384  rval = mb->tag_set_data( newTag, &vh, 1, &dvals[0] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
385  }
386  }
387  else // cell_data
388  {
389  for( size_t k = 0; k < recs; k++ )
390  {
391  EntityHandle ch = cGidHandle[k + 1]; // cell handle
392  size_t start_in_vals = k * size_tag, stride = 1;
393  for( size_t j = 0; j < size_tag; j++ )
394  dvals[j] = vals[start_in_vals + j * stride];
395  rval = mb->tag_set_data( newTag, &ch, 1, &dvals[0] );MB_CHK_SET_ERR( rval, "can't set tag on cell" );
396  }
397  }
398  }
399  }
400  }
401  else if( ( dims.size() == 3 ) && vertex_data && dimIndex == 2 && mbtype == MB_TYPE_DOUBLE &&
402  use_time ) // the last one is the vertex
403  {
404  // the case when the last dim is ncol (for homme type mesh))
405  size_t dim0, dim1; // dim 2 is ncol..
406  char recname0[NC_MAX_NAME + 1], recname1[NC_MAX_NAME + 1];
407  fail = nc_inq_dim( ncFile, dims[0], recname0, &dim0 );
408  if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
409  fail = nc_inq_dim( ncFile, dims[1], recname1, &dim1 );
410  if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
411  std::string name0( recname0 );
412  std::string name1( recname1 );
413  std::string timestr( "time" );
414  size_t start[3] = { 0, 0, 0 };
415  size_t count[3] = { 1, 0, 0 };
416  count[1] = dim1;
417  count[2] = nodes.size();
418  // create a few tags, with name inserted
419  if( name0.compare( "time" ) == 0 )
420  {
421 
422  std::vector< double > dvalues( dim1 * count[2] );
423  std::vector< float > fvalues( dim1 * count[2] );
424  std::vector< double > transp( dim1 * count[2] );
425  // get the variable time
426  for( size_t k = 0; k < dim0; k++ )
427  {
428  // create a tag for each time, and
429  std::stringstream tag_name;
430  tag_name << variable_name << "_t" << times[k];
431  std::vector< double > defvals( dim1, 0. );
432  rval = mb->tag_get_handle( tag_name.str().c_str(), (int)dim1, mbtype, newTag,
433  MB_TAG_CREAT | MB_TAG_DENSE, &defvals[0] );MB_CHK_SET_ERR( rval, "can't define new tag" );
434  start[0] = k;
435 
436  if( float_var )
437  {
438  fail = nc_get_vara_float( ncFile, nc_var, start, count, &fvalues[0] );
439  // now arrange them in a tag, transpose data
440  for( size_t ii = 0; ii < dim1; ii++ )
441  for( size_t j = 0; j < count[2]; j++ )
442  {
443  transp[j * dim1 + ii] = fvalues[ii * count[2] + j];
444  }
445  }
446  else // double
447  {
448  fail = nc_get_vara_double( ncFile, nc_var, start, count, &dvalues[0] );
449  // now arrange them in a tag, transpose data
450  for( size_t ii = 0; ii < dim1; ii++ )
451  for( size_t j = 0; j < count[2]; j++ )
452  {
453  transp[j * dim1 + ii] = dvalues[ii * count[2] + j];
454  }
455  }
456  for( size_t ii = 0; ii < nodes.size(); ii++ )
457  {
458  EntityHandle vh = vGidHandle[ii + 1];
459  rval = mb->tag_set_data( newTag, &vh, 1, &transp[ii * dim1] );MB_CHK_SET_ERR( rval, "can't set tag on nodes" );
460  }
461  }
462  }
463  }
464  rval = mb->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
465  std::cout << " wrote file " << outfile << "\n";
466 
467  // now, if s option, load the coarse mesh and put data on each element, according to a matrix
468  if( !sefile_name.empty() )
469  {
470  // load the file, check for GLOBAL_DOFS tag, and create a new tag with the data associated
471  Core* mb2 = new Core();
472  rval = mb2->load_file( sefile_name.c_str() );MB_CHK_SET_ERR( rval, "can't load spectral element file" );
473  std::cout << " loaded spectral file " << sefile_name << "\n";
474  // look for GLOBAL_DOFS tag
475  Tag gdofeTag;
476  rval = mb2->tag_get_handle( "GLOBAL_DOFS", gdofeTag );MB_CHK_SET_ERR( rval, "file does not have GLOBAL_DOFS file" );
477  int sizeTag;
478  rval = mb2->tag_get_length( gdofeTag, sizeTag );MB_CHK_SET_ERR( rval, "can't get size of tag" );
479  int np = (int)sqrt( 1.0 * sizeTag );
480  std::cout << " size of tag: " << sizeTag << " np = " << np << "\n";
481  std::vector< int > gdofs;
482  Range cells2;
483  rval = mb2->get_entities_by_dimension( 0, 2, cells2 );MB_CHK_SET_ERR( rval, "can't get cells on spectral mesh" );
484  gdofs.resize( cells2.size() * sizeTag );
485  rval = mb2->tag_get_data( gdofeTag, cells2, &gdofs[0] );MB_CHK_SET_ERR( rval, "can't get global dofs tag" );
486  // create a new tag for element data arranged as DFIELD
487 
488  std::vector< double > dfield;
489  dfield.resize( sizeTag, 0.0 );
490  Tag newTag2;
491  rval = mb2->tag_get_handle( variable_name.c_str(), (int)sizeTag, mbtype, newTag2, MB_TAG_CREAT | MB_TAG_DENSE,
492  &dfield[0] );MB_CHK_SET_ERR( rval, "can't define new tag" );
493 
494  int i1 = 0; // index in the gdofs array, per element
495 
496  // get the tag values from the other moab core, for newTag
497  int dataTagLen;
498  rval = mb->tag_get_length( newTag, dataTagLen );MB_CHK_SET_ERR( rval, "can't get size of newTag" );
499  //
500  std::vector< double > oldData;
501  oldData.resize( dataTagLen * nodes.size() ); //
502 
503  // get the "old" values
504  rval = mb->tag_get_data( newTag, nodes, &oldData[0] );MB_CHK_SET_ERR( rval, "can't get old values" );
505  for( Range::iterator it = cells2.begin(); it != cells2.end(); it++ )
506  {
507  EntityHandle cel = *it;
508  // gdofs per element are gdofs[i:i+np*np];
509  for( int k = 0; k < sizeTag; k++ )
510  {
511  int gdof = gdofs[i1 + k];
512  EntityHandle node = vGidHandle[gdof];
513  int index = nodes.index( node );
514  // get last layer of data
515  double val = oldData[index * dataTagLen + dataTagLen - 1]; //
516  dfield[k] = val;
517  }
518  i1 = i1 + sizeTag;
519  rval = mb2->tag_set_data( newTag2, &cel, 1, &dfield[0] );MB_CHK_SET_ERR( rval, "can't set new tag" );
520  }
521 
522  // write the appended file with the new field:
523  rval = mb2->write_file( "atm2.h5m" );MB_CHK_SET_ERR( rval, "can't write new spectral file" );
524  std::cout << " wrote file atm2.h5m \n";
525  }
526  return 0;
527 }