Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadNC.cpp
Go to the documentation of this file.
1 #include "ReadNC.hpp"
2 #include "NCHelper.hpp"
3 
4 #include "moab/ReadUtilIface.hpp"
5 #include "MBTagConventions.hpp"
6 #include "moab/FileOptions.hpp"
7 
8 namespace moab
9 {
10 
12 {
13  return new ReadNC( iface );
14 }
15 
17  : mbImpl( impl ), fileId( -1 ), mGlobalIdTag( 0 ), mpFileIdTag( NULL ), dbgOut( stderr ), isParallel( false ),
18  partMethod( ScdParData::ALLJORKORI ), scdi( NULL ),
19 #ifdef MOAB_HAVE_MPI
20  myPcomm( NULL ),
21 #endif
22  noMesh( false ), noVars( false ), spectralMesh( false ), noMixedElements( false ), noEdges( false ), culling(true),
23  repartition(false), gatherSetRank( -1 ), tStepBase( -1 ), trivialPartitionShift( 0 ), myHelper( NULL )
24 {
25  assert( impl != NULL );
27 }
28 
30 {
32  if( myHelper != NULL ) delete myHelper;
33 }
34 
35 ErrorCode ReadNC::load_file( const char* file_name,
36  const EntityHandle* file_set,
37  const FileOptions& opts,
38  const ReaderIface::SubsetList* /*subset_list*/,
39  const Tag* file_id_tag )
40 {
41  // See if opts has variable(s) specified
42  std::vector< std::string > var_names;
43  std::vector< int > tstep_nums;
44  std::vector< double > tstep_vals;
45 
46  // Get and cache predefined tag handles
48  // Store the pointer to the tag; if not null, set when global id tag
49  // is set too, with the same data, duplicated
50  mpFileIdTag = file_id_tag;
51 
52  ErrorCode rval = parse_options( opts, var_names, tstep_nums, tstep_vals );MB_CHK_SET_ERR( rval, "Trouble parsing option string" );
53 
54  // Open the file
55  dbgOut.tprintf( 1, "Opening file %s\n", file_name );
56  fileName = std::string( file_name );
57  int success;
58 
59 #ifdef MOAB_HAVE_PNETCDF
60  if( isParallel )
61  success = NCFUNC( open )( myPcomm->proc_config().proc_comm(), file_name, 0, MPI_INFO_NULL, &fileId );
62  else
63  success = NCFUNC( open )( MPI_COMM_SELF, file_name, 0, MPI_INFO_NULL, &fileId );
64 #else
65  success = NCFUNC( open )( file_name, 0, &fileId );
66 #endif
67  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble opening file " << file_name );
68 
69  // Read the header (num dimensions, dimensions, num variables, global attribs)
70  rval = read_header();MB_CHK_SET_ERR( rval, "Trouble reading file header" );
71 
72  // Make sure there's a file set to put things in
73  EntityHandle tmp_set;
74  if( noMesh && !file_set )
75  {
76  MB_SET_ERR( MB_FAILURE, "NOMESH option requires non-NULL file set on input" );
77  }
78  else if( !file_set || ( file_set && *file_set == 0 ) )
79  {
80  rval = mbImpl->create_meshset( MESHSET_SET, tmp_set );MB_CHK_SET_ERR( rval, "Trouble creating file set" );
81  }
82  else
83  tmp_set = *file_set;
84 
85  // Get the scd interface
86  scdi = NULL;
87  rval = mbImpl->query_interface( scdi );
88  if( NULL == scdi ) return MB_FAILURE;
89 
90  if( NULL != myHelper ) delete myHelper;
91 
92  // Get appropriate NC helper instance based on information read from the header
93  myHelper = NCHelper::get_nc_helper( this, fileId, opts, tmp_set );
94  if( NULL == myHelper )
95  {
96  MB_SET_ERR( MB_FAILURE, "Failed to get NCHelper class instance" );
97  }
98 
99  // Initialize mesh values
100  rval = myHelper->init_mesh_vals();MB_CHK_SET_ERR( rval, "Trouble initializing mesh values" );
101 
102  // Check existing mesh from last read
103  if( noMesh && !noVars )
104  {
105  rval = myHelper->check_existing_mesh();MB_CHK_SET_ERR( rval, "Trouble checking mesh from last read" );
106  }
107 
108  // Create some conventional tags, e.g. __NUM_DIMS
109  // For multiple reads to a specified file set, we assume a single file, or a series of
110  // files with separated timesteps. Keep a flag on the file set to prevent conventional
111  // tags from being created again on a second read
112  Tag convTagsCreated = 0;
113  int def_val = 0;
114  rval = mbImpl->tag_get_handle( "__CONV_TAGS_CREATED", 1, MB_TYPE_INTEGER, convTagsCreated,
115  MB_TAG_SPARSE | MB_TAG_CREAT, &def_val );MB_CHK_SET_ERR( rval, "Trouble getting _CONV_TAGS_CREATED tag" );
116  int create_conv_tags_flag = 0;
117  rval = mbImpl->tag_get_data( convTagsCreated, &tmp_set, 1, &create_conv_tags_flag );
118  // The first read to the file set
119  if( 0 == create_conv_tags_flag )
120  {
121  // Read dimensions (coordinate variables) by default to create tags like __<var_name>_DIMS
122  // This is done only once (assume that all files read to the file set have the same
123  // dimensions)
124  rval = myHelper->read_variables( dimNames, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading dimensions" );
125 
126  rval = myHelper->create_conventional_tags( tstep_nums );MB_CHK_SET_ERR( rval, "Trouble creating NC conventional tags" );
127 
128  create_conv_tags_flag = 1;
129  rval = mbImpl->tag_set_data( convTagsCreated, &tmp_set, 1, &create_conv_tags_flag );MB_CHK_SET_ERR( rval, "Trouble setting data to _CONV_TAGS_CREATED tag" );
130  }
131  // Another read to the file set
132  else
133  {
134  if( tStepBase > -1 )
135  {
136  // If timesteps spread across files, merge time values read
137  // from current file to existing time tag
138  rval = myHelper->update_time_tag_vals();MB_CHK_SET_ERR( rval, "Trouble updating time tag values" );
139  }
140  }
141 
142  // Create mesh vertex/edge/face sequences
143  Range faces;
144  if( !noMesh )
145  {
146  rval = myHelper->create_mesh( faces );MB_CHK_SET_ERR( rval, "Trouble creating mesh" );
147  }
148 
149  // Read specified variables onto grid
150  if( !noVars )
151  {
152  if( var_names.empty() )
153  {
154  // If VARIABLE option is missing, read all variables
155  rval = myHelper->read_variables( var_names, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading all variables" );
156  }
157  else
158  {
159  // Exclude dimensions that are read to the file set by default
160  std::vector< std::string > non_dim_var_names;
161  for( unsigned int i = 0; i < var_names.size(); i++ )
162  {
163  if( std::find( dimNames.begin(), dimNames.end(), var_names[i] ) == dimNames.end() )
164  non_dim_var_names.push_back( var_names[i] );
165  }
166 
167  if( !non_dim_var_names.empty() )
168  {
169  rval = myHelper->read_variables( non_dim_var_names, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading specified variables" );
170  }
171  }
172  }
173 
174 #ifdef MOAB_HAVE_MPI
175  // Create partition set, and populate with elements
176  if( isParallel )
177  {
178  // Write partition tag name on partition set
179  Tag part_tag = myPcomm->partition_tag();
180  int dum_rank = myPcomm->proc_config().proc_rank();
181  // the tmp_set is the file_set
182  rval = mbImpl->tag_set_data( part_tag, &tmp_set, 1, &dum_rank );MB_CHK_SET_ERR( rval, "Trouble writing partition tag name on partition set" );
183  }
184 #endif
185 
187  scdi = NULL;
188 
189  // Close the file
190  success = NCFUNC( close )( fileId );
191  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
192 
193  return MB_SUCCESS;
194 }
195 
197  std::vector< std::string >& var_names,
198  std::vector< int >& tstep_nums,
199  std::vector< double >& tstep_vals )
200 {
201  int tmpval;
202  if( MB_SUCCESS == opts.get_int_option( "DEBUG_IO", 1, tmpval ) )
203  {
204  dbgOut.set_verbosity( tmpval );
205  dbgOut.set_prefix( "NC " );
206  }
207 
208  ErrorCode rval = opts.get_strs_option( "VARIABLE", var_names );
209  if( MB_TYPE_OUT_OF_RANGE == rval )
210  noVars = true;
211  else
212  noVars = false;
213 
214  opts.get_ints_option( "TIMESTEP", tstep_nums );
215  opts.get_reals_option( "TIMEVAL", tstep_vals );
216 
217  rval = opts.get_null_option( "NOMESH" );
218  if( MB_SUCCESS == rval ) noMesh = true;
219 
220  rval = opts.get_null_option( "SPECTRAL_MESH" );
221  if( MB_SUCCESS == rval ) spectralMesh = true;
222 
223  rval = opts.get_null_option( "NO_MIXED_ELEMENTS" );
224  if( MB_SUCCESS == rval ) noMixedElements = true;
225 
226  rval = opts.get_null_option( "NO_EDGES" );
227  if( MB_SUCCESS == rval ) noEdges = true;
228 
229  rval = opts.get_null_option( "NO_CULLING" ); // used now only for domain nc convention
230  if( MB_SUCCESS == rval ) culling = false;
231 
232  rval = opts.get_null_option( "REPARTITION" ); // used now only for domain nc, to repartition with zoltan
233  if( MB_SUCCESS == rval ) repartition = true;
234 
235  if( 2 <= dbgOut.get_verbosity() )
236  {
237  if( !var_names.empty() )
238  {
239  std::cerr << "Variables requested: ";
240  for( unsigned int i = 0; i < var_names.size(); i++ )
241  std::cerr << var_names[i];
242  std::cerr << std::endl;
243  }
244 
245  if( !tstep_nums.empty() )
246  {
247  std::cerr << "Timesteps requested: ";
248  for( unsigned int i = 0; i < tstep_nums.size(); i++ )
249  std::cerr << tstep_nums[i];
250  std::cerr << std::endl;
251  }
252 
253  if( !tstep_vals.empty() )
254  {
255  std::cerr << "Time vals requested: ";
256  for( unsigned int i = 0; i < tstep_vals.size(); i++ )
257  std::cerr << tstep_vals[i];
258  std::cerr << std::endl;
259  }
260  }
261 
262  rval = opts.get_int_option( "GATHER_SET", 0, gatherSetRank );
263  if( MB_TYPE_OUT_OF_RANGE == rval )
264  {
265  MB_SET_ERR( rval, "Invalid value for GATHER_SET option" );
266  }
267 
268  rval = opts.get_int_option( "TIMESTEPBASE", 0, tStepBase );
269  if( MB_TYPE_OUT_OF_RANGE == rval )
270  {
271  MB_SET_ERR( rval, "Invalid value for TIMESTEPBASE option" );
272  }
273 
274  rval = opts.get_int_option( "TRIVIAL_PARTITION_SHIFT", 1, trivialPartitionShift );
275  if( MB_TYPE_OUT_OF_RANGE == rval )
276  {
277  MB_SET_ERR( rval, "Invalid value for TRIVIAL_PARTITION_SHIFT option" );
278  }
279 
280 #ifdef MOAB_HAVE_MPI
281  isParallel = ( opts.match_option( "PARALLEL", "READ_PART" ) != MB_ENTITY_NOT_FOUND );
282 
283  if( !isParallel )
284  // Return success here, since rval still has _NOT_FOUND from not finding option
285  // in this case, myPcomm will be NULL, so it can never be used; always check for isParallel
286  // before any use for myPcomm
287  return MB_SUCCESS;
288 
289  int pcomm_no = 0;
290  rval = opts.get_int_option( "PARALLEL_COMM", pcomm_no );
291  if( MB_TYPE_OUT_OF_RANGE == rval )
292  {
293  MB_SET_ERR( rval, "Invalid value for PARALLEL_COMM option" );
294  }
295  myPcomm = ParallelComm::get_pcomm( mbImpl, pcomm_no );
296  if( 0 == myPcomm )
297  {
298  myPcomm = new ParallelComm( mbImpl, MPI_COMM_WORLD );
299  }
300  const int rank = myPcomm->proc_config().proc_rank();
301  dbgOut.set_rank( rank );
302 
303  int dum;
304  rval = opts.match_option( "PARTITION_METHOD", ScdParData::PartitionMethodNames, dum );
305  if( MB_FAILURE == rval )
306  {
307  MB_SET_ERR( rval, "Unknown partition method specified" );
308  }
309  else if( MB_ENTITY_NOT_FOUND == rval )
311  else
312  partMethod = dum;
313 #endif
314 
315  return MB_SUCCESS;
316 }
317 
319 {
320  dbgOut.tprint( 1, "Reading header...\n" );
321 
322  // Get the global attributes
323  int numgatts;
324  int success;
325  success = NCFUNC( inq_natts )( fileId, &numgatts );
326  if( success ) MB_SET_ERR( MB_FAILURE, "Couldn't get number of global attributes" );
327 
328  // Read attributes into globalAtts
329  ErrorCode result = get_attributes( NC_GLOBAL, numgatts, globalAtts );MB_CHK_SET_ERR( result, "Trouble getting global attributes" );
330  dbgOut.tprintf( 1, "Read %u attributes\n", (unsigned int)globalAtts.size() );
331 
332  // Read in dimensions into dimNames and dimLens
333  result = get_dimensions( fileId, dimNames, dimLens );MB_CHK_SET_ERR( result, "Trouble getting dimensions" );
334  dbgOut.tprintf( 1, "Read %u dimensions\n", (unsigned int)dimNames.size() );
335 
336  // Read in variables into varInfo
337  result = get_variables();MB_CHK_SET_ERR( result, "Trouble getting variables" );
338  dbgOut.tprintf( 1, "Read %u variables\n", (unsigned int)varInfo.size() );
339 
340  return MB_SUCCESS;
341 }
342 
343 ErrorCode ReadNC::get_attributes( int var_id, int num_atts, std::map< std::string, AttData >& atts, const char* prefix )
344 {
345  char dum_name[120];
346 
347  for( int i = 0; i < num_atts; i++ )
348  {
349  // Get the name
350  int success = NCFUNC( inq_attname )( fileId, var_id, i, dum_name );
351  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting attribute name" );
352 
353  AttData& data = atts[std::string( dum_name )];
354  data.attName = std::string( dum_name );
355  success = NCFUNC( inq_att )( fileId, var_id, dum_name, &data.attDataType, &data.attLen );
356  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting info for attribute " << data.attName );
357  data.attVarId = var_id;
358 
359  dbgOut.tprintf( 2, "%sAttribute %s: length=%u, varId=%d, type=%d\n", ( prefix ? prefix : "" ),
360  data.attName.c_str(), (unsigned int)data.attLen, data.attVarId, data.attDataType );
361  }
362 
363  return MB_SUCCESS;
364 }
365 
366 ErrorCode ReadNC::get_dimensions( int file_id, std::vector< std::string >& dim_names, std::vector< int >& dim_lens )
367 {
368  // Get the number of dimensions
369  int num_dims;
370  int success = NCFUNC( inq_ndims )( file_id, &num_dims );
371  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of dimensions" );
372 
373  if( num_dims > NC_MAX_DIMS )
374  {
375  MB_SET_ERR( MB_FAILURE,
376  "ReadNC: File contains " << num_dims << " dims but NetCDF library supports only " << NC_MAX_DIMS );
377  }
378 
379  char dim_name[NC_MAX_NAME + 1];
380  NCDF_SIZE dim_len;
381  dim_names.resize( num_dims );
382  dim_lens.resize( num_dims );
383 
384  for( int i = 0; i < num_dims; i++ )
385  {
386  success = NCFUNC( inq_dim )( file_id, i, dim_name, &dim_len );
387  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting dimension info" );
388 
389  dim_names[i] = std::string( dim_name );
390  dim_lens[i] = dim_len;
391 
392  dbgOut.tprintf( 2, "Dimension %s, length=%u\n", dim_name, (unsigned int)dim_len );
393  }
394 
395  return MB_SUCCESS;
396 }
397 
399 {
400  // First cache the number of time steps
401  std::vector< std::string >::iterator vit = std::find( dimNames.begin(), dimNames.end(), "time" );
402  if( vit == dimNames.end() ) vit = std::find( dimNames.begin(), dimNames.end(), "t" );
403 
404  int ntimes = 0;
405  if( vit != dimNames.end() ) ntimes = dimLens[vit - dimNames.begin()];
406  if( !ntimes ) ntimes = 1;
407 
408  // Get the number of variables
409  int num_vars;
410  int success = NCFUNC( inq_nvars )( fileId, &num_vars );
411  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of variables" );
412 
413  if( num_vars > NC_MAX_VARS )
414  {
415  MB_SET_ERR( MB_FAILURE,
416  "ReadNC: File contains " << num_vars << " vars but NetCDF library supports only " << NC_MAX_VARS );
417  }
418 
419  char var_name[NC_MAX_NAME + 1];
420  int var_ndims;
421 
422  for( int i = 0; i < num_vars; i++ )
423  {
424  // Get the name first, so we can allocate a map iterate for this var
425  success = NCFUNC( inq_varname )( fileId, i, var_name );
426  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting variable name" );
427  VarData& data = varInfo[std::string( var_name )];
428  data.varName = std::string( var_name );
429  data.varId = i;
430  data.varTags.resize( ntimes, 0 );
431 
432  // Get the data type
433  success = NCFUNC( inq_vartype )( fileId, i, &data.varDataType );
434  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting data type for variable " << data.varName );
435 
436  // Get the number of dimensions, then the dimensions
437  success = NCFUNC( inq_varndims )( fileId, i, &var_ndims );
438  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of dims for variable " << data.varName );
439  data.varDims.resize( var_ndims );
440 
441  success = NCFUNC( inq_vardimid )( fileId, i, &data.varDims[0] );
442  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting dimensions for variable " << data.varName );
443 
444  // Finally, get the number of attributes, then the attributes
445  success = NCFUNC( inq_varnatts )( fileId, i, &data.numAtts );
446  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of dims for variable " << data.varName );
447 
448  // Print debug info here so attribute info comes afterwards
449  dbgOut.tprintf( 2, "Variable %s: Id=%d, numAtts=%d, datatype=%d, num_dims=%u\n", data.varName.c_str(),
450  data.varId, data.numAtts, data.varDataType, (unsigned int)data.varDims.size() );
451 
452  ErrorCode rval = get_attributes( i, data.numAtts, data.varAtts, " " );MB_CHK_SET_ERR( rval, "Trouble getting attributes for variable " << data.varName );
453  }
454 
455  return MB_SUCCESS;
456 }
457 
459  const char*,
460  const FileOptions&,
461  std::vector< int >&,
462  const SubsetList* )
463 {
464  return MB_FAILURE;
465 }
466 
467 } // namespace moab