 |
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 "NCHelperFV.hpp"
2 #include "moab/FileOptions.hpp"
3
4 #include <cmath>
5 #include <sstream>
6
7 namespace moab
8 {
9
10 bool NCHelperFV::can_read_file( ReadNC* readNC, int fileId )
11 {
12 std::vector< std::string >& dimNames = readNC->dimNames;
13
14 // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it should be the FV grid
15 if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "lon" ) ) != dimNames.end() ) &&
16 ( std::find( dimNames.begin(), dimNames.end(), std::string( "lat" ) ) != dimNames.end() ) &&
17 ( std::find( dimNames.begin(), dimNames.end(), std::string( "slon" ) ) != dimNames.end() ) &&
18 ( std::find( dimNames.begin(), dimNames.end(), std::string( "slat" ) ) != dimNames.end() ) )
19 {
20 // Make sure it is CAM grid
21 std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "source" );
22 if( attIt == readNC->globalAtts.end() ) return false;
23 unsigned int sz = attIt->second.attLen;
24 std::string att_data;
25 att_data.resize( sz + 1 );
26 att_data[sz] = '\000';
27 int success =
28 NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
29 if( success ) return false;
30 if( att_data.find( "CAM" ) == std::string::npos ) return false;
31
32 return true;
33 }
34
35 return false;
36 }
37
38 ErrorCode NCHelperFV::init_mesh_vals()
39 {
40 Interface*& mbImpl = _readNC->mbImpl;
41 std::vector< std::string >& dimNames = _readNC->dimNames;
42 std::vector< int >& dimLens = _readNC->dimLens;
43 std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
44 DebugOutput& dbgOut = _readNC->dbgOut;
45 bool& isParallel = _readNC->isParallel;
46 int& partMethod = _readNC->partMethod;
47 ScdParData& parData = _readNC->parData;
48
49 // Look for names of i/j dimensions
50 // First i
51 std::vector< std::string >::iterator vit;
52 unsigned int idx;
53 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "slon" ) ) != dimNames.end() )
54 idx = vit - dimNames.begin();
55 else
56 {
57 MB_SET_ERR( MB_FAILURE, "Couldn't find 'slon' variable" );
58 }
59 iDim = idx;
60 gDims[0] = 0;
61 gDims[3] = dimLens[idx] - 1;
62
63 // Then j
64 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "slat" ) ) != dimNames.end() )
65 idx = vit - dimNames.begin();
66 else
67 {
68 MB_SET_ERR( MB_FAILURE, "Couldn't find 'slat' variable" );
69 }
70 jDim = idx;
71 gDims[1] = 0;
72 gDims[4] = dimLens[idx] - 1 + 2; // Add 2 for the pole points
73
74 // Look for names of center i/j dimensions
75 // First i
76 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lon" ) ) != dimNames.end() )
77 idx = vit - dimNames.begin();
78 else
79 {
80 MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' variable" );
81 }
82 iCDim = idx;
83 gCDims[0] = 0;
84 gCDims[3] = dimLens[idx] - 1;
85
86 // Check i periodicity and set globallyPeriodic[0]
87 std::vector< double > til_vals( 2 );
88 ErrorCode rval = read_coordinate( "lon", gCDims[3] - 1, gCDims[3], til_vals );MB_CHK_SET_ERR( rval, "Trouble reading 'lon' variable" );
89 if( std::fabs( 2 * til_vals[1] - til_vals[0] - 360 ) < 0.001 ) globallyPeriodic[0] = 1;
90 if( globallyPeriodic[0] )
91 assert( "Number of vertices and edges should be same" && gDims[3] == gCDims[3] );
92 else
93 assert( "Number of vertices should equal to number of edges plus one" && gDims[3] == gCDims[3] + 1 );
94
95 // Then j
96 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lat" ) ) != dimNames.end() )
97 idx = vit - dimNames.begin();
98 else
99 {
100 MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' dimension" );
101 }
102 jCDim = idx;
103 gCDims[1] = 0;
104 gCDims[4] = dimLens[idx] - 1;
105
106 // For FV models, will always be non-periodic in j
107 assert( gDims[4] == gCDims[4] + 1 );
108
109 // Try a truly 2D mesh
110 gDims[2] = -1;
111 gDims[5] = -1;
112
113 // Look for time dimension
114 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() )
115 idx = vit - dimNames.begin();
116 else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "t" ) ) != dimNames.end() )
117 idx = vit - dimNames.begin();
118 else
119 {
120 MB_SET_ERR( MB_FAILURE, "Couldn't find 'time' or 't' dimension" );
121 }
122 tDim = idx;
123 nTimeSteps = dimLens[idx];
124
125 // Get number of levels
126 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lev" ) ) != dimNames.end() )
127 idx = vit - dimNames.begin();
128 else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ilev" ) ) != dimNames.end() )
129 idx = vit - dimNames.begin();
130 else
131 {
132 MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension" );
133 }
134 levDim = idx;
135 nLevels = dimLens[idx];
136
137 // Parse options to get subset
138 int rank = 0, procs = 1;
139 #ifdef MOAB_HAVE_MPI
140 if( isParallel )
141 {
142 ParallelComm*& myPcomm = _readNC->myPcomm;
143 rank = myPcomm->proc_config().proc_rank();
144 procs = myPcomm->proc_config().proc_size();
145 }
146 #endif
147 if( procs > 1 )
148 {
149 for( int i = 0; i < 6; i++ )
150 parData.gDims[i] = gDims[i];
151 for( int i = 0; i < 3; i++ )
152 parData.gPeriodic[i] = globallyPeriodic[i];
153 parData.partMethod = partMethod;
154 int pdims[3];
155
156 rval = ScdInterface::compute_partition( procs, rank, parData, lDims, locallyPeriodic, pdims );MB_CHK_ERR( rval );
157 for( int i = 0; i < 3; i++ )
158 parData.pDims[i] = pdims[i];
159
160 dbgOut.tprintf( 1, "Partition: %dx%d (out of %dx%d)\n", lDims[3] - lDims[0] + 1, lDims[4] - lDims[1] + 1,
161 gDims[3] - gDims[0] + 1, gDims[4] - gDims[1] + 1 );
162 if( 0 == rank )
163 dbgOut.tprintf( 1, "Contiguous chunks of size %d bytes.\n",
164 8 * ( lDims[3] - lDims[0] + 1 ) * ( lDims[4] - lDims[1] + 1 ) );
165 }
166 else
167 {
168 for( int i = 0; i < 6; i++ )
169 lDims[i] = gDims[i];
170 locallyPeriodic[0] = globallyPeriodic[0];
171 }
172
173 _opts.get_int_option( "IMIN", lDims[0] );
174 _opts.get_int_option( "IMAX", lDims[3] );
175 _opts.get_int_option( "JMIN", lDims[1] );
176 _opts.get_int_option( "JMAX", lDims[4] );
177
178 // Now get actual coordinate values for vertices and cell centers
179 lCDims[0] = lDims[0];
180 if( locallyPeriodic[0] )
181 // If locally periodic, doesn't matter what global periodicity is, # vertex coords = # elem
182 // coords
183 lCDims[3] = lDims[3];
184 else
185 lCDims[3] = lDims[3] - 1;
186
187 // For FV models, will always be non-periodic in j
188 lCDims[1] = lDims[1];
189 lCDims[4] = lDims[4] - 1;
190
191 // Resize vectors to store values later
192 if( -1 != lDims[0] ) ilVals.resize( lDims[3] - lDims[0] + 1 );
193 if( -1 != lCDims[0] ) ilCVals.resize( lCDims[3] - lCDims[0] + 1 );
194 if( -1 != lDims[1] ) jlVals.resize( lDims[4] - lDims[1] + 1 );
195 if( -1 != lCDims[1] ) jlCVals.resize( lCDims[4] - lCDims[1] + 1 );
196 if( nTimeSteps > 0 ) tVals.resize( nTimeSteps );
197
198 // Now read coord values
199 std::map< std::string, ReadNC::VarData >::iterator vmit;
200 if( -1 != lCDims[0] )
201 {
202 if( ( vmit = varInfo.find( "lon" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
203 {
204 rval = read_coordinate( "lon", lCDims[0], lCDims[3], ilCVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lon' variable" );
205 }
206 else
207 {
208 MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' variable" );
209 }
210 }
211
212 if( -1 != lCDims[1] )
213 {
214 if( ( vmit = varInfo.find( "lat" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
215 {
216 rval = read_coordinate( "lat", lCDims[1], lCDims[4], jlCVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lat' variable" );
217 }
218 else
219 {
220 MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' variable" );
221 }
222 }
223
224 if( -1 != lDims[0] )
225 {
226 if( ( vmit = varInfo.find( "slon" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
227 {
228 // Last column
229 if( !locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3] )
230 {
231 assert( lDims[3] == gDims[3] + 1 );
232 std::vector< double > dummyVar( lDims[3] - lDims[0] );
233 rval = read_coordinate( "slon", lDims[0], lDims[3] - 1, dummyVar );
234 double dif = dummyVar[1] - dummyVar[0];
235 std::size_t i;
236 for( i = 0; i != dummyVar.size(); i++ )
237 ilVals[i] = dummyVar[i];
238 ilVals[i] = ilVals[i - 1] + dif;
239 }
240 else
241 {
242 rval = read_coordinate( "slon", lDims[0], lDims[3], ilVals );MB_CHK_SET_ERR( rval, "Trouble reading 'slon' variable" );
243 }
244 }
245 else
246 {
247 MB_SET_ERR( MB_FAILURE, "Couldn't find 'slon' variable" );
248 }
249 }
250
251 if( -1 != lDims[1] )
252 {
253 if( ( vmit = varInfo.find( "slat" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
254 {
255 if( !isParallel || ( ( gDims[4] - gDims[1] ) == ( lDims[4] - lDims[1] ) ) )
256 {
257 std::vector< double > dummyVar( lDims[4] - lDims[1] - 1 );
258 rval = read_coordinate( "slat", lDims[1], lDims[4] - 2, dummyVar );MB_CHK_SET_ERR( rval, "Trouble reading 'slat' variable" );
259 // Copy the correct piece
260 jlVals[0] = -90.0;
261 std::size_t i = 0;
262 for( i = 1; i != dummyVar.size() + 1; i++ )
263 jlVals[i] = dummyVar[i - 1];
264 jlVals[i] = 90.0; // Using value of i after loop exits.
265 }
266 else
267 {
268 // If this is the first row
269 // Need to read one less then available and read it into a dummy var
270 if( lDims[1] == gDims[1] )
271 {
272 std::vector< double > dummyVar( lDims[4] - lDims[1] );
273 rval = read_coordinate( "slat", lDims[1], lDims[4] - 1, dummyVar );MB_CHK_SET_ERR( rval, "Trouble reading 'slat' variable" );
274 // Copy the correct piece
275 jlVals[0] = -90.0;
276 for( int i = 1; i < lDims[4] + 1; i++ )
277 jlVals[i] = dummyVar[i - 1];
278 }
279 // Or if it's the last row
280 else if( lDims[4] == gDims[4] )
281 {
282 std::vector< double > dummyVar( lDims[4] - lDims[1] );
283 rval = read_coordinate( "slat", lDims[1] - 1, lDims[4] - 2, dummyVar );MB_CHK_SET_ERR( rval, "Trouble reading 'slat' variable" );
284 // Copy the correct piece
285 std::size_t i = 0;
286 for( i = 0; i != dummyVar.size(); i++ )
287 jlVals[i] = dummyVar[i];
288 jlVals[i] = 90.0; // Using value of i after loop exits.
289 }
290 // It's in the middle
291 else
292 {
293 rval = read_coordinate( "slat", lDims[1] - 1, lDims[4] - 1, jlVals );MB_CHK_SET_ERR( rval, "Trouble reading 'slat' variable" );
294 }
295 }
296 }
297 else
298 {
299 MB_SET_ERR( MB_FAILURE, "Couldn't find 'slat' variable" );
300 }
301 }
302
303 // Store time coordinate values in tVals
304 if( nTimeSteps > 0 )
305 {
306 if( ( vmit = varInfo.find( "time" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
307 {
308 rval = read_coordinate( "time", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 'time' variable" );
309 }
310 else if( ( vmit = varInfo.find( "t" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
311 {
312 rval = read_coordinate( "t", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 't' variable" );
313 }
314 else
315 {
316 // If expected time variable is not available, set dummy time coordinate values to tVals
317 for( int t = 0; t < nTimeSteps; t++ )
318 tVals.push_back( (double)t );
319 }
320 }
321
322 dbgOut.tprintf( 1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4] );
323 dbgOut.tprintf( 1, "%d elements, %d vertices\n", ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ),
324 ( lDims[3] - lDims[0] + 1 ) * ( lDims[4] - lDims[1] + 1 ) );
325
326 // For each variable, determine the entity location type and number of levels
327 std::map< std::string, ReadNC::VarData >::iterator mit;
328 for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
329 {
330 ReadNC::VarData& vd = ( *mit ).second;
331
332 // Default entLoc is ENTLOCSET
333 if( std::find( vd.varDims.begin(), vd.varDims.end(), tDim ) != vd.varDims.end() )
334 {
335 if( ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) &&
336 ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) )
337 vd.entLoc = ReadNC::ENTLOCFACE;
338 else if( ( std::find( vd.varDims.begin(), vd.varDims.end(), jDim ) != vd.varDims.end() ) &&
339 ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) )
340 vd.entLoc = ReadNC::ENTLOCNSEDGE;
341 else if( ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) &&
342 ( std::find( vd.varDims.begin(), vd.varDims.end(), iDim ) != vd.varDims.end() ) )
343 vd.entLoc = ReadNC::ENTLOCEWEDGE;
344 }
345
346 // Default numLev is 0
347 if( std::find( vd.varDims.begin(), vd.varDims.end(), levDim ) != vd.varDims.end() ) vd.numLev = nLevels;
348 }
349
350 std::vector< std::string > ijdimNames( 4 );
351 ijdimNames[0] = "__slon";
352 ijdimNames[1] = "__slat";
353 ijdimNames[2] = "__lon";
354 ijdimNames[3] = "__lat";
355
356 std::string tag_name;
357 Tag tagh;
358
359 // __<dim_name>_LOC_MINMAX (for slon, slat, lon and lat)
360 for( unsigned int i = 0; i != ijdimNames.size(); i++ )
361 {
362 std::vector< int > val( 2, 0 );
363 if( ijdimNames[i] == "__slon" )
364 {
365 val[0] = lDims[0];
366 val[1] = lDims[3];
367 }
368 else if( ijdimNames[i] == "__slat" )
369 {
370 val[0] = lDims[1];
371 val[1] = lDims[4];
372 }
373 else if( ijdimNames[i] == "__lon" )
374 {
375 val[0] = lCDims[0];
376 val[1] = lCDims[3];
377 }
378 else if( ijdimNames[i] == "__lat" )
379 {
380 val[0] = lCDims[1];
381 val[1] = lCDims[4];
382 }
383 std::stringstream ss_tag_name;
384 ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
385 tag_name = ss_tag_name.str();
386 rval = mbImpl->tag_get_handle( tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
387 rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
388 if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
389 }
390
391 // __<dim_name>_LOC_VALS (for slon, slat, lon and lat)
392 // Assume all have the same data type as lon (expected type is float or double)
393 switch( varInfo["lon"].varDataType )
394 {
395 case NC_FLOAT:
396 case NC_DOUBLE:
397 break;
398 default:
399 MB_SET_ERR( MB_FAILURE, "Unexpected variable data type for 'lon'" );
400 }
401
402 for( unsigned int i = 0; i != ijdimNames.size(); i++ )
403 {
404 void* val = NULL;
405 int val_len = 0;
406 if( ijdimNames[i] == "__slon" )
407 {
408 val = &ilVals[0];
409 val_len = ilVals.size();
410 }
411 else if( ijdimNames[i] == "__slat" )
412 {
413 val = &jlVals[0];
414 val_len = jlVals.size();
415 }
416 else if( ijdimNames[i] == "__lon" )
417 {
418 val = &ilCVals[0];
419 val_len = ilCVals.size();
420 }
421 else if( ijdimNames[i] == "__lat" )
422 {
423 val = &jlCVals[0];
424 val_len = jlCVals.size();
425 }
426
427 std::stringstream ss_tag_name;
428 ss_tag_name << ijdimNames[i] << "_LOC_VALS";
429 tag_name = ss_tag_name.str();
430 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_DOUBLE, tagh,
431 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
432 rval = mbImpl->tag_set_by_ptr( tagh, &_fileSet, 1, &val, &val_len );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
433 if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
434 }
435
436 // __<dim_name>_GLOBAL_MINMAX (for slon, slat, lon and lat)
437 for( unsigned int i = 0; i != ijdimNames.size(); i++ )
438 {
439 std::vector< int > val( 2, 0 );
440 if( ijdimNames[i] == "__slon" )
441 {
442 val[0] = gDims[0];
443 val[1] = gDims[3];
444 }
445 else if( ijdimNames[i] == "__slat" )
446 {
447 val[0] = gDims[1];
448 val[1] = gDims[4];
449 }
450 else if( ijdimNames[i] == "__lon" )
451 {
452 val[0] = gCDims[0];
453 val[1] = gCDims[3];
454 }
455 else if( ijdimNames[i] == "__lat" )
456 {
457 val[0] = gCDims[1];
458 val[1] = gCDims[4];
459 }
460 std::stringstream ss_tag_name;
461 ss_tag_name << ijdimNames[i] << "_GLOBAL_MINMAX";
462 tag_name = ss_tag_name.str();
463 rval = mbImpl->tag_get_handle( tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
464 rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
465 if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
466 }
467
468 // Hack: create dummy variables, if needed, for dimensions with no corresponding coordinate
469 // variables
470 rval = create_dummy_variables();MB_CHK_SET_ERR( rval, "Failed to create dummy variables" );
471
472 return MB_SUCCESS;
473 }
474
475 } // namespace moab
1.9.1.