MOAB: Mesh Oriented datABase  (version 5.2.1)
TempestOnlineMapIO.cpp
Go to the documentation of this file.
00001 /*
00002  * =====================================================================================
00003  *
00004  *       Filename:  TempestOnlineMapIO.cpp
00005  *
00006  *    Description:  All I/O implementations related to TempestOnlineMap
00007  *
00008  *        Version:  1.0
00009  *        Created:  02/06/2021 02:35:41
00010  *
00011  *         Author:  Vijay S. Mahadevan (vijaysm), mahadevan@anl.gov
00012  *        Company:  Argonne National Lab
00013  *
00014  * =====================================================================================
00015  */
00016 
00017 #include "FiniteElementTools.h"
00018 #include "moab/Remapping/TempestOnlineMap.hpp"
00019 
00020 #ifdef MOAB_HAVE_NETCDFPAR
00021 #include "netcdfcpp_par.hpp"
00022 #else
00023 #include "netcdfcpp.h"
00024 #endif
00025 
00026 #ifdef MOAB_HAVE_MPI
00027 
00028 #define MPI_CHK_ERR( err )                                          \
00029     if( err )                                                       \
00030     {                                                               \
00031         std::cout << "MPI Failure. ErrorCode (" << ( err ) << ") "; \
00032         std::cout << "\nMPI Aborting... \n";                        \
00033         return moab::MB_FAILURE;                                    \
00034     }
00035 
00036 int moab::TempestOnlineMap::rearrange_arrays_by_dofs(
00037     const std::vector< unsigned int >& gdofmap, DataArray1D< double >& vecFaceArea, DataArray1D< double >& dCenterLon,
00038     DataArray1D< double >& dCenterLat, DataArray2D< double >& dVertexLon, DataArray2D< double >& dVertexLat,
00039     unsigned& N,  // will have the local, after
00040     int nv, int& maxdof )
00041 {
00042     // first decide maxdof, for partitioning
00043     unsigned int localmax = 0;
00044     for( unsigned i = 0; i < N; i++ )
00045         if( gdofmap[i] > localmax ) localmax = gdofmap[i];
00046 
00047     // decide partitioning based on maxdof/size
00048     MPI_Allreduce( &localmax, &maxdof, 1, MPI_INT, MPI_MAX, m_pcomm->comm() );
00049     // maxdof is 0 based, so actual number is +1
00050     // maxdof
00051     int size_per_task = ( maxdof + 1 ) / size;  // based on this, processor to process dof x is x/size_per_task
00052     // so we decide to reorder by actual dof, such that task 0 has dofs from [0 to size_per_task), etc
00053     moab::TupleList tl;
00054     unsigned numr = 2 * nv + 3;         //  doubles: area, centerlon, center lat, nv (vertex lon, vertex lat)
00055     tl.initialize( 2, 0, 0, numr, N );  // to proc, dof, then
00056     tl.enableWriteAccess();
00057     // populate
00058     for( unsigned i = 0; i < N; i++ )
00059     {
00060         int gdof    = gdofmap[i];
00061         int to_proc = gdof / size_per_task;
00062         if( to_proc >= size ) to_proc = size - 1;  // the last ones got to last proc
00063         int n                  = tl.get_n();
00064         tl.vi_wr[2 * n]        = to_proc;
00065         tl.vi_wr[2 * n + 1]    = gdof;
00066         tl.vr_wr[n * numr]     = vecFaceArea[i];
00067         tl.vr_wr[n * numr + 1] = dCenterLon[i];
00068         tl.vr_wr[n * numr + 2] = dCenterLat[i];
00069         for( int j = 0; j < nv; j++ )
00070         {
00071             tl.vr_wr[n * numr + 3 + j]      = dVertexLon[i][j];
00072             tl.vr_wr[n * numr + 3 + nv + j] = dVertexLat[i][j];
00073         }
00074         tl.inc_n();
00075     }
00076 
00077     // now do the heavy communication
00078     ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tl, 0 );
00079 
00080     // after communication, on each processor we should have tuples coming in
00081     // still need to order by global dofs; then rearrange input vectors
00082     moab::TupleList::buffer sort_buffer;
00083     sort_buffer.buffer_init( tl.get_n() );
00084     tl.sort( 1, &sort_buffer );
00085     // count how many are unique, and collapse
00086     int nb_unique = 1;
00087     for( unsigned i = 0; i < tl.get_n() - 1; i++ )
00088     {
00089         if( tl.vi_wr[2 * i + 1] != tl.vi_wr[2 * i + 3] ) nb_unique++;
00090     }
00091     vecFaceArea.Allocate( nb_unique );
00092     dCenterLon.Allocate( nb_unique );
00093     dCenterLat.Allocate( nb_unique );
00094     dVertexLon.Allocate( nb_unique, nv );
00095     dVertexLat.Allocate( nb_unique, nv );
00096     int current_size = 1;
00097     vecFaceArea[0]   = tl.vr_wr[0];
00098     dCenterLon[0]    = tl.vr_wr[1];
00099     dCenterLat[0]    = tl.vr_wr[2];
00100     for( int j = 0; j < nv; j++ )
00101     {
00102         dVertexLon[0][j] = tl.vr_wr[3 + j];
00103         dVertexLat[0][j] = tl.vr_wr[3 + nv + j];
00104     }
00105     for( unsigned i = 0; i < tl.get_n() - 1; i++ )
00106     {
00107         if( tl.vi_wr[2 * i + 1] != tl.vi_wr[2 * i + 3] )
00108         {
00109             vecFaceArea[current_size] = tl.vr_wr[i * numr];
00110             dCenterLon[current_size]  = tl.vr_wr[i * numr + 1];
00111             dCenterLat[current_size]  = tl.vr_wr[i * numr + 2];
00112             for( int j = 0; j < nv; j++ )
00113             {
00114                 dVertexLon[current_size][j] = tl.vr_wr[i * numr + 3 + j];
00115                 dVertexLat[current_size][j] = tl.vr_wr[i * numr + 3 + nv + j];
00116             }
00117             current_size++;
00118         }
00119     }
00120 
00121     N = current_size;  // or nb_unique, should be the same
00122     return 0;
00123 }
00124 #endif
00125 
00126 ///////////////////////////////////////////////////////////////////////////////
00127 
00128 moab::ErrorCode moab::TempestOnlineMap::WriteParallelMap( const std::string& strFilename )
00129 {
00130     moab::ErrorCode rval;
00131 
00132     size_t lastindex      = strFilename.find_last_of( "." );
00133     std::string extension = strFilename.substr( lastindex + 1, strFilename.size() );
00134 
00135     // Write the map file to disk in parallel
00136     if( extension == "nc" )
00137     {
00138         /* Invoke the actual call to write the parallel map to disk in SCRIP format */
00139         rval = this->WriteSCRIPMapFile( strFilename.c_str() );MB_CHK_ERR( rval );
00140     }
00141     else
00142     {
00143         /* Write to the parallel H5M format */
00144         rval = this->WriteHDF5MapFile( strFilename.c_str() );MB_CHK_ERR( rval );
00145     }
00146 
00147     return rval;
00148 }
00149 
00150 ///////////////////////////////////////////////////////////////////////////////
00151 
00152 moab::ErrorCode moab::TempestOnlineMap::WriteSCRIPMapFile( const std::string& strFilename )
00153 {
00154     NcError error( NcError::silent_nonfatal );
00155 
00156 #ifdef MOAB_HAVE_NETCDFPAR
00157     bool is_independent = true;
00158     ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcFile::Replace, NcFile::Netcdf4 );
00159     // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
00160 #else
00161     NcFile ncMap( strFilename.c_str(), NcFile::Replace );
00162 #endif
00163 
00164     if( !ncMap.is_valid() ) { _EXCEPTION1( "Unable to open output map file \"%s\"", strFilename.c_str() ); }
00165 
00166     // Attributes
00167     ncMap.add_att( "Title", "MOAB-TempestRemap Online Regridding Weight Generator" );
00168 
00169     /**
00170      * Need to get the global maximum of number of vertices per element
00171      * Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
00172      *dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
00173      *when writing this out, other processes may have a different size for the same array. This is
00174      *hence a mess to consolidate in h5mtoscrip eventually.
00175      **/
00176 
00177     /* Let us compute all relevant data for the current original source mesh on the process */
00178     DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
00179     DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
00180     DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
00181     if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
00182     {
00183         this->InitializeCoordinatesFromMeshFV( *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon,
00184                                                dSourceVertexLat, false /* fLatLon = false */,
00185                                                m_remapper->max_source_edges );
00186 
00187         vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
00188         for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
00189             vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
00190     }
00191     else
00192     {
00193         DataArray3D< double > dataGLLJacobianSrc;
00194         this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
00195                                                dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
00196 
00197         // Generate the continuous Jacobian for input mesh
00198         GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
00199 
00200         if( m_srcDiscType == DiscretizationType_CGLL )
00201         {
00202             GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, m_meshInput->vecFaceArea );
00203         }
00204         else
00205         {
00206             GenerateDiscontinuousJacobian( dataGLLJacobianSrc, m_meshInput->vecFaceArea );
00207         }
00208 
00209         vecSourceFaceArea.Allocate( m_nTotDofs_Src );
00210         int offset = 0;
00211         for( size_t e = 0; e < m_meshInput->faces.size(); e++ )
00212         {
00213             for( int s = 0; s < m_nDofsPEl_Src; s++ )
00214             {
00215                 for( int t = 0; t < m_nDofsPEl_Src; t++ )
00216                 {
00217                     vecSourceFaceArea[srccol_dtoc_dofmap[offset + s * m_nDofsPEl_Src + t]] =
00218                         dataGLLJacobianSrc[s][t][e];
00219                 }
00220             }
00221             offset += m_nDofsPEl_Src * m_nDofsPEl_Src;
00222         }
00223     }
00224 
00225     if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
00226     {
00227         this->InitializeCoordinatesFromMeshFV( *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon,
00228                                                dTargetVertexLat, false /* fLatLon = false */,
00229                                                m_remapper->max_target_edges );
00230 
00231         vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
00232         for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
00233         {
00234             vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
00235         }
00236     }
00237     else
00238     {
00239         DataArray3D< double > dataGLLJacobianDest;
00240         this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
00241                                                dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
00242 
00243         // Generate the continuous Jacobian for input mesh
00244         GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
00245 
00246         if( m_destDiscType == DiscretizationType_CGLL )
00247         {
00248             GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, m_meshOutput->vecFaceArea );
00249         }
00250         else
00251         {
00252             GenerateDiscontinuousJacobian( dataGLLJacobianDest, m_meshOutput->vecFaceArea );
00253         }
00254 
00255         vecTargetFaceArea.Allocate( m_nTotDofs_Dest );
00256         int offset = 0;
00257         for( size_t e = 0; e < m_meshOutput->faces.size(); e++ )
00258         {
00259             for( int s = 0; s < m_nDofsPEl_Dest; s++ )
00260             {
00261                 for( int t = 0; t < m_nDofsPEl_Dest; t++ )
00262                 {
00263                     vecTargetFaceArea[row_dtoc_dofmap[offset + s * m_nDofsPEl_Dest + t]] = dataGLLJacobianDest[s][t][e];
00264                 }
00265             }
00266             offset += m_nDofsPEl_Dest * m_nDofsPEl_Dest;
00267         }
00268     }
00269 
00270     // Map dimensions
00271     unsigned nA = ( vecSourceFaceArea.GetRows() );
00272     unsigned nB = ( vecTargetFaceArea.GetRows() );
00273 
00274     // Number of nodes per Face
00275     int nSourceNodesPerFace = dSourceVertexLon.GetColumns();
00276     int nTargetNodesPerFace = dTargetVertexLon.GetColumns();
00277     // first move data if in parallel
00278 #if defined( MOAB_HAVE_MPI )
00279     // if (size > 1)
00280     {
00281         int maxdof;  // output; arrays will be re-distributed in chunks [maxdof/size]
00282         int ierr = rearrange_arrays_by_dofs( srccol_gdofmap, vecSourceFaceArea, dSourceCenterLon, dSourceCenterLat,
00283                                              dSourceVertexLon, dSourceVertexLat, nA, nSourceNodesPerFace,
00284                                              maxdof );  // now nA will be close to maxdof/size
00285         if( ierr != 0 ) { _EXCEPTION1( "Unable to arrange source data %d ", nA ); }
00286         // rearrange target data: (nB)
00287         //
00288         ierr = rearrange_arrays_by_dofs( row_gdofmap, vecTargetFaceArea, dTargetCenterLon, dTargetCenterLat,
00289                                          dTargetVertexLon, dTargetVertexLat, nB, nTargetNodesPerFace,
00290                                          maxdof );  // now nA will be close to maxdof/size
00291         if( ierr != 0 ) { _EXCEPTION1( "Unable to arrange target data %d ", nB ); }
00292     }
00293 #endif
00294 
00295     // Number of non-zeros in the remap matrix operator
00296     int nS = m_weightMatrix.nonZeros();
00297 
00298 #if defined( MOAB_HAVE_MPI ) && defined( MOAB_HAVE_NETCDFPAR )
00299     int locbuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
00300     int offbuf[3] = { 0, 0, 0 };
00301     int globuf[5] = { 0, 0, 0, 0, 0 };
00302     MPI_Scan( locbuf, offbuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
00303     MPI_Allreduce( locbuf, globuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
00304     MPI_Allreduce( &locbuf[3], &globuf[3], 2, MPI_INT, MPI_MAX, m_pcomm->comm() );
00305 
00306     // MPI_Scan is inclusive of data in current rank; modify accordingly.
00307     offbuf[0] -= nA;
00308     offbuf[1] -= nB;
00309     offbuf[2] -= nS;
00310 
00311 #else
00312     int offbuf[3] = { 0, 0, 0 };
00313     int globuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
00314 #endif
00315 
00316     // Write output dimensions entries
00317     unsigned nSrcGridDims = ( m_vecSourceDimSizes.size() );
00318     unsigned nDstGridDims = ( m_vecTargetDimSizes.size() );
00319 
00320     NcDim* dimSrcGridRank = ncMap.add_dim( "src_grid_rank", nSrcGridDims );
00321     NcDim* dimDstGridRank = ncMap.add_dim( "dst_grid_rank", nDstGridDims );
00322 
00323     NcVar* varSrcGridDims = ncMap.add_var( "src_grid_dims", ncInt, dimSrcGridRank );
00324     NcVar* varDstGridDims = ncMap.add_var( "dst_grid_dims", ncInt, dimDstGridRank );
00325 
00326 #ifdef MOAB_HAVE_NETCDFPAR
00327     ncMap.enable_var_par_access( varSrcGridDims, is_independent );
00328     ncMap.enable_var_par_access( varDstGridDims, is_independent );
00329 #endif
00330 
00331     char szDim[64];
00332     if( ( nSrcGridDims == 1 ) && ( m_vecSourceDimSizes[0] != (int)nA ) )
00333     {
00334         varSrcGridDims->put( &globuf[0], 1 );
00335         varSrcGridDims->add_att( "name0", "num_dof" );
00336     }
00337     else
00338     {
00339         for( unsigned i = 0; i < m_vecSourceDimSizes.size(); i++ )
00340         {
00341             int tmp = ( i == 0 ? globuf[0] : m_vecSourceDimSizes[i] );
00342             varSrcGridDims->set_cur( nSrcGridDims - i - 1 );
00343             varSrcGridDims->put( &( tmp ), 1 );
00344         }
00345 
00346         for( unsigned i = 0; i < m_vecSourceDimSizes.size(); i++ )
00347         {
00348             sprintf( szDim, "name%i", i );
00349             varSrcGridDims->add_att( szDim, m_vecSourceDimNames[nSrcGridDims - i - 1].c_str() );
00350         }
00351     }
00352 
00353     if( ( nDstGridDims == 1 ) && ( m_vecTargetDimSizes[0] != (int)nB ) )
00354     {
00355         varDstGridDims->put( &globuf[1], 1 );
00356         varDstGridDims->add_att( "name0", "num_dof" );
00357     }
00358     else
00359     {
00360         for( unsigned i = 0; i < m_vecTargetDimSizes.size(); i++ )
00361         {
00362             int tmp = ( i == 0 ? globuf[1] : m_vecTargetDimSizes[i] );
00363             varDstGridDims->set_cur( nDstGridDims - i - 1 );
00364             varDstGridDims->put( &( tmp ), 1 );
00365         }
00366 
00367         for( unsigned i = 0; i < m_vecTargetDimSizes.size(); i++ )
00368         {
00369             sprintf( szDim, "name%i", i );
00370             varDstGridDims->add_att( szDim, m_vecTargetDimNames[nDstGridDims - i - 1].c_str() );
00371         }
00372     }
00373 
00374     // Source and Target mesh resolutions
00375     NcDim* dimNA = ncMap.add_dim( "n_a", globuf[0] );
00376     NcDim* dimNB = ncMap.add_dim( "n_b", globuf[1] );
00377 
00378     // Number of nodes per Face
00379     NcDim* dimNVA = ncMap.add_dim( "nv_a", globuf[3] );
00380     NcDim* dimNVB = ncMap.add_dim( "nv_b", globuf[4] );
00381 
00382     // Write coordinates
00383     NcVar* varYCA = ncMap.add_var( "yc_a", ncDouble, dimNA );
00384     NcVar* varYCB = ncMap.add_var( "yc_b", ncDouble, dimNB );
00385 
00386     NcVar* varXCA = ncMap.add_var( "xc_a", ncDouble, dimNA );
00387     NcVar* varXCB = ncMap.add_var( "xc_b", ncDouble, dimNB );
00388 
00389     NcVar* varYVA = ncMap.add_var( "yv_a", ncDouble, dimNA, dimNVA );
00390     NcVar* varYVB = ncMap.add_var( "yv_b", ncDouble, dimNB, dimNVB );
00391 
00392     NcVar* varXVA = ncMap.add_var( "xv_a", ncDouble, dimNA, dimNVA );
00393     NcVar* varXVB = ncMap.add_var( "xv_b", ncDouble, dimNB, dimNVB );
00394 
00395 #ifdef MOAB_HAVE_NETCDFPAR
00396     ncMap.enable_var_par_access( varYCA, is_independent );
00397     ncMap.enable_var_par_access( varYCB, is_independent );
00398     ncMap.enable_var_par_access( varXCA, is_independent );
00399     ncMap.enable_var_par_access( varXCB, is_independent );
00400     ncMap.enable_var_par_access( varYVA, is_independent );
00401     ncMap.enable_var_par_access( varYVB, is_independent );
00402     ncMap.enable_var_par_access( varXVA, is_independent );
00403     ncMap.enable_var_par_access( varXVB, is_independent );
00404 #endif
00405 
00406     varYCA->add_att( "units", "degrees" );
00407     varYCB->add_att( "units", "degrees" );
00408 
00409     varXCA->add_att( "units", "degrees" );
00410     varXCB->add_att( "units", "degrees" );
00411 
00412     varYVA->add_att( "units", "degrees" );
00413     varYVB->add_att( "units", "degrees" );
00414 
00415     varXVA->add_att( "units", "degrees" );
00416     varXVB->add_att( "units", "degrees" );
00417 
00418     // Verify dimensionality
00419     if( dSourceCenterLon.GetRows() != nA ) { _EXCEPTIONT( "Mismatch between dSourceCenterLon and nA" ); }
00420     if( dSourceCenterLat.GetRows() != nA ) { _EXCEPTIONT( "Mismatch between dSourceCenterLat and nA" ); }
00421     if( dTargetCenterLon.GetRows() != nB ) { _EXCEPTIONT( "Mismatch between dTargetCenterLon and nB" ); }
00422     if( dTargetCenterLat.GetRows() != nB ) { _EXCEPTIONT( "Mismatch between dTargetCenterLat and nB" ); }
00423     if( dSourceVertexLon.GetRows() != nA ) { _EXCEPTIONT( "Mismatch between dSourceVertexLon and nA" ); }
00424     if( dSourceVertexLat.GetRows() != nA ) { _EXCEPTIONT( "Mismatch between dSourceVertexLat and nA" ); }
00425     if( dTargetVertexLon.GetRows() != nB ) { _EXCEPTIONT( "Mismatch between dTargetVertexLon and nB" ); }
00426     if( dTargetVertexLat.GetRows() != nB ) { _EXCEPTIONT( "Mismatch between dTargetVertexLat and nB" ); }
00427 
00428     varYCA->set_cur( (long)offbuf[0] );
00429     varYCA->put( &( dSourceCenterLat[0] ), nA );
00430     varYCB->set_cur( (long)offbuf[1] );
00431     varYCB->put( &( dTargetCenterLat[0] ), nB );
00432 
00433     varXCA->set_cur( (long)offbuf[0] );
00434     varXCA->put( &( dSourceCenterLon[0] ), nA );
00435     varXCB->set_cur( (long)offbuf[1] );
00436     varXCB->put( &( dTargetCenterLon[0] ), nB );
00437 
00438     varYVA->set_cur( (long)offbuf[0] );
00439     varYVA->put( &( dSourceVertexLat[0][0] ), nA, nSourceNodesPerFace );
00440     varYVB->set_cur( (long)offbuf[1] );
00441     varYVB->put( &( dTargetVertexLat[0][0] ), nB, nTargetNodesPerFace );
00442 
00443     varXVA->set_cur( (long)offbuf[0] );
00444     varXVA->put( &( dSourceVertexLon[0][0] ), nA, nSourceNodesPerFace );
00445     varXVB->set_cur( (long)offbuf[1] );
00446     varXVB->put( &( dTargetVertexLon[0][0] ), nB, nTargetNodesPerFace );
00447 
00448     // Write areas
00449     NcVar* varAreaA = ncMap.add_var( "area_a", ncDouble, dimNA );
00450 #ifdef MOAB_HAVE_NETCDFPAR
00451     ncMap.enable_var_par_access( varAreaA, is_independent );
00452 #endif
00453     varAreaA->set_cur( (long)offbuf[0] );
00454     varAreaA->put( &( vecSourceFaceArea[0] ), nA );
00455 
00456     NcVar* varAreaB = ncMap.add_var( "area_b", ncDouble, dimNB );
00457 #ifdef MOAB_HAVE_NETCDFPAR
00458     ncMap.enable_var_par_access( varAreaB, is_independent );
00459 #endif
00460     varAreaB->set_cur( (long)offbuf[1] );
00461     varAreaB->put( &( vecTargetFaceArea[0] ), nB );
00462 
00463     // Write frac
00464     DataArray1D< double > dFracA( nA );
00465     for( unsigned i = 0; i < nA; i++ )
00466     {
00467         dFracA[i] = 1.0;
00468     }
00469     NcVar* varFracA = ncMap.add_var( "frac_a", ncDouble, dimNA );
00470 #ifdef MOAB_HAVE_NETCDFPAR
00471     ncMap.enable_var_par_access( varFracA, is_independent );
00472 #endif
00473     varFracA->set_cur( (long)offbuf[0] );
00474     varFracA->put( &( dFracA[0] ), nA );
00475 
00476     DataArray1D< double > dFracB( nB );
00477     for( unsigned i = 0; i < nB; i++ )
00478     {
00479         dFracB[i] = 1.0;
00480     }
00481     NcVar* varFracB = ncMap.add_var( "frac_b", ncDouble, dimNB );
00482 #ifdef MOAB_HAVE_NETCDFPAR
00483     ncMap.enable_var_par_access( varFracB, is_independent );
00484 #endif
00485     varFracB->set_cur( (long)offbuf[1] );
00486     varFracB->put( &( dFracB[0] ), nB );
00487 
00488     // Write SparseMatrix entries
00489     DataArray1D< int > vecRow( nS );
00490     DataArray1D< int > vecCol( nS );
00491     DataArray1D< double > vecS( nS );
00492 
00493     int offset = 0;
00494     for( int i = 0; i < m_weightMatrix.outerSize(); ++i )
00495     {
00496         for( WeightMatrix::InnerIterator it( m_weightMatrix, i ); it; ++it )
00497         {
00498             vecRow[offset] = 1 + this->GetRowGlobalDoF( it.row() );  // row index
00499             vecCol[offset] = 1 + this->GetColGlobalDoF( it.col() );  // col index
00500             vecS[offset]   = it.value();                             // value
00501             offset++;
00502         }
00503     }
00504 
00505     // Load in data
00506     NcDim* dimNS = ncMap.add_dim( "n_s", globuf[2] );
00507 
00508     NcVar* varRow = ncMap.add_var( "row", ncInt, dimNS );
00509     NcVar* varCol = ncMap.add_var( "col", ncInt, dimNS );
00510     NcVar* varS   = ncMap.add_var( "S", ncDouble, dimNS );
00511 #ifdef MOAB_HAVE_NETCDFPAR
00512     ncMap.enable_var_par_access( varRow, is_independent );
00513     ncMap.enable_var_par_access( varCol, is_independent );
00514     ncMap.enable_var_par_access( varS, is_independent );
00515 #endif
00516 
00517     varRow->set_cur( (long)offbuf[2] );
00518     varRow->put( vecRow, nS );
00519 
00520     varCol->set_cur( (long)offbuf[2] );
00521     varCol->put( vecCol, nS );
00522 
00523     varS->set_cur( (long)offbuf[2] );
00524     varS->put( &( vecS[0] ), nS );
00525 
00526     // Add global attributes
00527     // std::map<std::string, std::string>::const_iterator iterAttributes =
00528     //     mapAttributes.begin();
00529     // for (; iterAttributes != mapAttributes.end(); iterAttributes++) {
00530     //     ncMap.add_att(
00531     //         iterAttributes->first.c_str(),
00532     //         iterAttributes->second.c_str());
00533     // }
00534 
00535     ncMap.close();
00536 
00537     return moab::MB_SUCCESS;
00538 }
00539 
00540 ///////////////////////////////////////////////////////////////////////////////
00541 
00542 moab::ErrorCode moab::TempestOnlineMap::WriteHDF5MapFile( const std::string& strOutputFile )
00543 {
00544     moab::ErrorCode rval;
00545 
00546     /**
00547      * Need to get the global maximum of number of vertices per element
00548      * Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
00549      *dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
00550      *when writing this out, other processes may have a different size for the same array. This is
00551      *hence a mess to consolidate in h5mtoscrip eventually.
00552      **/
00553 
00554     /* Let us compute all relevant data for the current original source mesh on the process */
00555     DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
00556     DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
00557     DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
00558     if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
00559     {
00560         this->InitializeCoordinatesFromMeshFV( *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon,
00561                                                dSourceVertexLat, false /* fLatLon = false */,
00562                                                m_remapper->max_source_edges );
00563 
00564         vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
00565         for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
00566             vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
00567     }
00568     else
00569     {
00570         DataArray3D< double > dataGLLJacobianSrc;
00571         this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
00572                                                dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
00573 
00574         // Generate the continuous Jacobian for input mesh
00575         GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
00576 
00577         if( m_srcDiscType == DiscretizationType_CGLL )
00578         {
00579             GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, m_meshInput->vecFaceArea );
00580         }
00581         else
00582         {
00583             GenerateDiscontinuousJacobian( dataGLLJacobianSrc, m_meshInput->vecFaceArea );
00584         }
00585 
00586         vecSourceFaceArea.Allocate( m_meshInput->faces.size() * m_nDofsPEl_Src * m_nDofsPEl_Src );
00587         int offset = 0;
00588         for( size_t e = 0; e < m_meshInput->faces.size(); e++ )
00589         {
00590             for( int s = 0; s < m_nDofsPEl_Src; s++ )
00591             {
00592                 for( int t = 0; t < m_nDofsPEl_Src; t++ )
00593                 {
00594                     vecSourceFaceArea[srccol_dtoc_dofmap[offset + s * m_nDofsPEl_Src + t]] =
00595                         dataGLLJacobianSrc[s][t][e];
00596                 }
00597             }
00598             offset += m_nDofsPEl_Src * m_nDofsPEl_Src;
00599         }
00600     }
00601 
00602     if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
00603     {
00604         this->InitializeCoordinatesFromMeshFV( *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon,
00605                                                dTargetVertexLat, false /* fLatLon = false */,
00606                                                m_remapper->max_target_edges );
00607 
00608         vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
00609         for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
00610             vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
00611     }
00612     else
00613     {
00614         DataArray3D< double > dataGLLJacobianDest;
00615         this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
00616                                                dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
00617 
00618         // Generate the continuous Jacobian for input mesh
00619         GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
00620 
00621         if( m_destDiscType == DiscretizationType_CGLL )
00622         {
00623             GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, m_meshOutput->vecFaceArea );
00624         }
00625         else
00626         {
00627             GenerateDiscontinuousJacobian( dataGLLJacobianDest, m_meshOutput->vecFaceArea );
00628         }
00629 
00630         vecTargetFaceArea.Allocate( m_meshOutput->faces.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest );
00631         int offset = 0;
00632         for( size_t e = 0; e < m_meshOutput->faces.size(); e++ )
00633         {
00634             for( int s = 0; s < m_nDofsPEl_Dest; s++ )
00635             {
00636                 for( int t = 0; t < m_nDofsPEl_Dest; t++ )
00637                 {
00638                     vecTargetFaceArea[row_dtoc_dofmap[offset + s * m_nDofsPEl_Dest + t]] = dataGLLJacobianDest[s][t][e];
00639                 }
00640             }
00641             offset += m_nDofsPEl_Dest * m_nDofsPEl_Dest;
00642         }
00643     }
00644 
00645     moab::EntityHandle& m_meshOverlapSet = m_remapper->m_overlap_set;
00646     int tot_src_ents                     = m_remapper->m_source_entities.size();
00647     int tot_tgt_ents                     = m_remapper->m_target_entities.size();
00648     int tot_src_size                     = dSourceCenterLon.GetRows();
00649     int tot_tgt_size                     = m_dTargetCenterLon.GetRows();
00650     int tot_vsrc_size                    = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
00651     int tot_vtgt_size                    = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
00652 
00653     const int weightMatNNZ = m_weightMatrix.nonZeros();
00654     moab::Tag tagMapMetaData, tagMapIndexRow, tagMapIndexCol, tagMapValues, srcEleIDs, tgtEleIDs;
00655     rval = m_interface->tag_get_handle( "SMAT_DATA", 13, moab::MB_TYPE_INTEGER, tagMapMetaData,
00656                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00657     rval = m_interface->tag_get_handle( "SMAT_ROWS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexRow,
00658                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00659     rval = m_interface->tag_get_handle( "SMAT_COLS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexCol,
00660                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00661     rval = m_interface->tag_get_handle( "SMAT_VALS", weightMatNNZ, moab::MB_TYPE_DOUBLE, tagMapValues,
00662                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00663     rval = m_interface->tag_get_handle( "SourceGIDS", tot_src_size, moab::MB_TYPE_INTEGER, srcEleIDs,
00664                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00665     rval = m_interface->tag_get_handle( "TargetGIDS", tot_tgt_size, moab::MB_TYPE_INTEGER, tgtEleIDs,
00666                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00667     moab::Tag srcAreaValues, tgtAreaValues;
00668     rval = m_interface->tag_get_handle( "SourceAreas", tot_src_size, moab::MB_TYPE_DOUBLE, srcAreaValues,
00669                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00670     rval = m_interface->tag_get_handle( "TargetAreas", tot_tgt_size, moab::MB_TYPE_DOUBLE, tgtAreaValues,
00671                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00672     moab::Tag tagSrcCoordsCLon, tagSrcCoordsCLat, tagTgtCoordsCLon, tagTgtCoordsCLat;
00673     rval = m_interface->tag_get_handle( "SourceCoordCenterLon", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLon,
00674                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00675     rval = m_interface->tag_get_handle( "SourceCoordCenterLat", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLat,
00676                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00677     rval = m_interface->tag_get_handle( "TargetCoordCenterLon", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLon,
00678                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00679     rval = m_interface->tag_get_handle( "TargetCoordCenterLat", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLat,
00680                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00681     moab::Tag tagSrcCoordsVLon, tagSrcCoordsVLat, tagTgtCoordsVLon, tagTgtCoordsVLat;
00682     rval = m_interface->tag_get_handle( "SourceCoordVertexLon", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLon,
00683                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00684     rval = m_interface->tag_get_handle( "SourceCoordVertexLat", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLat,
00685                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00686     rval = m_interface->tag_get_handle( "TargetCoordVertexLon", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLon,
00687                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00688     rval = m_interface->tag_get_handle( "TargetCoordVertexLat", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLat,
00689                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00690     moab::Tag srcMaskValues, tgtMaskValues;
00691     if( m_iSourceMask.IsAttached() )
00692     {
00693         rval = m_interface->tag_get_handle( "SourceMask", m_iSourceMask.GetRows(), moab::MB_TYPE_INTEGER, srcMaskValues,
00694                                             moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00695     }
00696     if( m_iTargetMask.IsAttached() )
00697     {
00698         rval = m_interface->tag_get_handle( "TargetMask", m_iTargetMask.GetRows(), moab::MB_TYPE_INTEGER, tgtMaskValues,
00699                                             moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
00700     }
00701 
00702     std::vector< int > smatrowvals( weightMatNNZ ), smatcolvals( weightMatNNZ );
00703     std::vector< double > smatvals( weightMatNNZ );
00704     // const double* smatvals = m_weightMatrix.valuePtr();
00705     // Loop over the matrix entries and find the max global ID for rows and columns
00706     for( int k = 0, offset = 0; k < m_weightMatrix.outerSize(); ++k )
00707     {
00708         for( moab::TempestOnlineMap::WeightMatrix::InnerIterator it( m_weightMatrix, k ); it; ++it, ++offset )
00709         {
00710             smatrowvals[offset] = this->GetRowGlobalDoF( it.row() );
00711             smatcolvals[offset] = this->GetColGlobalDoF( it.col() );
00712             smatvals[offset]    = it.value();
00713         }
00714     }
00715 
00716     /* Set the global IDs for the DoFs */
00717     ////
00718     // col_gdofmap [ col_ldofmap [ 0 : local_ndofs ] ] = GDOF
00719     // row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
00720     ////
00721     int maxrow = 0, maxcol = 0;
00722     std::vector< int > src_global_dofs( tot_src_size ), tgt_global_dofs( tot_tgt_size );
00723     for( int i = 0; i < tot_src_size; ++i )
00724     {
00725         src_global_dofs[i] = srccol_gdofmap[i];
00726         maxcol             = ( src_global_dofs[i] > maxcol ) ? src_global_dofs[i] : maxcol;
00727     }
00728 
00729     for( int i = 0; i < tot_tgt_size; ++i )
00730     {
00731         tgt_global_dofs[i] = row_gdofmap[i];
00732         maxrow             = ( tgt_global_dofs[i] > maxrow ) ? tgt_global_dofs[i] : maxrow;
00733     }
00734 
00735     ///////////////////////////////////////////////////////////////////////////
00736     // The metadata in H5M file contains the following data:
00737     //
00738     //   1. n_a: Total source entities: (number of elements in source mesh)
00739     //   2. n_b: Total target entities: (number of elements in target mesh)
00740     //   3. nv_a: Max edge size of elements in source mesh
00741     //   4. nv_b: Max edge size of elements in target mesh
00742     //   5. maxrows: Number of rows in remap weight matrix
00743     //   6. maxcols: Number of cols in remap weight matrix
00744     //   7. nnz: Number of total nnz in sparse remap weight matrix
00745     //   8. np_a: The order of the field description on the source mesh: >= 1
00746     //   9. np_b: The order of the field description on the target mesh: >= 1
00747     //   10. method_a: The type of discretization for field on source mesh: [0 = FV, 1 = cGLL, 2 =
00748     //   dGLL]
00749     //   11. method_b: The type of discretization for field on target mesh: [0 = FV, 1 = cGLL, 2 =
00750     //   dGLL]
00751     //   12. conserved: Flag to specify whether the remap operator has conservation constraints: [0,
00752     //   1]
00753     //   13. monotonicity: Flags to specify whether the remap operator has monotonicity constraints:
00754     //   [0, 1, 2]
00755     //
00756     ///////////////////////////////////////////////////////////////////////////
00757     int map_disc_details[6];
00758     map_disc_details[0] = m_nDofsPEl_Src;
00759     map_disc_details[1] = m_nDofsPEl_Dest;
00760     map_disc_details[2] = ( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD
00761                                 ? 0
00762                                 : ( m_srcDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
00763     map_disc_details[3] = ( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD
00764                                 ? 0
00765                                 : ( m_destDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
00766     map_disc_details[4] = ( m_bConserved ? 1 : 0 );
00767     map_disc_details[5] = m_iMonotonicity;
00768 
00769 #ifdef MOAB_HAVE_MPI
00770     int loc_smatmetadata[13] = { tot_src_ents,
00771                                  tot_tgt_ents,
00772                                  m_remapper->max_source_edges,
00773                                  m_remapper->max_target_edges,
00774                                  maxrow + 1,
00775                                  maxcol + 1,
00776                                  weightMatNNZ,
00777                                  map_disc_details[0],
00778                                  map_disc_details[1],
00779                                  map_disc_details[2],
00780                                  map_disc_details[3],
00781                                  map_disc_details[4],
00782                                  map_disc_details[5] };
00783     rval                     = m_interface->tag_set_data( tagMapMetaData, &m_meshOverlapSet, 1, &loc_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00784     int glb_smatmetadata[13] = { 0,
00785                                  0,
00786                                  0,
00787                                  0,
00788                                  0,
00789                                  0,
00790                                  0,
00791                                  map_disc_details[0],
00792                                  map_disc_details[1],
00793                                  map_disc_details[2],
00794                                  map_disc_details[3],
00795                                  map_disc_details[4],
00796                                  map_disc_details[5] };
00797     int loc_buf[7]           = {
00798         tot_src_ents, tot_tgt_ents, weightMatNNZ, m_remapper->max_source_edges, m_remapper->max_target_edges,
00799         maxrow,       maxcol
00800     };
00801     int glb_buf[4] = { 0, 0, 0, 0 };
00802     MPI_Reduce( &loc_buf[0], &glb_buf[0], 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
00803     glb_smatmetadata[0] = glb_buf[0];
00804     glb_smatmetadata[1] = glb_buf[1];
00805     glb_smatmetadata[6] = glb_buf[2];
00806     MPI_Reduce( &loc_buf[3], &glb_buf[0], 4, MPI_INT, MPI_MAX, 0, m_pcomm->comm() );
00807     glb_smatmetadata[2] = glb_buf[0];
00808     glb_smatmetadata[3] = glb_buf[1];
00809     glb_smatmetadata[4] = glb_buf[2];
00810     glb_smatmetadata[5] = glb_buf[3];
00811 #else
00812     int glb_smatmetadata[13] = { tot_src_ents,
00813                                  tot_tgt_ents,
00814                                  m_remapper->max_source_edges,
00815                                  m_remapper->max_target_edges,
00816                                  maxrow,
00817                                  maxcol,
00818                                  weightMatNNZ,
00819                                  map_disc_details[0],
00820                                  map_disc_details[1],
00821                                  map_disc_details[2],
00822                                  map_disc_details[3],
00823                                  map_disc_details[4],
00824                                  map_disc_details[5] };
00825 #endif
00826     // These values represent number of rows and columns. So should be 1-based.
00827     glb_smatmetadata[4]++;
00828     glb_smatmetadata[5]++;
00829 
00830     if( this->is_root )
00831     {
00832         std::cout << "  " << this->rank << "  Writing remap weights with size [" << glb_smatmetadata[4] << " X "
00833                   << glb_smatmetadata[5] << "] and NNZ = " << glb_smatmetadata[6] << std::endl;
00834         EntityHandle root_set = 0;
00835         rval                  = m_interface->tag_set_data( tagMapMetaData, &root_set, 1, &glb_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00836     }
00837 
00838     int dsize;
00839     const int numval          = weightMatNNZ;
00840     const void* smatrowvals_d = smatrowvals.data();
00841     const void* smatcolvals_d = smatcolvals.data();
00842     const void* smatvals_d    = smatvals.data();
00843     rval = m_interface->tag_set_by_ptr( tagMapIndexRow, &m_meshOverlapSet, 1, &smatrowvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00844     rval = m_interface->tag_set_by_ptr( tagMapIndexCol, &m_meshOverlapSet, 1, &smatcolvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00845     rval = m_interface->tag_set_by_ptr( tagMapValues, &m_meshOverlapSet, 1, &smatvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00846 
00847     /* Set the global IDs for the DoFs */
00848     const void* srceleidvals_d = src_global_dofs.data();
00849     const void* tgteleidvals_d = tgt_global_dofs.data();
00850     dsize                      = src_global_dofs.size();
00851     rval = m_interface->tag_set_by_ptr( srcEleIDs, &m_meshOverlapSet, 1, &srceleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00852     dsize = tgt_global_dofs.size();
00853     rval  = m_interface->tag_set_by_ptr( tgtEleIDs, &m_meshOverlapSet, 1, &tgteleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00854 
00855     /* Set the source and target areas */
00856     const void* srcareavals_d = vecSourceFaceArea;
00857     const void* tgtareavals_d = vecTargetFaceArea;
00858     dsize                     = tot_src_size;
00859     rval = m_interface->tag_set_by_ptr( srcAreaValues, &m_meshOverlapSet, 1, &srcareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00860     dsize = tot_tgt_size;
00861     rval  = m_interface->tag_set_by_ptr( tgtAreaValues, &m_meshOverlapSet, 1, &tgtareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00862 
00863     /* Set the coordinates for source and target center vertices */
00864     const void* srccoordsclonvals_d = &dSourceCenterLon[0];
00865     const void* srccoordsclatvals_d = &dSourceCenterLat[0];
00866     dsize                           = dSourceCenterLon.GetRows();
00867     rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLon, &m_meshOverlapSet, 1, &srccoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00868     rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLat, &m_meshOverlapSet, 1, &srccoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00869     const void* tgtcoordsclonvals_d = &m_dTargetCenterLon[0];
00870     const void* tgtcoordsclatvals_d = &m_dTargetCenterLat[0];
00871     dsize                           = vecTargetFaceArea.GetRows();
00872     rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLon, &m_meshOverlapSet, 1, &tgtcoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00873     rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLat, &m_meshOverlapSet, 1, &tgtcoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00874 
00875     /* Set the coordinates for source and target element vertices */
00876     const void* srccoordsvlonvals_d = &( dSourceVertexLon[0][0] );
00877     const void* srccoordsvlatvals_d = &( dSourceVertexLat[0][0] );
00878     dsize                           = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
00879     rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLon, &m_meshOverlapSet, 1, &srccoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00880     rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLat, &m_meshOverlapSet, 1, &srccoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00881     const void* tgtcoordsvlonvals_d = &( m_dTargetVertexLon[0][0] );
00882     const void* tgtcoordsvlatvals_d = &( m_dTargetVertexLat[0][0] );
00883     dsize                           = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
00884     rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLon, &m_meshOverlapSet, 1, &tgtcoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00885     rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLat, &m_meshOverlapSet, 1, &tgtcoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00886 
00887     /* Set the masks for source and target meshes if available */
00888     if( m_iSourceMask.IsAttached() )
00889     {
00890         const void* srcmaskvals_d = m_iSourceMask;
00891         dsize                     = m_iSourceMask.GetRows();
00892         rval = m_interface->tag_set_by_ptr( srcMaskValues, &m_meshOverlapSet, 1, &srcmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00893     }
00894 
00895     if( m_iTargetMask.IsAttached() )
00896     {
00897         const void* tgtmaskvals_d = m_iTargetMask;
00898         dsize                     = m_iTargetMask.GetRows();
00899         rval = m_interface->tag_set_by_ptr( tgtMaskValues, &m_meshOverlapSet, 1, &tgtmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
00900     }
00901 
00902 #ifdef MOAB_HAVE_MPI
00903     const char* writeOptions = ( this->size > 1 ? "PARALLEL=WRITE_PART" : "" );
00904 #else
00905     const char* writeOptions = "";
00906 #endif
00907 
00908     // EntityHandle sets[3] = {m_remapper->m_source_set, m_remapper->m_target_set, m_remapper->m_overlap_set};
00909     EntityHandle sets[1] = { m_remapper->m_overlap_set };
00910     rval                 = m_interface->write_file( strOutputFile.c_str(), NULL, writeOptions, sets, 1 );MB_CHK_ERR( rval );
00911 
00912 #ifdef WRITE_SCRIP_FILE
00913     sstr.str( "" );
00914     sstr << ctx.outFilename.substr( 0, lastindex ) << "_" << proc_id << ".nc";
00915     std::map< std::string, std::string > mapAttributes;
00916     mapAttributes["Creator"] = "MOAB mbtempest workflow";
00917     if( !ctx.proc_id ) std::cout << "Writing offline map to file: " << sstr.str() << std::endl;
00918     this->Write( strOutputFile.c_str(), mapAttributes, NcFile::Netcdf4 );
00919     sstr.str( "" );
00920 #endif
00921 
00922     return moab::MB_SUCCESS;
00923 }
00924 
00925 ///////////////////////////////////////////////////////////////////////////////
00926 
00927 static void print_progress( const int barWidth, const float progress, const char* message )
00928 {
00929     std::cout << message << " [";
00930     int pos = barWidth * progress;
00931     for( int i = 0; i < barWidth; ++i )
00932     {
00933         if( i < pos )
00934             std::cout << "=";
00935         else if( i == pos )
00936             std::cout << ">";
00937         else
00938             std::cout << " ";
00939     }
00940     std::cout << "] " << int( progress * 100.0 ) << " %\r";
00941     std::cout.flush();
00942 }
00943 
00944 ///////////////////////////////////////////////////////////////////////////////
00945 
00946 moab::ErrorCode moab::TempestOnlineMap::ReadParallelMap( const char* strSource, const std::vector< int >& owned_dof_ids,
00947                                                          bool /* row_major_ownership */ )
00948 {
00949     NcError error( NcError::silent_nonfatal );
00950 
00951     NcVar *varRow = NULL, *varCol = NULL, *varS = NULL;
00952     int nS = 0, nB = 0;
00953 #ifdef MOAB_HAVE_NETCDFPAR
00954     bool is_independent = true;
00955     ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
00956     // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
00957 #else
00958     NcFile ncMap( strSource, NcFile::ReadOnly );
00959 #endif
00960 
00961     // Read SparseMatrix entries
00962 
00963     NcDim* dimNS = ncMap.get_dim( "n_s" );
00964     if( dimNS == NULL ) { _EXCEPTION1( "Map file \"%s\" does not contain dimension \"n_s\"", strSource ); }
00965 
00966     NcDim* dimNB = ncMap.get_dim( "n_b" );
00967     if( dimNB == NULL ) { _EXCEPTION1( "Map file \"%s\" does not contain dimension \"nB\"", strSource ); }
00968 
00969     // store total number of nonzeros
00970     nS = dimNS->size();
00971     nB = dimNB->size();
00972 
00973     varRow = ncMap.get_var( "row" );
00974     if( varRow == NULL ) { _EXCEPTION1( "Map file \"%s\" does not contain variable \"row\"", strSource ); }
00975 
00976     varCol = ncMap.get_var( "col" );
00977     if( varCol == NULL ) { _EXCEPTION1( "Map file \"%s\" does not contain variable \"col\"", strSource ); }
00978 
00979     varS = ncMap.get_var( "S" );
00980     if( varS == NULL ) { _EXCEPTION1( "Map file \"%s\" does not contain variable \"S\"", strSource ); }
00981 
00982 #ifdef MOAB_HAVE_NETCDFPAR
00983     ncMap.enable_var_par_access( varRow, is_independent );
00984     ncMap.enable_var_par_access( varCol, is_independent );
00985     ncMap.enable_var_par_access( varS, is_independent );
00986 #endif
00987 
00988     std::vector< int > rowOwnership;
00989     // if owned_dof_ids = NULL, use the default trivial partitioning scheme
00990     if( owned_dof_ids.size() == 0 )
00991     {
00992         // assert(row_major_ownership == true); // this block is valid only for row-based partitioning
00993         rowOwnership.resize( size );
00994         int nGRowPerPart   = nB / size;
00995         int nGRowRemainder = nB % size;  // Keep the remainder in root
00996         rowOwnership[0]    = nGRowPerPart + nGRowRemainder;
00997         for( int ip = 1, roffset = rowOwnership[0]; ip < size; ++ip )
00998         {
00999             roffset += nGRowPerPart;
01000             rowOwnership[ip] = roffset;
01001         }
01002     }
01003 
01004     // Let us declare the map object for every process
01005     SparseMatrix< double >& sparseMatrix = this->GetSparseMatrix();
01006 
01007     std::map< int, int > rowMap, colMap;
01008     int rindexMax = 0, cindexMax = 0;
01009 
01010     int localSize   = nS / size;
01011     long offsetRead = rank * localSize;
01012     // leftovers on last rank
01013     if( rank == size - 1 ) { localSize += nS % size; }
01014 
01015     std::vector< int > vecRow, vecCol;
01016     std::vector< double > vecS;
01017     vecRow.resize( localSize );
01018     vecCol.resize( localSize );
01019     vecS.resize( localSize );
01020 
01021     varRow->set_cur( (long)( offsetRead ) );
01022     varRow->get( &( vecRow[0] ), localSize );
01023 
01024     varCol->set_cur( (long)( offsetRead ) );
01025     varCol->get( &( vecCol[0] ), localSize );
01026 
01027     varS->set_cur( (long)( offsetRead ) );
01028     varS->get( &( vecS[0] ), localSize );
01029 
01030     ncMap.close();
01031 
01032     // bother with tuple list only if size > 1
01033     // otherwise, just fill the sparse matrix
01034 #ifdef MOAB_HAVE_MPI
01035     if( size > 1 )
01036     {
01037         // send to
01038         moab::TupleList tl;
01039         unsigned numr = 1;                          //
01040         tl.initialize( 3, 0, 0, numr, localSize );  // to proc, row, col, value
01041         tl.enableWriteAccess();
01042         // populate
01043         for( int i = 0; i < localSize; i++ )
01044         {
01045             int rowval  = vecRow[i] - 1;  // dofs are 1 based in the file
01046             int colval  = vecCol[i] - 1;
01047             int to_proc = -1;
01048             //
01049 
01050             if( rowOwnership[0] > rowval )
01051                 to_proc = 0;
01052             else
01053             {
01054                 for( int ip = 1; ip < size; ++ip )
01055                 {
01056                     if( rowOwnership[ip - 1] <= rowval && rowOwnership[ip] > rowval )
01057                     {
01058                         to_proc = ip;
01059                         break;
01060                     }
01061                 }
01062             }
01063 
01064             int n               = tl.get_n();
01065             tl.vi_wr[3 * n]     = to_proc;
01066             tl.vi_wr[3 * n + 1] = rowval;
01067             tl.vi_wr[3 * n + 2] = colval;
01068             tl.vr_wr[n]         = vecS[i];
01069 
01070             tl.inc_n();
01071         }
01072 
01073 
01074         // now do the heavy communication
01075         ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tl, 0 );
01076         // populate the sparsematrix, using rowMap and colMap; what is the need for them?
01077         int n = tl.get_n();
01078         for( int i = 0; i < n; i++ )
01079         {
01080 
01081             int rindex, cindex;
01082             const int& vecRowValue = tl.vi_wr[3 * i + 1];
01083             const int& vecColValue = tl.vi_wr[3 * i + 2];
01084 
01085             std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
01086             if( riter == rowMap.end() )
01087             {
01088                 rowMap[vecRowValue] = rindexMax;
01089                 rindex              = rindexMax;
01090                 rindexMax++;
01091             }
01092             else
01093                 rindex = riter->second;
01094 
01095             std::map< int, int >::iterator citer = colMap.find( vecColValue );
01096             if( citer == colMap.end() )
01097             {
01098                 colMap[vecColValue] = cindexMax;
01099                 cindex              = cindexMax;
01100                 cindexMax++;
01101             }
01102             else
01103                 cindex = citer->second;
01104 
01105             sparseMatrix( rindex, cindex ) = tl.vr_wr[i];
01106         }
01107     }
01108     else
01109 #endif
01110     {
01111         for( int i = 0; i < nS; i++ )
01112         {
01113 
01114             int rindex, cindex;
01115             const int& vecRowValue = vecRow[i];
01116             const int& vecColValue = vecCol[i];
01117 
01118             std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
01119             if( riter == rowMap.end() )
01120             {
01121                 rowMap[vecRowValue] = rindexMax;
01122                 rindex              = rindexMax;
01123                 rindexMax++;
01124             }
01125             else
01126                 rindex = riter->second;
01127 
01128             std::map< int, int >::iterator citer = colMap.find( vecColValue );
01129             if( citer == colMap.end() )
01130             {
01131                 colMap[vecColValue] = cindexMax;
01132                 cindex              = cindexMax;
01133                 cindexMax++;
01134             }
01135             else
01136                 cindex = citer->second;
01137 
01138             sparseMatrix( rindex, cindex ) = vecS[i];
01139         }
01140     }
01141 
01142     m_nTotDofs_SrcCov = sparseMatrix.GetColumns();
01143     m_nTotDofs_Dest   = sparseMatrix.GetRows();
01144 
01145 #ifdef MOAB_HAVE_EIGEN3
01146     this->copy_tempest_sparsemat_to_eigen3();
01147 #endif
01148 
01149     return moab::MB_SUCCESS;
01150 }
01151 
01152 ///////////////////////////////////////////////////////////////////////////////
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines