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