Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
WriteNC.cpp
Go to the documentation of this file.
1 #include "WriteNC.hpp"
2 #include "moab/CN.hpp"
3 #include "MBTagConventions.hpp"
5 #include "moab/Interface.hpp"
6 #include "moab/Range.hpp"
8 #include "moab/FileOptions.hpp"
9 #include "NCWriteHelper.hpp"
10 
11 #include <fstream>
12 #include <map>
13 #include <set>
14 
15 #include <iostream>
16 #include <sstream>
17 
18 #ifdef WIN32
19 #ifdef size_t
20 #undef size_t
21 #endif
22 #endif
23 
24 namespace moab
25 {
26 
28 {
29  return new WriteNC( iface );
30 }
31 
33  : mbImpl( impl ), dbgOut( stderr ),
34 #ifdef MOAB_HAVE_MPI
35  myPcomm( NULL ),
36 #endif
37  noMesh( false ), noVars( false ), append( false ), mGlobalIdTag( 0 ), isParallel( false ), myHelper( NULL )
38 {
39  assert( impl != NULL );
41 }
42 
44 {
46  if( myHelper != NULL ) delete myHelper;
47 }
48 
49 //! Writes out a file
50 ErrorCode WriteNC::write_file( const char* file_name,
51  const bool overwrite,
52  const FileOptions& options,
53  const EntityHandle* file_set,
54  const int num_set,
55  const std::vector< std::string >&,
56  const Tag*,
57  int,
58  int )
59 {
60  ErrorCode rval;
61  // See if opts has variable(s) specified
62  std::vector< std::string > var_names;
63  std::vector< std::string > desired_names;
64  std::vector< int > tstep_nums;
65  std::vector< double > tstep_vals;
66 
67  // Get and cache predefined tag handles
69 
70  // num set has to be 1, we will write only one set, the original file set used to load
71  if( num_set != 1 ) MB_SET_ERR( MB_FAILURE, "We should write only one set (the file set used to read data into)" );
72 
73  rval = parse_options( options, var_names, desired_names, tstep_nums, tstep_vals );MB_CHK_SET_ERR( rval, "Trouble parsing option string" );
74 
75  // Important to create some data that will be used to write the file; dimensions, variables, etc
76  // new variables still need to have some way of defining their dimensions
77  // maybe it will be passed as write options
78  rval = process_conventional_tags( *file_set );MB_CHK_SET_ERR( rval, "Trouble processing conventional tags" );
79 
80  // Create or append the file
81  if( append )
82  dbgOut.tprintf( 1, "opening file %s for appending \n", file_name );
83  else
84  dbgOut.tprintf( 1, "creating file %s\n", file_name );
85  fileName = file_name;
86  int success;
87 
88  if( append )
89  {
90  int omode = NC_WRITE;
91 #ifdef MOAB_HAVE_PNETCDF
92  if( isParallel )
93  success = NCFUNC( open )( myPcomm->proc_config().proc_comm(), file_name, omode, MPI_INFO_NULL, &fileId );
94  else
95  success = NCFUNC( open )( MPI_COMM_SELF, file_name, omode, MPI_INFO_NULL, &fileId );
96 #else
97  // This is a regular netcdf file, open in write mode
98  success = NCFUNC( open )( file_name, omode, &fileId );
99 #endif
100  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble opening file " << file_name << " for appending" );
101  }
102  else
103  { // Case when the file is new, will be overwritten, most likely
104  int cmode = overwrite ? NC_CLOBBER : NC_NOCLOBBER;
105 #ifdef MOAB_HAVE_PNETCDF
106  if( isParallel )
107  success = NCFUNC( create )( myPcomm->proc_config().proc_comm(), file_name, cmode, MPI_INFO_NULL, &fileId );
108  else
109  success = NCFUNC( create )( MPI_COMM_SELF, file_name, cmode, MPI_INFO_NULL, &fileId );
110 #else
111  // This is a regular netcdf file
112  success = NCFUNC( create )( file_name, cmode, &fileId );
113 #endif
114  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble creating file " << file_name << " for writing" );
115  }
116 
117  if( NULL != myHelper ) delete myHelper;
118 
119  // Get appropriate helper instance for WriteNC class based on some info in the file set
120  myHelper = NCWriteHelper::get_nc_helper( this, fileId, options, *file_set );
121  if( NULL == myHelper )
122  {
123  MB_SET_ERR( MB_FAILURE, "Failed to get NCWriteHelper class instance" );
124  }
125 
126  rval = myHelper->collect_mesh_info();MB_CHK_SET_ERR( rval, "Trouble collecting mesh information" );
127 
128  rval = myHelper->collect_variable_data( var_names, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble collecting variable data" );
129 
130  rval = myHelper->init_file( var_names, desired_names, append );MB_CHK_SET_ERR( rval, "Trouble initializing file" );
131 
132  rval = myHelper->write_values( var_names, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble writing values to file" );
133 
134  success = NCFUNC( close )( fileId );
135  if( success ) MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
136 
137  return MB_SUCCESS;
138 }
139 
141  std::vector< std::string >& var_names,
142  std::vector< std::string >& desired_names,
143  std::vector< int >& tstep_nums,
144  std::vector< double >& tstep_vals )
145 {
146  int tmpval;
147  if( MB_SUCCESS == opts.get_int_option( "DEBUG_IO", 1, tmpval ) )
148  {
149  dbgOut.set_verbosity( tmpval );
150  dbgOut.set_prefix( "NCWrite" );
151  }
152 
153  ErrorCode rval = opts.get_strs_option( "VARIABLE", var_names );
154  if( MB_TYPE_OUT_OF_RANGE == rval )
155  noVars = true;
156  else
157  noVars = false;
158 
159  rval = opts.get_strs_option( "RENAME", desired_names );
160  if( MB_ENTITY_NOT_FOUND == rval )
161  {
162  if( !noVars )
163  {
164  desired_names.resize( var_names.size() );
165  std::copy( var_names.begin(), var_names.end(), desired_names.begin() );
166  }
167  }
168  // Either way
169  assert( desired_names.size() == var_names.size() );
170 
171  opts.get_ints_option( "TIMESTEP", tstep_nums );
172  opts.get_reals_option( "TIMEVAL", tstep_vals );
173  rval = opts.get_null_option( "NOMESH" );
174  if( MB_SUCCESS == rval ) noMesh = true;
175 
176  rval = opts.get_null_option( "APPEND" );
177  if( MB_SUCCESS == rval ) append = true;
178 
179  if( 2 <= dbgOut.get_verbosity() )
180  {
181  if( !var_names.empty() )
182  {
183  std::cerr << "Variables requested: ";
184  for( unsigned int i = 0; i < var_names.size(); i++ )
185  std::cerr << var_names[i];
186  std::cerr << std::endl;
187  }
188  if( !tstep_nums.empty() )
189  {
190  std::cerr << "Timesteps requested: ";
191  for( unsigned int i = 0; i < tstep_nums.size(); i++ )
192  std::cerr << tstep_nums[i];
193  std::cerr << std::endl;
194  }
195  if( !tstep_vals.empty() )
196  {
197  std::cerr << "Time vals requested: ";
198  for( unsigned int i = 0; i < tstep_vals.size(); i++ )
199  std::cerr << tstep_vals[i];
200  std::cerr << std::endl;
201  }
202  }
203 
204 // FIXME: copied from ReadNC, may need revise
205 #ifdef MOAB_HAVE_MPI
206  isParallel = ( opts.match_option( "PARALLEL", "WRITE_PART" ) != MB_ENTITY_NOT_FOUND );
207 
208  if( !isParallel )
209  // Return success here, since rval still has _NOT_FOUND from not finding option
210  // in this case, myPcomm will be NULL, so it can never be used; always check for isParallel
211  // before any use for myPcomm
212  return MB_SUCCESS;
213 
214  int pcomm_no = 0;
215  rval = opts.get_int_option( "PARALLEL_COMM", pcomm_no );
216  if( MB_TYPE_OUT_OF_RANGE == rval )
217  {
218  MB_SET_ERR( rval, "Invalid value for PARALLEL_COMM option" );
219  }
220 
221  myPcomm = ParallelComm::get_pcomm( mbImpl, pcomm_no );
222  if( 0 == myPcomm )
223  {
224  myPcomm = new ParallelComm( mbImpl, MPI_COMM_WORLD );
225  }
226 
227 #ifndef MOAB_HAVE_PNETCDF
228  const int procs = myPcomm->proc_config().proc_size();
229  if( procs > 1 )
230  {
231  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Attempt to launch NC writer in parallel without pnetcdf support" );
232  }
233 #endif
234 
235  const int rank = myPcomm->proc_config().proc_rank();
236  dbgOut.set_rank( rank );
237 #endif
238 
239  return MB_SUCCESS;
240 }
241 
242 // This is the inverse process to create conventional tags
243 // Will look at <pargal_source>/src/core/fileinfo.cpp, init dim, vars, atts
245 {
246  ErrorCode rval;
247 
248  // Start copy
249  Tag dimNamesTag = 0;
250  std::string tag_name = "__DIM_NAMES";
251  const void* data = NULL;
252  int dimNamesSz = 0;
253  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, dimNamesTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
254  rval = mbImpl->tag_get_by_ptr( dimNamesTag, &fileSet, 1, &data, &dimNamesSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
255  const char* p = static_cast< const char* >( data );
256  dbgOut.tprintf( 1, "__DIM_NAMES tag has string length %d\n", dimNamesSz );
257 
258  std::size_t start = 0;
259 
260  Tag dimLensTag = 0;
261  tag_name = "__DIM_LENS";
262  data = NULL;
263  int dimLensSz = 0;
264  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_INTEGER, dimLensTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
265  rval = mbImpl->tag_get_by_ptr( dimLensTag, &fileSet, 1, &data, &dimLensSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
266  const int* int_p = static_cast< const int* >( data );
267  dbgOut.tprintf( 1, "__DIM_LENS tag has %d values\n", dimLensSz );
268 
269  int idxDim = 0;
270  // Dim names are separated by '\0' in the string of __DIM_NAMES tag
271  for( std::size_t i = 0; i != static_cast< std::size_t >( dimNamesSz ); i++ )
272  {
273  if( p[i] == '\0' )
274  {
275  std::string dim_name( &p[start], i - start );
276  int len = int_p[idxDim];
277  dimNames.push_back( dim_name );
278  dimLens.push_back( len );
279  dbgOut.tprintf( 2, "Dimension %s has length %d\n", dim_name.c_str(), len );
280  // FIXME: Need info from moab to set unlimited dimension
281  // Currently assume each file has the same number of time dimensions
282  /*if ((dim_name == "time") || (dim_name == "Time"))
283  insert(dim_name, *(new pcdim(dim_name, len * m_file_names.size(), true)));
284  else
285  insert(dim_name, *(new pcdim(dim_name, len)));*/
286  start = i + 1;
287  idxDim++;
288  }
289  }
290 
291  Tag meshTypeTag = 0;
292  tag_name = "__MESH_TYPE";
293  data = NULL;
294  int meshTypeSz = 0;
295  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, meshTypeTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
296  rval = mbImpl->tag_get_by_ptr( meshTypeTag, &fileSet, 1, &data, &meshTypeSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
297  p = static_cast< const char* >( data );
298  grid_type = std::string( &p[0], meshTypeSz );
299  dbgOut.tprintf( 2, "Mesh type: %s\n", grid_type.c_str() );
300 
301  // Read <__VAR_NAMES_LOCATIONS> tag
302  Tag varNamesLocsTag = 0;
303  tag_name = "__VAR_NAMES_LOCATIONS";
304  data = NULL;
305  int varNamesLocsSz = 0;
306  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_INTEGER, varNamesLocsTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
307  rval = mbImpl->tag_get_by_ptr( varNamesLocsTag, &fileSet, 1, &data, &varNamesLocsSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
308  int_p = static_cast< const int* >( data );
309  std::vector< int > varNamesLocs( varNamesLocsSz );
310  std::copy( int_p, int_p + varNamesLocsSz, varNamesLocs.begin() );
311 
312  Tag varNamesTag = 0;
313  tag_name = "__VAR_NAMES";
314  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, varNamesTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
315  data = NULL;
316  int varNamesSz = 0;
317  rval = mbImpl->tag_get_by_ptr( varNamesTag, &fileSet, 1, &data, &varNamesSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
318  dbgOut.tprintf( 2, "__VAR_NAMES tag has string length %d\n", varNamesSz );
319  p = static_cast< const char* >( data );
320 
321  start = 0;
322  int idxVar = 0;
323  int sz;
324  // Var names are separated by '\0' in the string of __VAR_NAMES tag
325  for( std::size_t i = 0; i != static_cast< std::size_t >( varNamesSz ); i++ )
326  {
327  if( p[i] == '\0' )
328  {
329  std::string var_name( &p[start], i - start );
330 
331  dbgOut.tprintf( 2, "var name: %s index %d \n", var_name.c_str(), idxVar );
332  // Process var name:
333  // This will create/initiate map; we will populate variableDataStruct with info about
334  // dims, tags, etc reference & is important; otherwise variableDataStruct will go out of
335  // scope, and deleted :(
336  VarData& variableDataStruct = varInfo[var_name];
337  variableDataStruct.varName = var_name;
338  variableDataStruct.entLoc = varNamesLocs[idxVar];
339 
340  dbgOut.tprintf( 2, "at var name %s varInfo size %d \n", var_name.c_str(), (int)varInfo.size() );
341 
342  sz = 0;
343  Tag dims_tag = 0;
344  std::string dim_names = "__" + var_name + "_DIMS";
345  rval = mbImpl->tag_get_handle( dim_names.c_str(), 0, MB_TYPE_OPAQUE, dims_tag, MB_TAG_ANY );
346  if( MB_SUCCESS != rval )
347  {
348  if( MB_TAG_NOT_FOUND == rval )
349  {
350  dbgOut.tprintf( 2, "tag : %s not found, continue \n", dim_names.c_str() );
351  start = i + 1;
352  idxVar++;
353  continue;
354  }
355  MB_SET_ERR( rval, "Trouble getting conventional tag " << dim_names );
356  }
357  rval = mbImpl->tag_get_length( dims_tag, sz );MB_CHK_SET_ERR( rval, "Trouble getting size of dimensions for variable " << var_name );
358  sz /= sizeof( Tag ); // The type is MB_TYPE_OPAQUE, but it is a list of tags, so we
359  // need to divide by the size of Tag
360  // sz is used for number of dimension tags in this list
361  dbgOut.tprintf( 2, "var name: %s has %d dimensions \n", var_name.c_str(), sz );
362 
363  variableDataStruct.varDims.resize( sz );
364  const void* ptr = NULL;
365  rval = mbImpl->tag_get_by_ptr( dims_tag, &fileSet, 1, &ptr );
366 
367  const Tag* ptags = static_cast< const moab::Tag* >( ptr );
368  for( std::size_t j = 0; j != static_cast< std::size_t >( sz ); j++ )
369  {
370  std::string dim_name;
371  rval = mbImpl->tag_get_name( ptags[j], dim_name );MB_CHK_SET_ERR( rval, "Trouble getting dimension of variable " << var_name );
372  dbgOut.tprintf( 2, "var name: %s has %s as dimension \n", var_name.c_str(), dim_name.c_str() );
373  std::vector< std::string >::iterator vit = std::find( dimNames.begin(), dimNames.end(), dim_name );
374  if( vit == dimNames.end() )
375  MB_SET_ERR( MB_FAILURE, "Dimension " << dim_name << " not found for variable " << var_name );
376  variableDataStruct.varDims[j] = (int)( vit - dimNames.begin() ); // Will be used for writing
377  // This will have to change to actual file dimension, for writing
378  }
379 
380  // Attributes for this variable
381  std::stringstream ssTagName;
382  ssTagName << "__" << var_name << "_ATTRIBS";
383  tag_name = ssTagName.str();
384  Tag varAttTag = 0;
385  rval =
386  mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, varAttTag, MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
387  const void* varAttPtr = NULL;
388  int varAttSz = 0;
389  rval = mbImpl->tag_get_by_ptr( varAttTag, &fileSet, 1, &varAttPtr, &varAttSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
390  if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Tag retrieved for variable %s\n", tag_name.c_str() );
391 
392  std::string attribString( (char*)varAttPtr, (char*)varAttPtr + varAttSz );
393  if( attribString == "NO_ATTRIBS" )
394  {
395  // This variable has no attributes
396  variableDataStruct.numAtts = 0;
397  }
398  else if( attribString == "DUMMY_VAR" )
399  {
400  // This variable is a dummy coordinate variable
401  variableDataStruct.numAtts = 0;
402  dummyVarNames.insert( variableDataStruct.varName );
403  }
404  else
405  {
406  ssTagName << "_LEN";
407  tag_name = ssTagName.str();
408  Tag varAttLenTag = 0;
409  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_INTEGER, varAttLenTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
410  int varAttLenSz = 0;
411  rval = mbImpl->tag_get_length( varAttLenTag, varAttLenSz );MB_CHK_SET_ERR( rval, "Trouble getting length of conventional tag " << tag_name );
412  std::vector< int > varAttLen( varAttLenSz );
413  rval = mbImpl->tag_get_data( varAttLenTag, &fileSet, 1, &varAttLen[0] );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
414 
415  rval = process_concatenated_attribute( varAttPtr, varAttSz, varAttLen, variableDataStruct.varAtts );MB_CHK_SET_ERR( rval, "Trouble processing attributes of variable " << var_name );
416 
417  if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Tag metadata for variable %s\n", tag_name.c_str() );
418  }
419  // End attribute
420 
421  start = i + 1;
422  idxVar++;
423  } // if (p[i] == '\0')
424  }
425 
426  // Global attributes
427  tag_name = "__GLOBAL_ATTRIBS";
428  Tag globalAttTag = 0;
429  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, globalAttTag, MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
430  std::vector< int > gattLen;
431 
432  const void* gattptr = NULL;
433  int globalAttSz = 0;
434  rval = mbImpl->tag_get_by_ptr( globalAttTag, &fileSet, 1, &gattptr, &globalAttSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
435 
436  if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Tag value retrieved for %s size %d\n", tag_name.c_str(), globalAttSz );
437 
438  // <__GLOBAL_ATTRIBS_LEN>
439  tag_name = "__GLOBAL_ATTRIBS_LEN";
440  Tag globalAttLenTag = 0;
441 
442  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_INTEGER, globalAttLenTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
443  int sizeGAtt = 0;
444  rval = mbImpl->tag_get_length( globalAttLenTag, sizeGAtt );MB_CHK_SET_ERR( rval, "Trouble getting length of conventional tag " << tag_name );
445  gattLen.resize( sizeGAtt );
446  rval = mbImpl->tag_get_data( globalAttLenTag, &fileSet, 1, &gattLen[0] );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
447  if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Tag retrieved for variable %s\n", tag_name.c_str() );
448 
449  rval = process_concatenated_attribute( gattptr, globalAttSz, gattLen, globalAtts );MB_CHK_SET_ERR( rval, "Trouble processing global attributes" );
450 
451  return MB_SUCCESS;
452 }
453 
454 // Reverse process from create_attrib_string
456  int attSz,
457  std::vector< int >& attLen,
458  std::map< std::string, AttData >& attributes )
459 {
460  std::size_t start = 0;
461  std::size_t att_counter = 0;
462  std::string concatString( (char*)attPtr, (char*)attPtr + attSz );
463 
464  for( std::size_t i = 0; i != (size_t)attSz; i++ )
465  {
466  if( concatString[i] == '\0' )
467  {
468  std::string att_name( &concatString[start], i - start );
469  start = i + 1;
470  while( concatString[i] != ';' )
471  ++i;
472  std::string data_type( &concatString[start], i - start );
473  ++i;
474  start = i;
475  i = attLen[att_counter];
476  if( concatString[i] != ';' ) MB_SET_ERR( MB_FAILURE, "Error parsing attributes" );
477 
478  std::string data_val( &concatString[start], i - start );
479  start = i + 1;
480 
481  AttData& attrib = attributes[att_name];
482  attrib.attValue = data_val;
483  attrib.attLen = data_val.size();
484 
485  if( data_type == "char" )
486  attrib.attDataType = NC_CHAR;
487  else if( data_type == "double" )
488  attrib.attDataType = NC_DOUBLE;
489  else if( data_type == "float" )
490  attrib.attDataType = NC_FLOAT;
491  else if( data_type == "int" )
492  attrib.attDataType = NC_INT;
493  else if( data_type == "short" )
494  attrib.attDataType = NC_SHORT;
495 
496  ++att_counter;
497  dbgOut.tprintf( 2, " Process attribute %s with value %s \n", att_name.c_str(), data_val.c_str() );
498  }
499  }
500 
501  return MB_SUCCESS;
502 }
503 
504 } // namespace moab