MOAB: Mesh Oriented datABase  (version 5.4.1)
NCHelperFV.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines