Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
NCHelperHOMME.cpp
Go to the documentation of this file.
1 #include "NCHelperHOMME.hpp"
2 #include "moab/ReadUtilIface.hpp"
3 #include "moab/FileOptions.hpp"
5 
6 #include <cmath>
7 
8 namespace moab
9 {
10 
11 NCHelperHOMME::NCHelperHOMME( ReadNC* readNC, int fileId, const FileOptions& opts, EntityHandle fileSet )
12  : UcdNCHelper( readNC, fileId, opts, fileSet ), _spectralOrder( -1 ), connectId( -1 ), isConnFile( false )
13 {
14  // Calculate spectral order
15  std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "np" );
16  if( attIt != readNC->globalAtts.end() )
17  {
18  int success = NCFUNC( get_att_int )( readNC->fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
19  &_spectralOrder );
20  if( 0 == success ) _spectralOrder--; // Spectral order is one less than np
21  }
22  else
23  {
24  // As can_read_file() returns true and there is no global attribute "np", it should be a
25  // connectivity file
26  isConnFile = true;
27  _spectralOrder = 3; // Assume np is 4
28  }
29 }
30 
31 bool NCHelperHOMME::can_read_file( ReadNC* readNC, int fileId )
32 {
33  // If global attribute "np" exists then it should be the HOMME grid
34  if( readNC->globalAtts.find( "np" ) != readNC->globalAtts.end() )
35  {
36  // Make sure it is CAM grid
37  std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "source" );
38  if( attIt == readNC->globalAtts.end() ) return false;
39  unsigned int sz = attIt->second.attLen;
40  std::string att_data;
41  att_data.resize( sz + 1 );
42  att_data[sz] = '\000';
43  int success =
44  NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
45  if( success ) return false;
46  if( att_data.find( "CAM" ) == std::string::npos ) return false;
47 
48  return true;
49  }
50  else
51  {
52  // If dimension names "ncol" AND "ncorners" AND "ncells" exist, then it should be the HOMME
53  // connectivity file In this case, the mesh can still be created
54  std::vector< std::string >& dimNames = readNC->dimNames;
55  if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "ncol" ) ) != dimNames.end() ) &&
56  ( std::find( dimNames.begin(), dimNames.end(), std::string( "ncorners" ) ) != dimNames.end() ) &&
57  ( std::find( dimNames.begin(), dimNames.end(), std::string( "ncells" ) ) != dimNames.end() ) )
58  return true;
59  }
60 
61  return false;
62 }
63 
65 {
66  std::vector< std::string >& dimNames = _readNC->dimNames;
67  std::vector< int >& dimLens = _readNC->dimLens;
68  std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
69 
70  ErrorCode rval;
71  unsigned int idx;
72  std::vector< std::string >::iterator vit;
73 
74  // Look for time dimension
75  if( isConnFile )
76  {
77  // Connectivity file might not have time dimension
78  }
79  else
80  {
81  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() )
82  idx = vit - dimNames.begin();
83  else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "t" ) ) != dimNames.end() )
84  idx = vit - dimNames.begin();
85  else
86  {
87  MB_SET_ERR( MB_FAILURE, "Couldn't find 'time' or 't' dimension" );
88  }
89  tDim = idx;
90  nTimeSteps = dimLens[idx];
91  }
92 
93  // Get number of vertices (labeled as number of columns)
94  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ncol" ) ) != dimNames.end() )
95  idx = vit - dimNames.begin();
96  else
97  {
98  MB_SET_ERR( MB_FAILURE, "Couldn't find 'ncol' dimension" );
99  }
100  vDim = idx;
101  nVertices = dimLens[idx];
102 
103  // Set number of cells
104  nCells = nVertices - 2;
105 
106  // Get number of levels
107  if( isConnFile )
108  {
109  // Connectivity file might not have level dimension
110  }
111  else
112  {
113  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lev" ) ) != dimNames.end() )
114  idx = vit - dimNames.begin();
115  else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ilev" ) ) != dimNames.end() )
116  idx = vit - dimNames.begin();
117  else
118  {
119  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension" );
120  }
121  levDim = idx;
122  nLevels = dimLens[idx];
123  }
124 
125  // Store lon values in xVertVals
126  std::map< std::string, ReadNC::VarData >::iterator vmit;
127  if( ( vmit = varInfo.find( "lon" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
128  {
129  rval = read_coordinate( "lon", 0, nVertices - 1, xVertVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lon' variable" );
130  }
131  else
132  {
133  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' variable" );
134  }
135 
136  // Store lat values in yVertVals
137  if( ( vmit = varInfo.find( "lat" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
138  {
139  rval = read_coordinate( "lat", 0, nVertices - 1, yVertVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lat' variable" );
140  }
141  else
142  {
143  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' variable" );
144  }
145 
146  // Store lev values in levVals
147  if( isConnFile )
148  {
149  // Connectivity file might not have level variable
150  }
151  else
152  {
153  if( ( vmit = varInfo.find( "lev" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
154  {
155  rval = read_coordinate( "lev", 0, nLevels - 1, levVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lev' variable" );
156 
157  // Decide whether down is positive
158  char posval[10] = { 0 };
159  int success = NCFUNC( get_att_text )( _fileId, ( *vmit ).second.varId, "positive", posval );
160  if( 0 == success && !strcmp( posval, "down" ) )
161  {
162  for( std::vector< double >::iterator dvit = levVals.begin(); dvit != levVals.end(); ++dvit )
163  ( *dvit ) *= -1.0;
164  }
165  }
166  else
167  {
168  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' variable" );
169  }
170  }
171 
172  // Store time coordinate values in tVals
173  if( isConnFile )
174  {
175  // Connectivity file might not have time variable
176  }
177  else
178  {
179  if( ( vmit = varInfo.find( "time" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
180  {
181  rval = read_coordinate( "time", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 'time' variable" );
182  }
183  else if( ( vmit = varInfo.find( "t" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
184  {
185  rval = read_coordinate( "t", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 't' variable" );
186  }
187  else
188  {
189  // If expected time variable does not exist, set dummy values to tVals
190  for( int t = 0; t < nTimeSteps; t++ )
191  tVals.push_back( (double)t );
192  }
193  }
194 
195  // For each variable, determine the entity location type and number of levels
196  std::map< std::string, ReadNC::VarData >::iterator mit;
197  for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
198  {
199  ReadNC::VarData& vd = ( *mit ).second;
200 
201  // Default entLoc is ENTLOCSET
202  if( std::find( vd.varDims.begin(), vd.varDims.end(), tDim ) != vd.varDims.end() )
203  {
204  if( std::find( vd.varDims.begin(), vd.varDims.end(), vDim ) != vd.varDims.end() )
206  }
207 
208  // Default numLev is 0
209  if( std::find( vd.varDims.begin(), vd.varDims.end(), levDim ) != vd.varDims.end() ) vd.numLev = nLevels;
210  }
211 
212  // Hack: create dummy variables for dimensions (like ncol) with no corresponding coordinate
213  // variables
214  rval = create_dummy_variables();MB_CHK_SET_ERR( rval, "Failed to create dummy variables" );
215 
216  return MB_SUCCESS;
217 }
218 
219 // When noMesh option is used on this read, the old ReadNC class instance for last read can get out
220 // of scope (and deleted). The old instance initialized localGidVerts properly when the mesh was
221 // created, but it is now lost. The new instance (will not create the mesh with noMesh option) has
222 // to restore it based on the existing mesh from last read
224 {
225  Interface*& mbImpl = _readNC->mbImpl;
226  Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
227  bool& noMesh = _readNC->noMesh;
228 
229  if( noMesh && localGidVerts.empty() )
230  {
231  // We need to populate localGidVerts range with the gids of vertices from current file set
232  // localGidVerts is important in reading the variable data into the nodes
233  // Also, for our purposes, localGidVerts is truly the GLOBAL_ID tag data, not other
234  // file_id tags that could get passed around in other scenarios for parallel reading
235 
236  // Get all vertices from current file set (it is the input set in no_mesh scenario)
237  Range local_verts;
238  ErrorCode rval = mbImpl->get_entities_by_dimension( _fileSet, 0, local_verts );MB_CHK_SET_ERR( rval, "Trouble getting local vertices in current file set" );
239 
240  if( !local_verts.empty() )
241  {
242  std::vector< int > gids( local_verts.size() );
243 
244  // !IMPORTANT : this has to be the GLOBAL_ID tag
245  rval = mbImpl->tag_get_data( mGlobalIdTag, local_verts, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting local gid values of vertices" );
246 
247  // Restore localGidVerts
248  std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidVerts ) );
250  }
251  }
252 
253  return MB_SUCCESS;
254 }
255 
257 {
258  Interface*& mbImpl = _readNC->mbImpl;
259  std::string& fileName = _readNC->fileName;
260  Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
261  const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
262  DebugOutput& dbgOut = _readNC->dbgOut;
263  bool& spectralMesh = _readNC->spectralMesh;
264  int& gatherSetRank = _readNC->gatherSetRank;
265  int& trivialPartitionShift = _readNC->trivialPartitionShift;
266 
267  int rank = 0;
268  int procs = 1;
269 #ifdef MOAB_HAVE_MPI
270  bool& isParallel = _readNC->isParallel;
271  if( isParallel )
272  {
273  ParallelComm*& myPcomm = _readNC->myPcomm;
274  rank = myPcomm->proc_config().proc_rank();
275  procs = myPcomm->proc_config().proc_size();
276  }
277 #endif
278 
279  ErrorCode rval;
280  int success = 0;
281 
282  // Need to get/read connectivity data before creating elements
283  std::string conn_fname;
284 
285  if( isConnFile )
286  {
287  // Connectivity file has already been read
289  }
290  else
291  {
292  // Try to open the connectivity file through CONN option, if used
293  rval = _opts.get_str_option( "CONN", conn_fname );
294  if( MB_SUCCESS != rval )
295  {
296  // Default convention for reading HOMME is a file HommeMapping.nc in same dir as data
297  // file
298  conn_fname = std::string( fileName );
299  size_t idx = conn_fname.find_last_of( "/" );
300  if( idx != std::string::npos )
301  conn_fname = conn_fname.substr( 0, idx ).append( "/HommeMapping.nc" );
302  else
303  conn_fname = "HommeMapping.nc";
304  }
305 #ifdef MOAB_HAVE_PNETCDF
306 #ifdef MOAB_HAVE_MPI
307  if( isParallel )
308  {
309  ParallelComm*& myPcomm = _readNC->myPcomm;
310  success =
311  NCFUNC( open )( myPcomm->proc_config().proc_comm(), conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId );
312  }
313  else
314  success = NCFUNC( open )( MPI_COMM_SELF, conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId );
315 #endif
316 #else
317  success = NCFUNC( open )( conn_fname.c_str(), 0, &connectId );
318 #endif
319  if( success ) MB_SET_ERR( MB_FAILURE, "Failed on open" );
320  }
321 
322  std::vector< std::string > conn_names;
323  std::vector< int > conn_vals;
324  rval = _readNC->get_dimensions( connectId, conn_names, conn_vals );MB_CHK_SET_ERR( rval, "Failed to get dimensions for connectivity" );
325 
326  // Read connectivity into temporary variable
327  int num_fine_quads = 0;
328  int num_coarse_quads = 0;
329  int start_idx = 0;
330  std::vector< std::string >::iterator vit;
331  int idx = 0;
332  if( ( vit = std::find( conn_names.begin(), conn_names.end(), "ncells" ) ) != conn_names.end() )
333  idx = vit - conn_names.begin();
334  else if( ( vit = std::find( conn_names.begin(), conn_names.end(), "ncenters" ) ) != conn_names.end() )
335  idx = vit - conn_names.begin();
336  else
337  {
338  MB_SET_ERR( MB_FAILURE, "Failed to get number of quads" );
339  }
340  int num_quads = conn_vals[idx];
341  if( !isConnFile && num_quads != nCells )
342  {
343  dbgOut.tprintf( 1,
344  "Warning: number of quads from %s and cells from %s are inconsistent; "
345  "num_quads = %d, nCells = %d.\n",
346  conn_fname.c_str(), fileName.c_str(), num_quads, nCells );
347  }
348 
349  // Get the connectivity into tmp_conn2 and permute into tmp_conn
350  int cornerVarId;
351  success = NCFUNC( inq_varid )( connectId, "element_corners", &cornerVarId );
352  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of 'element_corners'" );
353  NCDF_SIZE tmp_starts[2] = { 0, 0 };
354  NCDF_SIZE tmp_counts[2] = { 4, static_cast< NCDF_SIZE >( num_quads ) };
355  std::vector< int > tmp_conn( 4 * num_quads ), tmp_conn2( 4 * num_quads );
356  success = NCFUNCAG( _vara_int )( connectId, cornerVarId, tmp_starts, tmp_counts, &tmp_conn2[0] );
357  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get temporary connectivity" );
358  if( isConnFile )
359  {
360  // This data/connectivity file will be closed later in ReadNC::load_file()
361  }
362  else
363  {
364  success = NCFUNC( close )( connectId );
365  if( success ) MB_SET_ERR( MB_FAILURE, "Failed on close" );
366  }
367  // Permute the connectivity
368  for( int i = 0; i < num_quads; i++ )
369  {
370  tmp_conn[4 * i] = tmp_conn2[i];
371  tmp_conn[4 * i + 1] = tmp_conn2[i + 1 * num_quads];
372  tmp_conn[4 * i + 2] = tmp_conn2[i + 2 * num_quads];
373  tmp_conn[4 * i + 3] = tmp_conn2[i + 3 * num_quads];
374  }
375 
376  // Need to know whether we'll be creating gather mesh later, to make sure
377  // we allocate enough space in one shot
378  bool create_gathers = false;
379  if( rank == gatherSetRank ) create_gathers = true;
380 
381  // Shift rank to obtain a rotated trivial partition
382  int shifted_rank = rank;
383  if( procs >= 2 && trivialPartitionShift > 0 ) shifted_rank = ( rank + trivialPartitionShift ) % procs;
384 
385  // Compute the number of local quads, accounting for coarse or fine representation
386  // spectral_unit is the # fine quads per coarse quad, or spectralOrder^2
387  int spectral_unit = ( spectralMesh ? _spectralOrder * _spectralOrder : 1 );
388  // num_coarse_quads is the number of quads instantiated in MOAB; if !spectralMesh,
389  // num_coarse_quads = num_fine_quads
390  num_coarse_quads = int( std::floor( 1.0 * num_quads / ( spectral_unit * procs ) ) );
391  // start_idx is the starting index in the HommeMapping connectivity list for this proc, before
392  // converting to coarse quad representation
393  start_idx = 4 * shifted_rank * num_coarse_quads * spectral_unit;
394  // iextra = # coarse quads extra after equal split over procs
395  int iextra = num_quads % ( procs * spectral_unit );
396  if( shifted_rank < iextra ) num_coarse_quads++;
397  start_idx += 4 * spectral_unit * std::min( shifted_rank, iextra );
398  // num_fine_quads is the number of quads in the connectivity list in HommeMapping file assigned
399  // to this proc
400  num_fine_quads = spectral_unit * num_coarse_quads;
401 
402  // Now create num_coarse_quads
403  EntityHandle* conn_arr;
404  EntityHandle start_vertex;
405  Range tmp_range;
406 
407  // Read connectivity into that space
408  EntityHandle* sv_ptr = NULL;
409  EntityHandle start_quad;
410  SpectralMeshTool smt( mbImpl, _spectralOrder );
411  if( !spectralMesh )
412  {
413  rval = _readNC->readMeshIface->get_element_connect( num_coarse_quads, 4, MBQUAD, 0, start_quad, conn_arr,
414  // Might have to create gather mesh later
415  ( create_gathers ? num_coarse_quads + num_quads
416  : num_coarse_quads ) );MB_CHK_SET_ERR( rval, "Failed to create local quads" );
417  tmp_range.insert( start_quad, start_quad + num_coarse_quads - 1 );
418  int* tmp_conn_end = ( &tmp_conn[start_idx + 4 * num_fine_quads - 1] ) + 1;
419  std::copy( &tmp_conn[start_idx], tmp_conn_end, conn_arr );
420  std::copy( conn_arr, conn_arr + 4 * num_fine_quads, range_inserter( localGidVerts ) );
421  }
422  else
423  {
424  rval = smt.create_spectral_elems( &tmp_conn[0], num_fine_quads, 2, tmp_range, start_idx, &localGidVerts );MB_CHK_SET_ERR( rval, "Failed to create spectral elements" );
425  int count, v_per_e;
426  rval = mbImpl->connect_iterate( tmp_range.begin(), tmp_range.end(), conn_arr, v_per_e, count );MB_CHK_SET_ERR( rval, "Failed to get connectivity of spectral elements" );
427  rval = mbImpl->tag_iterate( smt.spectral_vertices_tag( true ), tmp_range.begin(), tmp_range.end(), count,
428  (void*&)sv_ptr );MB_CHK_SET_ERR( rval, "Failed to get fine connectivity of spectral elements" );
429  }
430 
431  // Create vertices
433  std::vector< double* > arrays;
434  rval = _readNC->readMeshIface->get_node_coords( 3, nLocalVertices, 0, start_vertex, arrays,
435  // Might have to create gather mesh later
436  ( create_gathers ? nLocalVertices + nVertices : nLocalVertices ) );MB_CHK_SET_ERR( rval, "Failed to create local vertices" );
437 
438  // Set vertex coordinates
439  Range::iterator rit;
440  double* xptr = arrays[0];
441  double* yptr = arrays[1];
442  double* zptr = arrays[2];
443  int i;
444  for( i = 0, rit = localGidVerts.begin(); i < nLocalVertices; i++, ++rit )
445  {
446  assert( *rit < xVertVals.size() + 1 );
447  xptr[i] = xVertVals[( *rit ) - 1]; // lon
448  yptr[i] = yVertVals[( *rit ) - 1]; // lat
449  }
450 
451  // Convert lon/lat/rad to x/y/z
452  const double pideg = acos( -1.0 ) / 180.0;
453  double rad = ( isConnFile ) ? 8000.0 : 8000.0 + levVals[0];
454  for( i = 0; i < nLocalVertices; i++ )
455  {
456  double cosphi = cos( pideg * yptr[i] );
457  double zmult = sin( pideg * yptr[i] );
458  double xmult = cosphi * cos( xptr[i] * pideg );
459  double ymult = cosphi * sin( xptr[i] * pideg );
460  xptr[i] = rad * xmult;
461  yptr[i] = rad * ymult;
462  zptr[i] = rad * zmult;
463  }
464 
465  // Get ptr to gid memory for vertices
466  Range vert_range( start_vertex, start_vertex + nLocalVertices - 1 );
467  void* data;
468  int count;
469  rval = mbImpl->tag_iterate( mGlobalIdTag, vert_range.begin(), vert_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate global id tag on local vertices" );
470  assert( count == nLocalVertices );
471  int* gid_data = (int*)data;
472  std::copy( localGidVerts.begin(), localGidVerts.end(), gid_data );
473 
474  // Duplicate global id data, which will be used to resolve sharing
475  if( mpFileIdTag )
476  {
477  rval = mbImpl->tag_iterate( *mpFileIdTag, vert_range.begin(), vert_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate file id tag on local vertices" );
478  assert( count == nLocalVertices );
479  int bytes_per_tag = 4;
480  rval = mbImpl->tag_get_bytes( *mpFileIdTag, bytes_per_tag );MB_CHK_SET_ERR( rval, "Can't get number of bytes for file id tag" );
481  if( 4 == bytes_per_tag )
482  {
483  gid_data = (int*)data;
484  std::copy( localGidVerts.begin(), localGidVerts.end(), gid_data );
485  }
486  else if( 8 == bytes_per_tag )
487  { // Should be a handle tag on 64 bit machine?
488  long* handle_tag_data = (long*)data;
489  std::copy( localGidVerts.begin(), localGidVerts.end(), handle_tag_data );
490  }
491  }
492 
493  // Create map from file ids to vertex handles, used later to set connectivity
494  std::map< EntityHandle, EntityHandle > vert_handles;
495  for( rit = localGidVerts.begin(), i = 0; rit != localGidVerts.end(); ++rit, i++ )
496  vert_handles[*rit] = start_vertex + i;
497 
498  // Compute proper handles in connectivity using offset
499  for( int q = 0; q < 4 * num_coarse_quads; q++ )
500  {
501  conn_arr[q] = vert_handles[conn_arr[q]];
502  assert( conn_arr[q] );
503  }
504  if( spectralMesh )
505  {
506  int verts_per_quad = ( _spectralOrder + 1 ) * ( _spectralOrder + 1 );
507  for( int q = 0; q < verts_per_quad * num_coarse_quads; q++ )
508  {
509  sv_ptr[q] = vert_handles[sv_ptr[q]];
510  assert( sv_ptr[q] );
511  }
512  }
513 
514  // Add new vertices and quads to current file set
515  faces.merge( tmp_range );
516  tmp_range.insert( start_vertex, start_vertex + nLocalVertices - 1 );
517  rval = mbImpl->add_entities( _fileSet, tmp_range );MB_CHK_SET_ERR( rval, "Failed to add new vertices and quads to current file set" );
518 
519  // Mark the set with the spectral order
520  Tag sporder;
521  rval = mbImpl->tag_get_handle( "SPECTRAL_ORDER", 1, MB_TYPE_INTEGER, sporder, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating SPECTRAL_ORDER tag" );
522  rval = mbImpl->tag_set_data( sporder, &_fileSet, 1, &_spectralOrder );MB_CHK_SET_ERR( rval, "Trouble setting data to SPECTRAL_ORDER tag" );
523 
524  if( create_gathers )
525  {
526  EntityHandle gather_set;
527  rval = _readNC->readMeshIface->create_gather_set( gather_set );MB_CHK_SET_ERR( rval, "Failed to create gather set" );
528 
529  // Create vertices
530  arrays.clear();
531  // Don't need to specify allocation number here, because we know enough verts were created
532  // before
533  rval = _readNC->readMeshIface->get_node_coords( 3, nVertices, 0, start_vertex, arrays );MB_CHK_SET_ERR( rval, "Failed to create gather set vertices" );
534 
535  xptr = arrays[0];
536  yptr = arrays[1];
537  zptr = arrays[2];
538  for( i = 0; i < nVertices; i++ )
539  {
540  double cosphi = cos( pideg * yVertVals[i] );
541  double zmult = sin( pideg * yVertVals[i] );
542  double xmult = cosphi * cos( xVertVals[i] * pideg );
543  double ymult = cosphi * sin( xVertVals[i] * pideg );
544  xptr[i] = rad * xmult;
545  yptr[i] = rad * ymult;
546  zptr[i] = rad * zmult;
547  }
548 
549  // Get ptr to gid memory for vertices
550  Range gather_set_verts_range( start_vertex, start_vertex + nVertices - 1 );
551  rval = mbImpl->tag_iterate( mGlobalIdTag, gather_set_verts_range.begin(), gather_set_verts_range.end(), count,
552  data );MB_CHK_SET_ERR( rval, "Failed to iterate global id tag on gather set vertices" );
553  assert( count == nVertices );
554  gid_data = (int*)data;
555  for( int j = 1; j <= nVertices; j++ )
556  gid_data[j - 1] = j;
557  // Set the file id tag too, it should be bigger something not interfering with global id
558  if( mpFileIdTag )
559  {
560  rval = mbImpl->tag_iterate( *mpFileIdTag, gather_set_verts_range.begin(), gather_set_verts_range.end(),
561  count, data );MB_CHK_SET_ERR( rval, "Failed to iterate file id tag on gather set vertices" );
562  assert( count == nVertices );
563  int bytes_per_tag = 4;
564  rval = mbImpl->tag_get_bytes( *mpFileIdTag, bytes_per_tag );MB_CHK_SET_ERR( rval, "Can't get number of bytes for file id tag" );
565  if( 4 == bytes_per_tag )
566  {
567  gid_data = (int*)data;
568  for( int j = 1; j <= nVertices; j++ )
569  gid_data[j - 1] = nVertices + j; // Bigger than global id tag
570  }
571  else if( 8 == bytes_per_tag )
572  { // Should be a handle tag on 64 bit machine?
573  long* handle_tag_data = (long*)data;
574  for( int j = 1; j <= nVertices; j++ )
575  handle_tag_data[j - 1] = nVertices + j; // Bigger than global id tag
576  }
577  }
578 
579  rval = mbImpl->add_entities( gather_set, gather_set_verts_range );MB_CHK_SET_ERR( rval, "Failed to add vertices to the gather set" );
580 
581  // Create quads
582  Range gather_set_quads_range;
583  // Don't need to specify allocation number here, because we know enough quads were created
584  // before
585  rval = _readNC->readMeshIface->get_element_connect( num_quads, 4, MBQUAD, 0, start_quad, conn_arr );MB_CHK_SET_ERR( rval, "Failed to create gather set quads" );
586  gather_set_quads_range.insert( start_quad, start_quad + num_quads - 1 );
587  int* tmp_conn_end = ( &tmp_conn[4 * num_quads - 1] ) + 1;
588  std::copy( &tmp_conn[0], tmp_conn_end, conn_arr );
589  for( i = 0; i != 4 * num_quads; i++ )
590  conn_arr[i] += start_vertex - 1; // Connectivity array is shifted by where the gather verts start
591  rval = mbImpl->add_entities( gather_set, gather_set_quads_range );MB_CHK_SET_ERR( rval, "Failed to add quads to the gather set" );
592  }
593 
594  return MB_SUCCESS;
595 }
596 
597 ErrorCode NCHelperHOMME::read_ucd_variables_to_nonset_allocate( std::vector< ReadNC::VarData >& vdatas,
598  std::vector< int >& tstep_nums )
599 {
600  Interface*& mbImpl = _readNC->mbImpl;
601  std::vector< int >& dimLens = _readNC->dimLens;
602  DebugOutput& dbgOut = _readNC->dbgOut;
603 
604  Range* range = NULL;
605 
606  // Get vertices
607  Range verts;
608  ErrorCode rval = mbImpl->get_entities_by_dimension( _fileSet, 0, verts );MB_CHK_SET_ERR( rval, "Trouble getting vertices in current file set" );
609  assert( "Should only have a single vertex subrange, since they were read in one shot" && verts.psize() == 1 );
610 
611  for( unsigned int i = 0; i < vdatas.size(); i++ )
612  {
613  // Support non-set variables with 3 dimensions like (time, lev, ncol)
614  assert( 3 == vdatas[i].varDims.size() );
615 
616  // For a non-set variable, time should be the first dimension
617  assert( tDim == vdatas[i].varDims[0] );
618 
619  // Set up readStarts and readCounts
620  vdatas[i].readStarts.resize( 3 );
621  vdatas[i].readCounts.resize( 3 );
622 
623  // First: time
624  vdatas[i].readStarts[0] = 0; // This value is timestep dependent, will be set later
625  vdatas[i].readCounts[0] = 1;
626 
627  // Next: lev
628  vdatas[i].readStarts[1] = 0;
629  vdatas[i].readCounts[1] = vdatas[i].numLev;
630 
631  // Finally: ncol
632  switch( vdatas[i].entLoc )
633  {
634  case ReadNC::ENTLOCVERT:
635  // Vertices
636  // Start from the first localGidVerts
637  // Actually, this will be reset later on in a loop
638  vdatas[i].readStarts[2] = localGidVerts[0] - 1;
639  vdatas[i].readCounts[2] = nLocalVertices;
640  range = &verts;
641  break;
642  default:
643  MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << vdatas[i].varName );
644  }
645 
646  // Get variable size
647  vdatas[i].sz = 1;
648  for( std::size_t idx = 0; idx != 3; idx++ )
649  vdatas[i].sz *= vdatas[i].readCounts[idx];
650 
651  for( unsigned int t = 0; t < tstep_nums.size(); t++ )
652  {
653  dbgOut.tprintf( 2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t] );
654 
655  if( tstep_nums[t] >= dimLens[tDim] )
656  {
657  MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number " << tstep_nums[t] );
658  }
659 
660  // Get the tag to read into
661  if( !vdatas[i].varTags[t] )
662  {
663  rval = get_tag_to_nonset( vdatas[i], tstep_nums[t], vdatas[i].varTags[t], vdatas[i].numLev );MB_CHK_SET_ERR( rval, "Trouble getting tag for variable " << vdatas[i].varName );
664  }
665 
666  // Get ptr to tag space
667  void* data;
668  int count;
669  rval = mbImpl->tag_iterate( vdatas[i].varTags[t], range->begin(), range->end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate tag for variable " << vdatas[i].varName );
670  assert( (unsigned)count == range->size() );
671  vdatas[i].varDatas[t] = data;
672  }
673  }
674 
675  return rval;
676 }
677 
678 #ifdef MOAB_HAVE_PNETCDF
679 ErrorCode NCHelperHOMME::read_ucd_variables_to_nonset_async( std::vector< ReadNC::VarData >& vdatas,
680  std::vector< int >& tstep_nums )
681 {
682  DebugOutput& dbgOut = _readNC->dbgOut;
683 
684  ErrorCode rval = read_ucd_variables_to_nonset_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read non-set variables" );
685 
686  // Finally, read into that space
687  int success;
688 
689  for( unsigned int i = 0; i < vdatas.size(); i++ )
690  {
691  std::size_t sz = vdatas[i].sz;
692 
693  // A typical supported variable: float T(time, lev, ncol)
694  // For tag values, need transpose (lev, ncol) to (ncol, lev)
695  size_t ni = vdatas[i].readCounts[2]; // ncol
696  size_t nj = 1; // Here we should just set nj to 1
697  size_t nk = vdatas[i].readCounts[1]; // lev
698 
699  for( unsigned int t = 0; t < tstep_nums.size(); t++ )
700  {
701  // We will synchronize all these reads with the other processors,
702  // so the wait will be inside this double loop; is it too much?
703  size_t nb_reads = localGidVerts.psize();
704  std::vector< int > requests( nb_reads ), statuss( nb_reads );
705  size_t idxReq = 0;
706 
707  // Tag data for this timestep
708  void* data = vdatas[i].varDatas[t];
709 
710  // Set readStart for each timestep along time dimension
711  vdatas[i].readStarts[0] = tstep_nums[t];
712 
713  switch( vdatas[i].varDataType )
714  {
715  case NC_FLOAT:
716  case NC_DOUBLE: {
717  // Read float as double
718  std::vector< double > tmpdoubledata( sz );
719 
720  // In the case of ucd mesh, and on multiple proc,
721  // we need to read as many times as subranges we have in the
722  // localGidVerts range;
723  // basically, we have to give a different point
724  // for data to start, for every subrange :(
725  size_t indexInDoubleArray = 0;
726  size_t ic = 0;
727  for( Range::pair_iterator pair_iter = localGidVerts.pair_begin();
728  pair_iter != localGidVerts.pair_end(); ++pair_iter, ic++ )
729  {
730  EntityHandle starth = pair_iter->first;
731  EntityHandle endh = pair_iter->second; // Inclusive
732  vdatas[i].readStarts[2] = (NCDF_SIZE)( starth - 1 );
733  vdatas[i].readCounts[2] = (NCDF_SIZE)( endh - starth + 1 );
734 
735  // Do a partial read, in each subrange
736  // Wait outside this loop
737  success =
738  NCFUNCREQG( _vara_double )( _fileId, vdatas[i].varId, &( vdatas[i].readStarts[0] ),
739  &( vdatas[i].readCounts[0] ),
740  &( tmpdoubledata[indexInDoubleArray] ), &requests[idxReq++] );
741  if( success )
742  MB_SET_ERR( MB_FAILURE,
743  "Failed to read double data in a loop for variable " << vdatas[i].varName );
744  // We need to increment the index in double array for the
745  // next subrange
746  indexInDoubleArray += ( endh - starth + 1 ) * 1 * vdatas[i].numLev;
747  }
748  assert( ic == localGidVerts.psize() );
749 
750  success = ncmpi_wait_all( _fileId, requests.size(), &requests[0], &statuss[0] );
751  if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
752 
753  if( vdatas[i].numLev > 1 )
754  // Transpose (lev, ncol) to (ncol, lev)
755  kji_to_jik_stride( ni, nj, nk, data, &tmpdoubledata[0], localGidVerts );
756  else
757  {
758  for( std::size_t idx = 0; idx != tmpdoubledata.size(); idx++ )
759  ( (double*)data )[idx] = tmpdoubledata[idx];
760  }
761 
762  break;
763  }
764  default:
765  MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
766  }
767  }
768  }
769 
770  // Debug output, if requested
771  if( 1 == dbgOut.get_verbosity() )
772  {
773  dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
774  for( unsigned int i = 1; i < vdatas.size(); i++ )
775  dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
776  dbgOut.tprintf( 1, "\n" );
777  }
778 
779  return rval;
780 }
781 #else
782 ErrorCode NCHelperHOMME::read_ucd_variables_to_nonset( std::vector< ReadNC::VarData >& vdatas,
783  std::vector< int >& tstep_nums )
784 {
785  DebugOutput& dbgOut = _readNC->dbgOut;
786 
787  ErrorCode rval = read_ucd_variables_to_nonset_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read non-set variables" );
788 
789  // Finally, read into that space
790  int success;
791  for( unsigned int i = 0; i < vdatas.size(); i++ )
792  {
793  std::size_t sz = vdatas[i].sz;
794 
795  // A typical supported variable: float T(time, lev, ncol)
796  // For tag values, need transpose (lev, ncol) to (ncol, lev)
797  size_t ni = vdatas[i].readCounts[2]; // ncol
798  size_t nj = 1; // Here we should just set nj to 1
799  size_t nk = vdatas[i].readCounts[1]; // lev
800 
801  for( unsigned int t = 0; t < tstep_nums.size(); t++ )
802  {
803  // Tag data for this timestep
804  void* data = vdatas[i].varDatas[t];
805 
806  // Set readStart for each timestep along time dimension
807  vdatas[i].readStarts[0] = tstep_nums[t];
808 
809  switch( vdatas[i].varDataType )
810  {
811  case NC_FLOAT:
812  case NC_DOUBLE: {
813  // Read float as double
814  std::vector< double > tmpdoubledata( sz );
815 
816  // In the case of ucd mesh, and on multiple proc,
817  // we need to read as many times as subranges we have in the
818  // localGidVerts range;
819  // basically, we have to give a different point
820  // for data to start, for every subrange :(
821  size_t indexInDoubleArray = 0;
822  size_t ic = 0;
823  for( Range::pair_iterator pair_iter = localGidVerts.pair_begin();
824  pair_iter != localGidVerts.pair_end(); ++pair_iter, ic++ )
825  {
826  EntityHandle starth = pair_iter->first;
827  EntityHandle endh = pair_iter->second; // Inclusive
828  vdatas[i].readStarts[2] = (NCDF_SIZE)( starth - 1 );
829  vdatas[i].readCounts[2] = (NCDF_SIZE)( endh - starth + 1 );
830 
831  success = NCFUNCAG( _vara_double )( _fileId, vdatas[i].varId, &( vdatas[i].readStarts[0] ),
832  &( vdatas[i].readCounts[0] ),
833  &( tmpdoubledata[indexInDoubleArray] ) );
834  if( success )
835  MB_SET_ERR( MB_FAILURE,
836  "Failed to read double data in a loop for variable " << vdatas[i].varName );
837  // We need to increment the index in double array for the
838  // next subrange
839  indexInDoubleArray += ( endh - starth + 1 ) * 1 * vdatas[i].numLev;
840  }
841  assert( ic == localGidVerts.psize() );
842 
843  if( vdatas[i].numLev > 1 )
844  // Transpose (lev, ncol) to (ncol, lev)
845  kji_to_jik_stride( ni, nj, nk, data, &tmpdoubledata[0], localGidVerts );
846  else
847  {
848  for( std::size_t idx = 0; idx != tmpdoubledata.size(); idx++ )
849  ( (double*)data )[idx] = tmpdoubledata[idx];
850  }
851 
852  break;
853  }
854  default:
855  MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
856  }
857  }
858  }
859 
860  // Debug output, if requested
861  if( 1 == dbgOut.get_verbosity() )
862  {
863  dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
864  for( unsigned int i = 1; i < vdatas.size(); i++ )
865  dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
866  dbgOut.tprintf( 1, "\n" );
867  }
868 
869  return rval;
870 }
871 #endif
872 
873 } // namespace moab