MOAB: Mesh Oriented datABase  (version 5.4.1)
NCWriteMPAS.cpp
Go to the documentation of this file.
00001 /*
00002  * NCWriteMPAS.cpp
00003  *
00004  *  Created on: April 9, 2014
00005  */
00006 
00007 #include "NCWriteMPAS.hpp"
00008 #include "MBTagConventions.hpp"
00009 
00010 namespace moab
00011 {
00012 
00013 NCWriteMPAS::~NCWriteMPAS()
00014 {
00015     // TODO Auto-generated destructor stub
00016 }
00017 
00018 ErrorCode NCWriteMPAS::collect_mesh_info()
00019 {
00020     Interface*& mbImpl                   = _writeNC->mbImpl;
00021     std::vector< std::string >& dimNames = _writeNC->dimNames;
00022     std::vector< int >& dimLens          = _writeNC->dimLens;
00023     Tag& mGlobalIdTag                    = _writeNC->mGlobalIdTag;
00024 
00025     ErrorCode rval;
00026 
00027     // Look for time dimension
00028     std::vector< std::string >::iterator vecIt;
00029     if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "Time" ) ) != dimNames.end() )
00030         tDim = vecIt - dimNames.begin();
00031     else if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() )
00032         tDim = vecIt - dimNames.begin();
00033     else
00034     {
00035         MB_SET_ERR( MB_FAILURE, "Couldn't find 'Time' or 'time' dimension" );
00036     }
00037     nTimeSteps = dimLens[tDim];
00038 
00039     // Get number of levels
00040     if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "nVertLevels" ) ) != dimNames.end() )
00041         levDim = vecIt - dimNames.begin();
00042     else
00043     {
00044         MB_SET_ERR( MB_FAILURE, "Couldn't find 'nVertLevels' dimension" );
00045     }
00046     nLevels = dimLens[levDim];
00047 
00048     // Get local vertices
00049     rval = mbImpl->get_entities_by_dimension( _fileSet, 0, localVertsOwned );MB_CHK_SET_ERR( rval, "Trouble getting local vertices in current file set" );
00050     assert( !localVertsOwned.empty() );
00051 
00052     // Get local edges
00053     rval = mbImpl->get_entities_by_dimension( _fileSet, 1, localEdgesOwned );MB_CHK_SET_ERR( rval, "Trouble getting local edges in current file set" );
00054     // There are no edges if NO_EDGES read option is set
00055 
00056     // Get local cells
00057     rval = mbImpl->get_entities_by_dimension( _fileSet, 2, localCellsOwned );MB_CHK_SET_ERR( rval, "Trouble getting local cells in current file set" );
00058     assert( !localCellsOwned.empty() );
00059 
00060 #ifdef MOAB_HAVE_MPI
00061     bool& isParallel = _writeNC->isParallel;
00062     if( isParallel )
00063     {
00064         ParallelComm*& myPcomm = _writeNC->myPcomm;
00065         int rank               = myPcomm->proc_config().proc_rank();
00066         int procs              = myPcomm->proc_config().proc_size();
00067         if( procs > 1 )
00068         {
00069 #ifndef NDEBUG
00070             unsigned int num_local_verts = localVertsOwned.size();
00071 #endif
00072             rval = myPcomm->filter_pstatus( localVertsOwned, PSTATUS_NOT_OWNED, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Trouble getting owned vertices in current file set" );
00073 
00074             // Assume that PARALLEL_RESOLVE_SHARED_ENTS option is set
00075             // Verify that not all local vertices are owned by the last processor
00076             if( procs - 1 == rank )
00077                 assert( "PARALLEL_RESOLVE_SHARED_ENTS option is set" && localVertsOwned.size() < num_local_verts );
00078 
00079             if( !localEdgesOwned.empty() )
00080             {
00081                 rval = myPcomm->filter_pstatus( localEdgesOwned, PSTATUS_NOT_OWNED, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Trouble getting owned edges in current file set" );
00082             }
00083 
00084             rval = myPcomm->filter_pstatus( localCellsOwned, PSTATUS_NOT_OWNED, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Trouble getting owned cells in current file set" );
00085         }
00086     }
00087 #endif
00088 
00089     std::vector< int > gids( localVertsOwned.size() );
00090     rval = mbImpl->tag_get_data( mGlobalIdTag, localVertsOwned, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting global IDs on local vertices" );
00091 
00092     // Get localGidVertsOwned
00093     std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidVertsOwned ) );
00094 
00095     if( !localEdgesOwned.empty() )
00096     {
00097         gids.resize( localEdgesOwned.size() );
00098         rval = mbImpl->tag_get_data( mGlobalIdTag, localEdgesOwned, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting global IDs on local edges" );
00099 
00100         // Get localGidEdgesOwned
00101         std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidEdgesOwned ) );
00102     }
00103 
00104     gids.resize( localCellsOwned.size() );
00105     rval = mbImpl->tag_get_data( mGlobalIdTag, localCellsOwned, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting global IDs on local cells" );
00106 
00107     // Get localGidCellsOwned
00108     std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidCellsOwned ) );
00109 
00110     return MB_SUCCESS;
00111 }
00112 
00113 ErrorCode NCWriteMPAS::collect_variable_data( std::vector< std::string >& var_names, std::vector< int >& tstep_nums )
00114 {
00115     NCWriteHelper::collect_variable_data( var_names, tstep_nums );
00116 
00117     std::vector< std::string >& dimNames = _writeNC->dimNames;
00118     std::vector< int >& dimLens          = _writeNC->dimLens;
00119 
00120     // Dimension indices for other optional levels
00121     std::vector< unsigned int > opt_lev_dims;
00122 
00123     unsigned int lev_idx;
00124     std::vector< std::string >::iterator vecIt;
00125 
00126     // Get index of vertex levels P1
00127     if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "nVertLevelsP1" ) ) != dimNames.end() )
00128     {
00129         lev_idx = vecIt - dimNames.begin();
00130         opt_lev_dims.push_back( lev_idx );
00131     }
00132 
00133     // Get index of vertex levels P2
00134     if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "nVertLevelsP2" ) ) != dimNames.end() )
00135     {
00136         lev_idx = vecIt - dimNames.begin();
00137         opt_lev_dims.push_back( lev_idx );
00138     }
00139 
00140     // Get index of soil levels
00141     if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "nSoilLevels" ) ) != dimNames.end() )
00142     {
00143         lev_idx = vecIt - dimNames.begin();
00144         opt_lev_dims.push_back( lev_idx );
00145     }
00146 
00147     std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo;
00148 
00149     for( size_t i = 0; i < var_names.size(); i++ )
00150     {
00151         std::string varname                                     = var_names[i];
00152         std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( varname );
00153         if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find variable " << varname );
00154 
00155         WriteNC::VarData& currentVarData = vit->second;
00156         std::vector< int >& varDims      = currentVarData.varDims;
00157 
00158         // Skip edge variables, if there are no edges
00159         if( localEdgesOwned.empty() && currentVarData.entLoc == WriteNC::ENTLOCEDGE ) continue;
00160 
00161         // If nVertLevels dimension is not found, try other optional levels such as nVertLevelsP1
00162         if( std::find( varDims.begin(), varDims.end(), levDim ) == varDims.end() )
00163         {
00164             for( unsigned int j = 0; j < opt_lev_dims.size(); j++ )
00165             {
00166                 if( std::find( varDims.begin(), varDims.end(), opt_lev_dims[j] ) != varDims.end() )
00167                 {
00168                     currentVarData.numLev = dimLens[opt_lev_dims[j]];
00169                     break;
00170                 }
00171             }
00172         }
00173 
00174         // Skip set variables, which were already processed in
00175         // NCWriteHelper::collect_variable_data()
00176         if( WriteNC::ENTLOCSET == currentVarData.entLoc ) continue;
00177 
00178         // Set up writeStarts and writeCounts (maximum number of dimensions is 3)
00179         currentVarData.writeStarts.resize( 3 );
00180         currentVarData.writeCounts.resize( 3 );
00181         unsigned int dim_idx = 0;
00182 
00183         // First: Time
00184         if( currentVarData.has_tsteps )
00185         {
00186             // Non-set variables with timesteps
00187             // 3 dimensions like (Time, nCells, nVertLevels)
00188             // 2 dimensions like (Time, nCells)
00189             assert( 3 == varDims.size() || 2 == varDims.size() );
00190 
00191             // Time should be the first dimension
00192             assert( tDim == varDims[0] );
00193 
00194             currentVarData.writeStarts[dim_idx] = 0;  // This value is timestep dependent, will be set later
00195             currentVarData.writeCounts[dim_idx] = 1;
00196             dim_idx++;
00197         }
00198         else
00199         {
00200             // Non-set variables without timesteps
00201             // 2 dimensions like (nCells, nVertLevels)
00202             // 1 dimension like (nCells)
00203             assert( 2 == varDims.size() || 1 == varDims.size() );
00204         }
00205 
00206         // Next: nVertices / nCells / nEdges
00207         switch( currentVarData.entLoc )
00208         {
00209             case WriteNC::ENTLOCVERT:
00210                 // Vertices
00211                 // Start from the first localGidVerts
00212                 // Actually, this will be reset later for writing
00213                 currentVarData.writeStarts[dim_idx] = localGidVertsOwned[0] - 1;
00214                 currentVarData.writeCounts[dim_idx] = localGidVertsOwned.size();
00215                 break;
00216             case WriteNC::ENTLOCFACE:
00217                 // Faces
00218                 // Start from the first localGidCells
00219                 // Actually, this will be reset later for writing
00220                 currentVarData.writeStarts[dim_idx] = localGidCellsOwned[0] - 1;
00221                 currentVarData.writeCounts[dim_idx] = localGidCellsOwned.size();
00222                 break;
00223             case WriteNC::ENTLOCEDGE:
00224                 // Edges
00225                 // Start from the first localGidEdges
00226                 // Actually, this will be reset later for writing
00227                 currentVarData.writeStarts[dim_idx] = localGidEdgesOwned[0] - 1;
00228                 currentVarData.writeCounts[dim_idx] = localGidEdgesOwned.size();
00229                 break;
00230             default:
00231                 MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << varname );
00232         }
00233         dim_idx++;
00234 
00235         // Finally: nVertLevels or other optional levels, it is possible that there is no
00236         // level dimension (numLev is 0) for this non-set variable, e.g. (Time, nCells)
00237         if( currentVarData.numLev > 0 )
00238         {
00239             // Non-set variables with levels
00240             // 3 dimensions like (Time, nCells, nVertLevels)
00241             // 2 dimensions like (nCells, nVertLevels)
00242             assert( 3 == varDims.size() || 2 == varDims.size() );
00243 
00244             currentVarData.writeStarts[dim_idx] = 0;
00245             currentVarData.writeCounts[dim_idx] = currentVarData.numLev;
00246             dim_idx++;
00247         }
00248         else
00249         {
00250             // Non-set variables without levels
00251             // 2 dimensions like (Time, nCells)
00252             // 1 dimension like (nCells)
00253             assert( 2 == varDims.size() || 1 == varDims.size() );
00254         }
00255 
00256         // Get variable size
00257         currentVarData.sz = 1;
00258         for( std::size_t idx = 0; idx < dim_idx; idx++ )
00259             currentVarData.sz *= currentVarData.writeCounts[idx];
00260     }
00261 
00262     return MB_SUCCESS;
00263 }
00264 
00265 ErrorCode NCWriteMPAS::write_nonset_variables( std::vector< WriteNC::VarData >& vdatas, std::vector< int >& tstep_nums )
00266 {
00267     Interface*& mbImpl = _writeNC->mbImpl;
00268 
00269     int success;
00270 
00271     // For each variable tag in the indexed lists, write a time step data
00272     for( unsigned int i = 0; i < vdatas.size(); i++ )
00273     {
00274         WriteNC::VarData& variableData = vdatas[i];
00275 
00276         // Skip edge variables, if there are no edges
00277         if( localEdgesOwned.empty() && variableData.entLoc == WriteNC::ENTLOCEDGE ) continue;
00278 
00279         // Get local owned entities of this variable
00280         Range* pLocalEntsOwned    = NULL;
00281         Range* pLocalGidEntsOwned = NULL;
00282         switch( variableData.entLoc )
00283         {
00284             case WriteNC::ENTLOCVERT:
00285                 // Vertices
00286                 pLocalEntsOwned    = &localVertsOwned;
00287                 pLocalGidEntsOwned = &localGidVertsOwned;
00288                 break;
00289             case WriteNC::ENTLOCEDGE:
00290                 // Edges
00291                 pLocalEntsOwned    = &localEdgesOwned;
00292                 pLocalGidEntsOwned = &localGidEdgesOwned;
00293                 break;
00294             case WriteNC::ENTLOCFACE:
00295                 // Cells
00296                 pLocalEntsOwned    = &localCellsOwned;
00297                 pLocalGidEntsOwned = &localGidCellsOwned;
00298                 break;
00299             default:
00300                 MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << variableData.varName );
00301         }
00302 
00303         unsigned int num_timesteps;
00304         unsigned int ents_idx = 0;
00305         if( variableData.has_tsteps )
00306         {
00307             // Non-set variables with timesteps
00308             // 3 dimensions like (Time, nCells, nVertLevels)
00309             // 2 dimensions like (Time, nCells)
00310             num_timesteps = tstep_nums.size();
00311             ents_idx++;
00312         }
00313         else
00314         {
00315             // Non-set variables without timesteps
00316             // 2 dimensions like (nCells, nVertLevels)
00317             // 1 dimension like (nCells)
00318             num_timesteps = 1;
00319         }
00320 
00321         unsigned int num_lev;
00322         if( variableData.numLev > 0 )
00323         {
00324             // Non-set variables with levels
00325             // 3 dimensions like (Time, nCells, nVertLevels)
00326             // 2 dimensions like (nCells, nVertLevels)
00327             num_lev = variableData.numLev;
00328         }
00329         else
00330         {
00331             // Non-set variables without levels
00332             // 2 dimensions like (Time, nCells)
00333             // 1 dimension like (nCells)
00334             num_lev = 1;
00335         }
00336 
00337         for( unsigned int t = 0; t < num_timesteps; t++ )
00338         {
00339             // We will write one time step, and count will be one; start will be different
00340             // Use tag_get_data instead of tag_iterate to copy tag data, as localEntsOwned
00341             // might not be contiguous.
00342             if( tDim == variableData.varDims[0] ) variableData.writeStarts[0] = t;  // This is start for time
00343             std::vector< double > tag_data( pLocalEntsOwned->size() * num_lev );
00344             ErrorCode rval = mbImpl->tag_get_data( variableData.varTags[t], *pLocalEntsOwned, &tag_data[0] );MB_CHK_SET_ERR( rval, "Trouble getting tag data on owned entities" );
00345 
00346 #ifdef MOAB_HAVE_PNETCDF
00347             size_t nb_writes = pLocalGidEntsOwned->psize();
00348             std::vector< int > requests( nb_writes ), statuss( nb_writes );
00349             size_t idxReq = 0;
00350 #endif
00351 
00352             // Now write copied tag data
00353             // Use nonblocking put (request aggregation)
00354             switch( variableData.varDataType )
00355             {
00356                 case NC_DOUBLE: {
00357                     size_t indexInDoubleArray = 0;
00358                     size_t ic                 = 0;
00359                     for( Range::pair_iterator pair_iter = pLocalGidEntsOwned->pair_begin();
00360                          pair_iter != pLocalGidEntsOwned->pair_end(); ++pair_iter, ic++ )
00361                     {
00362                         EntityHandle starth                = pair_iter->first;
00363                         EntityHandle endh                  = pair_iter->second;
00364                         variableData.writeStarts[ents_idx] = (NCDF_SIZE)( starth - 1 );
00365                         variableData.writeCounts[ents_idx] = (NCDF_SIZE)( endh - starth + 1 );
00366 
00367                         // Do a partial write, in each subrange
00368 #ifdef MOAB_HAVE_PNETCDF
00369                         // Wait outside this loop
00370                         success =
00371                             NCFUNCREQP( _vara_double )( _fileId, variableData.varId, &( variableData.writeStarts[0] ),
00372                                                         &( variableData.writeCounts[0] ),
00373                                                         &( tag_data[indexInDoubleArray] ), &requests[idxReq++] );
00374 #else
00375                         success = NCFUNCAP(
00376                             _vara_double )( _fileId, variableData.varId, &( variableData.writeStarts[0] ),
00377                                             &( variableData.writeCounts[0] ), &( tag_data[indexInDoubleArray] ) );
00378 #endif
00379                         if( success )
00380                             MB_SET_ERR( MB_FAILURE,
00381                                         "Failed to write double data in a loop for variable " << variableData.varName );
00382                         // We need to increment the index in double array for the
00383                         // next subrange
00384                         indexInDoubleArray += ( endh - starth + 1 ) * num_lev;
00385                     }
00386                     assert( ic == pLocalGidEntsOwned->psize() );
00387 #ifdef MOAB_HAVE_PNETCDF
00388                     success = ncmpi_wait_all( _fileId, requests.size(), &requests[0], &statuss[0] );
00389                     if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
00390 #endif
00391                     break;
00392                 }
00393                 default:
00394                     MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing non-double data is not implemented yet" );
00395             }
00396         }
00397     }
00398 
00399     return MB_SUCCESS;
00400 }
00401 
00402 } /* namespace moab */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines