|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
00001 #include "NCHelperFV.hpp" 00002 #include "moab/FileOptions.hpp" 00003 00004 #include <cmath> 00005 #include <sstream> 00006 00007 namespace moab 00008 { 00009 00010 bool NCHelperFV::can_read_file( ReadNC* readNC, int fileId ) 00011 { 00012 std::vector< std::string >& dimNames = readNC->dimNames; 00013 00014 // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it should be the FV grid 00015 if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "lon" ) ) != dimNames.end() ) && 00016 ( std::find( dimNames.begin(), dimNames.end(), std::string( "lat" ) ) != dimNames.end() ) && 00017 ( std::find( dimNames.begin(), dimNames.end(), std::string( "slon" ) ) != dimNames.end() ) && 00018 ( std::find( dimNames.begin(), dimNames.end(), std::string( "slat" ) ) != dimNames.end() ) ) 00019 { 00020 // Make sure it is CAM grid 00021 std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "source" ); 00022 if( attIt == readNC->globalAtts.end() ) return false; 00023 unsigned int sz = attIt->second.attLen; 00024 std::string att_data; 00025 att_data.resize( sz + 1 ); 00026 att_data[sz] = '\000'; 00027 int success = 00028 NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] ); 00029 if( success ) return false; 00030 if( att_data.find( "CAM" ) == std::string::npos ) return false; 00031 00032 return true; 00033 } 00034 00035 return false; 00036 } 00037 00038 ErrorCode NCHelperFV::init_mesh_vals() 00039 { 00040 Interface*& mbImpl = _readNC->mbImpl; 00041 std::vector< std::string >& dimNames = _readNC->dimNames; 00042 std::vector< int >& dimLens = _readNC->dimLens; 00043 std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo; 00044 DebugOutput& dbgOut = _readNC->dbgOut; 00045 bool& isParallel = _readNC->isParallel; 00046 int& partMethod = _readNC->partMethod; 00047 ScdParData& parData = _readNC->parData; 00048 00049 // Look for names of i/j dimensions 00050 // First i 00051 std::vector< std::string >::iterator vit; 00052 unsigned int idx; 00053 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "slon" ) ) != dimNames.end() ) 00054 idx = vit - dimNames.begin(); 00055 else 00056 { 00057 MB_SET_ERR( MB_FAILURE, "Couldn't find 'slon' variable" ); 00058 } 00059 iDim = idx; 00060 gDims[0] = 0; 00061 gDims[3] = dimLens[idx] - 1; 00062 00063 // Then j 00064 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "slat" ) ) != dimNames.end() ) 00065 idx = vit - dimNames.begin(); 00066 else 00067 { 00068 MB_SET_ERR( MB_FAILURE, "Couldn't find 'slat' variable" ); 00069 } 00070 jDim = idx; 00071 gDims[1] = 0; 00072 gDims[4] = dimLens[idx] - 1 + 2; // Add 2 for the pole points 00073 00074 // Look for names of center i/j dimensions 00075 // First i 00076 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lon" ) ) != dimNames.end() ) 00077 idx = vit - dimNames.begin(); 00078 else 00079 { 00080 MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' variable" ); 00081 } 00082 iCDim = idx; 00083 gCDims[0] = 0; 00084 gCDims[3] = dimLens[idx] - 1; 00085 00086 // Check i periodicity and set globallyPeriodic[0] 00087 std::vector< double > til_vals( 2 ); 00088 ErrorCode rval = read_coordinate( "lon", gCDims[3] - 1, gCDims[3], til_vals );MB_CHK_SET_ERR( rval, "Trouble reading 'lon' variable" ); 00089 if( std::fabs( 2 * til_vals[1] - til_vals[0] - 360 ) < 0.001 ) globallyPeriodic[0] = 1; 00090 if( globallyPeriodic[0] ) 00091 assert( "Number of vertices and edges should be same" && gDims[3] == gCDims[3] ); 00092 else 00093 assert( "Number of vertices should equal to number of edges plus one" && gDims[3] == gCDims[3] + 1 ); 00094 00095 // Then j 00096 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lat" ) ) != dimNames.end() ) 00097 idx = vit - dimNames.begin(); 00098 else 00099 { 00100 MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' dimension" ); 00101 } 00102 jCDim = idx; 00103 gCDims[1] = 0; 00104 gCDims[4] = dimLens[idx] - 1; 00105 00106 // For FV models, will always be non-periodic in j 00107 assert( gDims[4] == gCDims[4] + 1 ); 00108 00109 // Try a truly 2D mesh 00110 gDims[2] = -1; 00111 gDims[5] = -1; 00112 00113 // Look for time dimension 00114 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() ) 00115 idx = vit - dimNames.begin(); 00116 else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "t" ) ) != dimNames.end() ) 00117 idx = vit - dimNames.begin(); 00118 else 00119 { 00120 MB_SET_ERR( MB_FAILURE, "Couldn't find 'time' or 't' dimension" ); 00121 } 00122 tDim = idx; 00123 nTimeSteps = dimLens[idx]; 00124 00125 // Get number of levels 00126 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lev" ) ) != dimNames.end() ) 00127 idx = vit - dimNames.begin(); 00128 else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ilev" ) ) != dimNames.end() ) 00129 idx = vit - dimNames.begin(); 00130 else 00131 { 00132 MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension" ); 00133 } 00134 levDim = idx; 00135 nLevels = dimLens[idx]; 00136 00137 // Parse options to get subset 00138 int rank = 0, procs = 1; 00139 #ifdef MOAB_HAVE_MPI 00140 if( isParallel ) 00141 { 00142 ParallelComm*& myPcomm = _readNC->myPcomm; 00143 rank = myPcomm->proc_config().proc_rank(); 00144 procs = myPcomm->proc_config().proc_size(); 00145 } 00146 #endif 00147 if( procs > 1 ) 00148 { 00149 for( int i = 0; i < 6; i++ ) 00150 parData.gDims[i] = gDims[i]; 00151 for( int i = 0; i < 3; i++ ) 00152 parData.gPeriodic[i] = globallyPeriodic[i]; 00153 parData.partMethod = partMethod; 00154 int pdims[3]; 00155 00156 rval = ScdInterface::compute_partition( procs, rank, parData, lDims, locallyPeriodic, pdims );MB_CHK_ERR( rval ); 00157 for( int i = 0; i < 3; i++ ) 00158 parData.pDims[i] = pdims[i]; 00159 00160 dbgOut.tprintf( 1, "Partition: %dx%d (out of %dx%d)\n", lDims[3] - lDims[0] + 1, lDims[4] - lDims[1] + 1, 00161 gDims[3] - gDims[0] + 1, gDims[4] - gDims[1] + 1 ); 00162 if( 0 == rank ) 00163 dbgOut.tprintf( 1, "Contiguous chunks of size %d bytes.\n", 00164 8 * ( lDims[3] - lDims[0] + 1 ) * ( lDims[4] - lDims[1] + 1 ) ); 00165 } 00166 else 00167 { 00168 for( int i = 0; i < 6; i++ ) 00169 lDims[i] = gDims[i]; 00170 locallyPeriodic[0] = globallyPeriodic[0]; 00171 } 00172 00173 _opts.get_int_option( "IMIN", lDims[0] ); 00174 _opts.get_int_option( "IMAX", lDims[3] ); 00175 _opts.get_int_option( "JMIN", lDims[1] ); 00176 _opts.get_int_option( "JMAX", lDims[4] ); 00177 00178 // Now get actual coordinate values for vertices and cell centers 00179 lCDims[0] = lDims[0]; 00180 if( locallyPeriodic[0] ) 00181 // If locally periodic, doesn't matter what global periodicity is, # vertex coords = # elem 00182 // coords 00183 lCDims[3] = lDims[3]; 00184 else 00185 lCDims[3] = lDims[3] - 1; 00186 00187 // For FV models, will always be non-periodic in j 00188 lCDims[1] = lDims[1]; 00189 lCDims[4] = lDims[4] - 1; 00190 00191 // Resize vectors to store values later 00192 if( -1 != lDims[0] ) ilVals.resize( lDims[3] - lDims[0] + 1 ); 00193 if( -1 != lCDims[0] ) ilCVals.resize( lCDims[3] - lCDims[0] + 1 ); 00194 if( -1 != lDims[1] ) jlVals.resize( lDims[4] - lDims[1] + 1 ); 00195 if( -1 != lCDims[1] ) jlCVals.resize( lCDims[4] - lCDims[1] + 1 ); 00196 if( nTimeSteps > 0 ) tVals.resize( nTimeSteps ); 00197 00198 // Now read coord values 00199 std::map< std::string, ReadNC::VarData >::iterator vmit; 00200 if( -1 != lCDims[0] ) 00201 { 00202 if( ( vmit = varInfo.find( "lon" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 ) 00203 { 00204 rval = read_coordinate( "lon", lCDims[0], lCDims[3], ilCVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lon' variable" ); 00205 } 00206 else 00207 { 00208 MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' variable" ); 00209 } 00210 } 00211 00212 if( -1 != lCDims[1] ) 00213 { 00214 if( ( vmit = varInfo.find( "lat" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 ) 00215 { 00216 rval = read_coordinate( "lat", lCDims[1], lCDims[4], jlCVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lat' variable" ); 00217 } 00218 else 00219 { 00220 MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' variable" ); 00221 } 00222 } 00223 00224 if( -1 != lDims[0] ) 00225 { 00226 if( ( vmit = varInfo.find( "slon" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 ) 00227 { 00228 // Last column 00229 if( !locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3] ) 00230 { 00231 assert( lDims[3] == gDims[3] + 1 ); 00232 std::vector< double > dummyVar( lDims[3] - lDims[0] ); 00233 rval = read_coordinate( "slon", lDims[0], lDims[3] - 1, dummyVar ); 00234 double dif = dummyVar[1] - dummyVar[0]; 00235 std::size_t i; 00236 for( i = 0; i != dummyVar.size(); i++ ) 00237 ilVals[i] = dummyVar[i]; 00238 ilVals[i] = ilVals[i - 1] + dif; 00239 } 00240 else 00241 { 00242 rval = read_coordinate( "slon", lDims[0], lDims[3], ilVals );MB_CHK_SET_ERR( rval, "Trouble reading 'slon' variable" ); 00243 } 00244 } 00245 else 00246 { 00247 MB_SET_ERR( MB_FAILURE, "Couldn't find 'slon' variable" ); 00248 } 00249 } 00250 00251 if( -1 != lDims[1] ) 00252 { 00253 if( ( vmit = varInfo.find( "slat" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 ) 00254 { 00255 if( !isParallel || ( ( gDims[4] - gDims[1] ) == ( lDims[4] - lDims[1] ) ) ) 00256 { 00257 std::vector< double > dummyVar( lDims[4] - lDims[1] - 1 ); 00258 rval = read_coordinate( "slat", lDims[1], lDims[4] - 2, dummyVar );MB_CHK_SET_ERR( rval, "Trouble reading 'slat' variable" ); 00259 // Copy the correct piece 00260 jlVals[0] = -90.0; 00261 std::size_t i = 0; 00262 for( i = 1; i != dummyVar.size() + 1; i++ ) 00263 jlVals[i] = dummyVar[i - 1]; 00264 jlVals[i] = 90.0; // Using value of i after loop exits. 00265 } 00266 else 00267 { 00268 // If this is the first row 00269 // Need to read one less then available and read it into a dummy var 00270 if( lDims[1] == gDims[1] ) 00271 { 00272 std::vector< double > dummyVar( lDims[4] - lDims[1] ); 00273 rval = read_coordinate( "slat", lDims[1], lDims[4] - 1, dummyVar );MB_CHK_SET_ERR( rval, "Trouble reading 'slat' variable" ); 00274 // Copy the correct piece 00275 jlVals[0] = -90.0; 00276 for( int i = 1; i < lDims[4] + 1; i++ ) 00277 jlVals[i] = dummyVar[i - 1]; 00278 } 00279 // Or if it's the last row 00280 else if( lDims[4] == gDims[4] ) 00281 { 00282 std::vector< double > dummyVar( lDims[4] - lDims[1] ); 00283 rval = read_coordinate( "slat", lDims[1] - 1, lDims[4] - 2, dummyVar );MB_CHK_SET_ERR( rval, "Trouble reading 'slat' variable" ); 00284 // Copy the correct piece 00285 std::size_t i = 0; 00286 for( i = 0; i != dummyVar.size(); i++ ) 00287 jlVals[i] = dummyVar[i]; 00288 jlVals[i] = 90.0; // Using value of i after loop exits. 00289 } 00290 // It's in the middle 00291 else 00292 { 00293 rval = read_coordinate( "slat", lDims[1] - 1, lDims[4] - 1, jlVals );MB_CHK_SET_ERR( rval, "Trouble reading 'slat' variable" ); 00294 } 00295 } 00296 } 00297 else 00298 { 00299 MB_SET_ERR( MB_FAILURE, "Couldn't find 'slat' variable" ); 00300 } 00301 } 00302 00303 // Store time coordinate values in tVals 00304 if( nTimeSteps > 0 ) 00305 { 00306 if( ( vmit = varInfo.find( "time" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 ) 00307 { 00308 rval = read_coordinate( "time", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 'time' variable" ); 00309 } 00310 else if( ( vmit = varInfo.find( "t" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 ) 00311 { 00312 rval = read_coordinate( "t", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 't' variable" ); 00313 } 00314 else 00315 { 00316 // If expected time variable is not available, set dummy time coordinate values to tVals 00317 for( int t = 0; t < nTimeSteps; t++ ) 00318 tVals.push_back( (double)t ); 00319 } 00320 } 00321 00322 dbgOut.tprintf( 1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4] ); 00323 dbgOut.tprintf( 1, "%d elements, %d vertices\n", ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ), 00324 ( lDims[3] - lDims[0] + 1 ) * ( lDims[4] - lDims[1] + 1 ) ); 00325 00326 // For each variable, determine the entity location type and number of levels 00327 std::map< std::string, ReadNC::VarData >::iterator mit; 00328 for( mit = varInfo.begin(); mit != varInfo.end(); ++mit ) 00329 { 00330 ReadNC::VarData& vd = ( *mit ).second; 00331 00332 // Default entLoc is ENTLOCSET 00333 if( std::find( vd.varDims.begin(), vd.varDims.end(), tDim ) != vd.varDims.end() ) 00334 { 00335 if( ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) && 00336 ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) ) 00337 vd.entLoc = ReadNC::ENTLOCFACE; 00338 else if( ( std::find( vd.varDims.begin(), vd.varDims.end(), jDim ) != vd.varDims.end() ) && 00339 ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) ) 00340 vd.entLoc = ReadNC::ENTLOCNSEDGE; 00341 else if( ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) && 00342 ( std::find( vd.varDims.begin(), vd.varDims.end(), iDim ) != vd.varDims.end() ) ) 00343 vd.entLoc = ReadNC::ENTLOCEWEDGE; 00344 } 00345 00346 // Default numLev is 0 00347 if( std::find( vd.varDims.begin(), vd.varDims.end(), levDim ) != vd.varDims.end() ) vd.numLev = nLevels; 00348 } 00349 00350 std::vector< std::string > ijdimNames( 4 ); 00351 ijdimNames[0] = "__slon"; 00352 ijdimNames[1] = "__slat"; 00353 ijdimNames[2] = "__lon"; 00354 ijdimNames[3] = "__lat"; 00355 00356 std::string tag_name; 00357 Tag tagh; 00358 00359 // __<dim_name>_LOC_MINMAX (for slon, slat, lon and lat) 00360 for( unsigned int i = 0; i != ijdimNames.size(); i++ ) 00361 { 00362 std::vector< int > val( 2, 0 ); 00363 if( ijdimNames[i] == "__slon" ) 00364 { 00365 val[0] = lDims[0]; 00366 val[1] = lDims[3]; 00367 } 00368 else if( ijdimNames[i] == "__slat" ) 00369 { 00370 val[0] = lDims[1]; 00371 val[1] = lDims[4]; 00372 } 00373 else if( ijdimNames[i] == "__lon" ) 00374 { 00375 val[0] = lCDims[0]; 00376 val[1] = lCDims[3]; 00377 } 00378 else if( ijdimNames[i] == "__lat" ) 00379 { 00380 val[0] = lCDims[1]; 00381 val[1] = lCDims[4]; 00382 } 00383 std::stringstream ss_tag_name; 00384 ss_tag_name << ijdimNames[i] << "_LOC_MINMAX"; 00385 tag_name = ss_tag_name.str(); 00386 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 ); 00387 rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name ); 00388 if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() ); 00389 } 00390 00391 // __<dim_name>_LOC_VALS (for slon, slat, lon and lat) 00392 // Assume all have the same data type as lon (expected type is float or double) 00393 switch( varInfo["lon"].varDataType ) 00394 { 00395 case NC_FLOAT: 00396 case NC_DOUBLE: 00397 break; 00398 default: 00399 MB_SET_ERR( MB_FAILURE, "Unexpected variable data type for 'lon'" ); 00400 } 00401 00402 for( unsigned int i = 0; i != ijdimNames.size(); i++ ) 00403 { 00404 void* val = NULL; 00405 int val_len = 0; 00406 if( ijdimNames[i] == "__slon" ) 00407 { 00408 val = &ilVals[0]; 00409 val_len = ilVals.size(); 00410 } 00411 else if( ijdimNames[i] == "__slat" ) 00412 { 00413 val = &jlVals[0]; 00414 val_len = jlVals.size(); 00415 } 00416 else if( ijdimNames[i] == "__lon" ) 00417 { 00418 val = &ilCVals[0]; 00419 val_len = ilCVals.size(); 00420 } 00421 else if( ijdimNames[i] == "__lat" ) 00422 { 00423 val = &jlCVals[0]; 00424 val_len = jlCVals.size(); 00425 } 00426 00427 std::stringstream ss_tag_name; 00428 ss_tag_name << ijdimNames[i] << "_LOC_VALS"; 00429 tag_name = ss_tag_name.str(); 00430 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_DOUBLE, tagh, 00431 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name ); 00432 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 ); 00433 if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() ); 00434 } 00435 00436 // __<dim_name>_GLOBAL_MINMAX (for slon, slat, lon and lat) 00437 for( unsigned int i = 0; i != ijdimNames.size(); i++ ) 00438 { 00439 std::vector< int > val( 2, 0 ); 00440 if( ijdimNames[i] == "__slon" ) 00441 { 00442 val[0] = gDims[0]; 00443 val[1] = gDims[3]; 00444 } 00445 else if( ijdimNames[i] == "__slat" ) 00446 { 00447 val[0] = gDims[1]; 00448 val[1] = gDims[4]; 00449 } 00450 else if( ijdimNames[i] == "__lon" ) 00451 { 00452 val[0] = gCDims[0]; 00453 val[1] = gCDims[3]; 00454 } 00455 else if( ijdimNames[i] == "__lat" ) 00456 { 00457 val[0] = gCDims[1]; 00458 val[1] = gCDims[4]; 00459 } 00460 std::stringstream ss_tag_name; 00461 ss_tag_name << ijdimNames[i] << "_GLOBAL_MINMAX"; 00462 tag_name = ss_tag_name.str(); 00463 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 ); 00464 rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name ); 00465 if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() ); 00466 } 00467 00468 // Hack: create dummy variables, if needed, for dimensions with no corresponding coordinate 00469 // variables 00470 rval = create_dummy_variables();MB_CHK_SET_ERR( rval, "Failed to create dummy variables" ); 00471 00472 return MB_SUCCESS; 00473 } 00474 00475 } // namespace moab
1.7.6.1 
