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