 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 #include "NCHelperHOMME.hpp"
2 #include "moab/ReadUtilIface.hpp"
3 #include "moab/FileOptions.hpp"
4 #include "moab/SpectralMeshTool.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
64 ErrorCode NCHelperHOMME::init_mesh_vals()
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() )
205 vd.entLoc = ReadNC::ENTLOCVERT;
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
223 ErrorCode NCHelperHOMME::check_existing_mesh()
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 ) );
249 nLocalVertices = localGidVerts.size();
250 }
251 }
252
253 return MB_SUCCESS;
254 }
255
256 ErrorCode NCHelperHOMME::create_mesh( Range& faces )
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
288 connectId = _readNC->fileId;
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
432 nLocalVertices = localGidVerts.size();
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
1.9.1.