MOAB: Mesh Oriented datABase  (version 5.4.0)
TempestOnlineMap.hpp
Go to the documentation of this file.
00001 ///////////////////////////////////////////////////////////////////////////////
00002 ///
00003 /// \file    TempestOnlineMap.h
00004 /// \author  Vijay Mahadevan
00005 /// \version November 20, 2017
00006 ///
00007 
00008 #ifndef _TEMPESTONLINEMAP_H_
00009 #define _TEMPESTONLINEMAP_H_
00010 
00011 #include "moab/MOABConfig.h"
00012 
00013 // Tempest includes
00014 #ifdef MOAB_HAVE_TEMPESTREMAP
00015 #include "moab/Remapping/TempestRemapper.hpp"
00016 #include "OfflineMap.h"
00017 #else
00018 #error Re-configure with TempestRemap
00019 #endif
00020 
00021 #include <string>
00022 #include <vector>
00023 
00024 #ifdef MOAB_HAVE_EIGEN3
00025 #include <Eigen/Sparse>
00026 #endif
00027 
00028 ///////////////////////////////////////////////////////////////////////////////
00029 
00030 #if !defined( RECTANGULAR_TRUNCATION ) && !defined( TRIANGULAR_TRUNCATION )
00031 #define RECTANGULAR_TRUNCATION
00032 // #define TRIANGULAR_TRUNCATION
00033 #endif
00034 
00035 ///////////////////////////////////////////////////////////////////////////////
00036 
00037 // Forward declarations
00038 class Mesh;
00039 
00040 ///////////////////////////////////////////////////////////////////////////////
00041 
00042 namespace moab
00043 {
00044 
00045 /// <summary>
00046 ///     An offline map between two Meshes.
00047 /// </summary>
00048 class TempestOnlineMap : public OfflineMap
00049 {
00050 
00051   public:
00052     /// <summary>
00053     ///     Generate the metadata associated with the offline map.
00054     /// </summary>
00055     TempestOnlineMap( moab::TempestRemapper* remapper );
00056 
00057     /// <summary>
00058     ///     Define a virtual destructor.
00059     /// </summary>
00060     virtual ~TempestOnlineMap();
00061 
00062   public:
00063     // Input / Output types
00064     enum DiscretizationType
00065     {
00066         DiscretizationType_FV,
00067         DiscretizationType_CGLL,
00068         DiscretizationType_DGLL,
00069         DiscretizationType_PCLOUD
00070     };
00071 
00072     /// <summary>
00073     ///     Generate the offline map, given the source and target mesh and discretization details.
00074     ///     This method generates the mapping between the two meshes based on the overlap and stores
00075     ///     the result in the SparseMatrix.
00076     /// </summary>
00077     moab::ErrorCode GenerateRemappingWeights( std::string strInputType,
00078                                               std::string strOutputType,
00079                                               const GenerateOfflineMapAlgorithmOptions& mapOptions,
00080                                               const std::string& srcDofTagName = "GLOBAL_ID",
00081                                               const std::string& tgtDofTagName = "GLOBAL_ID" );
00082 
00083     /// <summary>
00084     ///     Generate the metadata associated with the offline map.
00085     /// </summary>
00086     // moab::ErrorCode GenerateMetaData();
00087 
00088     /// <summary>
00089     ///     Read the OfflineMap from a NetCDF file.
00090     /// </summary>
00091     moab::ErrorCode ReadParallelMap( const char* strSource,
00092                                      const std::vector< int >& owned_dof_ids,
00093                                      bool row_major_ownership = true );
00094 
00095     /// <summary>
00096     ///     Write the TempestOnlineMap to a parallel NetCDF file.
00097     /// </summary>
00098     moab::ErrorCode WriteParallelMap( const std::string& strTarget );
00099 
00100     /// <summary>
00101     ///     Determine if the map is first-order accurate.
00102     /// </summary>
00103     virtual int IsConsistent( double dTolerance );
00104 
00105     /// <summary>
00106     ///     Determine if the map is conservative.
00107     /// </summary>
00108     virtual int IsConservative( double dTolerance );
00109 
00110     /// <summary>
00111     ///     Determine if the map is monotone.
00112     /// </summary>
00113     virtual int IsMonotone( double dTolerance );
00114 
00115     /// <summary>
00116     ///     If we computed the reduction, get the vector representing the source areas for all entities
00117     /// in the mesh
00118     /// </summary>
00119     const DataArray1D< double >& GetGlobalSourceAreas() const;
00120 
00121     /// <summary>
00122     ///     If we computed the reduction, get the vector representing the target areas for all entities
00123     /// in the mesh
00124     /// </summary>
00125     const DataArray1D< double >& GetGlobalTargetAreas() const;
00126 
00127   private:
00128     /// <summary>
00129     ///     Compute the remapping weights as a permutation matrix that relates DoFs on the source
00130     /// mesh
00131     ///     to DoFs on the target mesh.
00132     /// </summary>
00133     moab::ErrorCode LinearRemapNN_MOAB( bool use_GID_matching = false, bool strict_check = false );
00134 
00135     /// <summary>
00136     ///     Compute the remapping weights for a FV field defined on the source to a
00137     ///     FV field defined on the target mesh.
00138     /// </summary>
00139     void LinearRemapFVtoFV_Tempest_MOAB( int nOrder );
00140 
00141     /// <summary>
00142     ///     Generate the OfflineMap for linear conserative element-average
00143     ///     spectral element to element average remapping.
00144     /// </summary>
00145     void LinearRemapSE0_Tempest_MOAB( const DataArray3D< int >& dataGLLNodes,
00146                                       const DataArray3D< double >& dataGLLJacobian );
00147 
00148     /// <summary>
00149     ///     Generate the OfflineMap for cubic conserative element-average
00150     ///     spectral element to element average remapping.
00151     /// </summary>
00152     void LinearRemapSE4_Tempest_MOAB( const DataArray3D< int >& dataGLLNodes,
00153                                       const DataArray3D< double >& dataGLLJacobian,
00154                                       int nMonotoneType,
00155                                       bool fContinuousIn,
00156                                       bool fNoConservation );
00157 
00158     /// <summary>
00159     ///     Generate the OfflineMap for remapping from finite volumes to finite
00160     ///     elements.
00161     /// </summary>
00162     void LinearRemapFVtoGLL_MOAB( const DataArray3D< int >& dataGLLNodes,
00163                                   const DataArray3D< double >& dataGLLJacobian,
00164                                   const DataArray1D< double >& dataGLLNodalArea,
00165                                   int nOrder,
00166                                   int nMonotoneType,
00167                                   bool fContinuous,
00168                                   bool fNoConservation );
00169 
00170     /// <summary>
00171     ///     Generate the OfflineMap for remapping from finite elements to finite
00172     ///     elements.
00173     /// </summary>
00174     void LinearRemapGLLtoGLL2_MOAB( const DataArray3D< int >& dataGLLNodesIn,
00175                                     const DataArray3D< double >& dataGLLJacobianIn,
00176                                     const DataArray3D< int >& dataGLLNodesOut,
00177                                     const DataArray3D< double >& dataGLLJacobianOut,
00178                                     const DataArray1D< double >& dataNodalAreaOut,
00179                                     int nPin,
00180                                     int nPout,
00181                                     int nMonotoneType,
00182                                     bool fContinuousIn,
00183                                     bool fContinuousOut,
00184                                     bool fNoConservation );
00185 
00186     /// <summary>
00187     ///     Generate the OfflineMap for remapping from finite elements to finite
00188     ///     elements (pointwise interpolation).
00189     /// </summary>
00190     void LinearRemapGLLtoGLL2_Pointwise_MOAB( const DataArray3D< int >& dataGLLNodesIn,
00191                                               const DataArray3D< double >& dataGLLJacobianIn,
00192                                               const DataArray3D< int >& dataGLLNodesOut,
00193                                               const DataArray3D< double >& dataGLLJacobianOut,
00194                                               const DataArray1D< double >& dataNodalAreaOut,
00195                                               int nPin,
00196                                               int nPout,
00197                                               int nMonotoneType,
00198                                               bool fContinuousIn,
00199                                               bool fContinuousOut );
00200 
00201     /// <summary>
00202     ///     Copy the local matrix from Tempest SparseMatrix representation (ELL)
00203     ///     to the parallel CSR Eigen Matrix for scalable application of matvec
00204     ///     needed for projections.
00205     /// </summary>
00206 #ifdef MOAB_HAVE_EIGEN3
00207     void copy_tempest_sparsemat_to_eigen3();
00208 #endif
00209 
00210     /// <summary>
00211     ///     Parallel I/O with HDF5 to write out the remapping weights from multiple processors.
00212     /// </summary>
00213     moab::ErrorCode WriteSCRIPMapFile( const std::string& strOutputFile );
00214 
00215     /// <summary>
00216     ///     Parallel I/O with NetCDF to write out the SCRIP file from multiple processors.
00217     /// </summary>
00218     moab::ErrorCode WriteHDF5MapFile( const std::string& filename );
00219 
00220   public:
00221     /// <summary>
00222     ///     Store the tag names associated with global DoF ids for source and target meshes
00223     /// </summary>
00224     moab::ErrorCode SetDOFmapTags( const std::string srcDofTagName, const std::string tgtDofTagName );
00225 
00226     /// <summary>
00227     ///     Compute the association between the solution tag global DoF numbering and
00228     ///     the local matrix numbering so that matvec operations can be performed
00229     ///     consistently.
00230     /// </summary>
00231     moab::ErrorCode SetDOFmapAssociation( DiscretizationType srcType,
00232                                           bool isSrcContinuous,
00233                                           DataArray3D< int >* srcdataGLLNodes,
00234                                           DataArray3D< int >* srcdataGLLNodesSrc,
00235                                           DiscretizationType destType,
00236                                           bool isDestContinuous,
00237                                           DataArray3D< int >* tgtdataGLLNodes );
00238 
00239 #ifdef MOAB_HAVE_EIGEN3
00240 
00241     typedef Eigen::Matrix< double, 1, Eigen::Dynamic > WeightDRowVector;
00242     typedef Eigen::Matrix< double, Eigen::Dynamic, 1 > WeightDColVector;
00243     typedef Eigen::SparseVector< double > WeightSVector;
00244     typedef Eigen::SparseMatrix< double, Eigen::RowMajor > WeightRMatrix;
00245     typedef Eigen::SparseMatrix< double, Eigen::ColMajor > WeightCMatrix;
00246 
00247     typedef WeightDRowVector WeightRowVector;
00248     typedef WeightDColVector WeightColVector;
00249     typedef WeightRMatrix WeightMatrix;
00250 
00251     /// <summary>
00252     ///     Get the raw reference to the Eigen weight matrix representing the projection from source to
00253     /// destination mesh.
00254     /// </summary>
00255     WeightMatrix& GetWeightMatrix();
00256 
00257     /// <summary>
00258     ///     Get the row vector that is amenable for application of A*x operation.
00259     /// </summary>
00260     WeightRowVector& GetRowVector();
00261 
00262     /// <summary>
00263     ///     Get the column vector that is amenable for application of A^T*x operation.
00264     /// </summary>
00265     WeightColVector& GetColVector();
00266 
00267 #endif
00268 
00269     /// <summary>
00270     ///     Get the number of total Degrees-Of-Freedom defined on the source mesh.
00271     /// </summary>
00272     int GetSourceGlobalNDofs();
00273 
00274     /// <summary>
00275     ///     Get the number of total Degrees-Of-Freedom defined on the destination mesh.
00276     /// </summary>
00277     int GetDestinationGlobalNDofs();
00278 
00279     /// <summary>
00280     ///     Get the number of local Degrees-Of-Freedom defined on the source mesh.
00281     /// </summary>
00282     int GetSourceLocalNDofs();
00283 
00284     /// <summary>
00285     ///     Get the number of local Degrees-Of-Freedom defined on the destination mesh.
00286     /// </summary>
00287     int GetDestinationLocalNDofs();
00288 
00289     /// <summary>
00290     ///     Get the number of Degrees-Of-Freedom per element on the source mesh.
00291     /// </summary>
00292     int GetSourceNDofsPerElement();
00293 
00294     /// <summary>
00295     ///     Get the number of Degrees-Of-Freedom per element on the destination mesh.
00296     /// </summary>
00297     int GetDestinationNDofsPerElement();
00298 
00299     /// <summary>
00300     ///     Set the number of Degrees-Of-Freedom per element on the source mesh.
00301     /// </summary>
00302     void SetSourceNDofsPerElement( int ns );
00303 
00304     /// <summary>
00305     ///     Get the number of Degrees-Of-Freedom per element on the destination mesh.
00306     /// </summary>
00307     void SetDestinationNDofsPerElement( int nt );
00308 
00309     /// <summary>
00310     ///     Get the global Degrees-Of-Freedom ID on the destination mesh.
00311     /// </summary>
00312     int GetRowGlobalDoF( int localID ) const;
00313 
00314     /// <summary>
00315     ///     Get the index of globaRowDoF.
00316     /// </summary>
00317     inline int GetIndexOfRowGlobalDoF( int globalRowDoF ) const;
00318 
00319     /// <summary>
00320     ///     Get the global Degrees-Of-Freedom ID on the source mesh.
00321     /// </summary>
00322     int GetColGlobalDoF( int localID ) const;
00323 
00324     /// <summary>
00325     ///     Get the index of globaColDoF.
00326     /// </summary>
00327     inline int GetIndexOfColGlobalDoF( int globalColDoF ) const;
00328 
00329     /// <summary>
00330     ///     Apply the weight matrix onto the source vector provided as input, and return the column
00331     /// vector (solution projection) after the map application
00332     ///     Compute:        \p tgtVals = A(S->T) * \srcVals, or
00333     ///     if (transpose)  \p tgtVals = [A(T->S)]^T * \srcVals
00334     /// </summary>
00335     moab::ErrorCode ApplyWeights( std::vector< double >& srcVals,
00336                                   std::vector< double >& tgtVals,
00337                                   bool transpose = false );
00338 
00339     /// <summary>
00340     ///     Apply the weight matrix onto the source vector (tag) provided as input, and return the
00341     /// column vector (solution projection) in a tag, after the map application
00342     ///     Compute:        \p tgtVals = A(S->T) * \srcVals, or
00343     ///     if (transpose)  \p tgtVals = [A(T->S)]^T * \srcVals
00344     /// </summary>
00345     moab::ErrorCode ApplyWeights( moab::Tag srcSolutionTag, moab::Tag tgtSolutionTag, bool transpose = false );
00346 
00347     typedef double ( *sample_function )( double, double );
00348 
00349     /// <summary>
00350     ///     Define an analytical solution over the given (source or target) mesh, as specificed in
00351     ///     the context. This routine will define a tag that is compatible with the specified
00352     ///     discretization method type and order and sample the solution exactly using the
00353     ///     analytical function provided by the user.
00354     /// </summary>
00355     moab::ErrorCode DefineAnalyticalSolution( moab::Tag& exactSolnTag,
00356                                               const std::string& solnName,
00357                                               Remapper::IntersectionContext ctx,
00358                                               sample_function testFunction,
00359                                               moab::Tag* clonedSolnTag  = NULL,
00360                                               std::string cloneSolnName = "" );
00361 
00362     /// <summary>
00363     ///     Compute the error between a sampled (exact) solution and a projected solution in various
00364     ///     error norms.
00365     /// </summary>
00366     moab::ErrorCode ComputeMetrics( Remapper::IntersectionContext ctx,
00367                                     moab::Tag& exactTag,
00368                                     moab::Tag& approxTag,
00369                                     std::map< std::string, double >& metrics,
00370                                     bool verbose = true );
00371 
00372     moab::ErrorCode fill_row_ids( std::vector< int >& ids_of_interest )
00373     {
00374         ids_of_interest.reserve( row_gdofmap.size() );
00375         // need to add 1
00376         for( auto it = row_gdofmap.begin(); it != row_gdofmap.end(); it++ )
00377             ids_of_interest.push_back( *it + 1 );
00378 
00379         return moab::MB_SUCCESS;
00380     }
00381 
00382     moab::ErrorCode fill_col_ids( std::vector< int >& ids_of_interest )
00383     {
00384         ids_of_interest.reserve( col_gdofmap.size() );
00385         // need to add 1
00386         for( auto it = col_gdofmap.begin(); it != col_gdofmap.end(); it++ )
00387             ids_of_interest.push_back( *it + 1 );
00388         return moab::MB_SUCCESS;
00389     }
00390 
00391     moab::ErrorCode set_col_dc_dofs( std::vector< int >& values_entities );
00392 
00393     moab::ErrorCode set_row_dc_dofs( std::vector< int >& values_entities );
00394 
00395   private:
00396     void setup_sizes_dimensions();
00397 
00398 #ifdef MOAB_HAVE_MPI
00399     int rearrange_arrays_by_dofs( const std::vector< unsigned int >& gdofmap,
00400                                   DataArray1D< double >& vecFaceArea,
00401                                   DataArray1D< double >& dCenterLon,
00402                                   DataArray1D< double >& dCenterLat,
00403                                   DataArray2D< double >& dVertexLat,
00404                                   DataArray2D< double >& dVertexLon,
00405                                   std::vector< int >& masks,
00406                                   unsigned& N,  // this will be output too now
00407                                   int nv,
00408                                   int& maxdof );
00409 #endif
00410 
00411     /// <summary>
00412     ///     The fundamental remapping operator object.
00413     /// </summary>
00414     moab::TempestRemapper* m_remapper;
00415 
00416 #ifdef MOAB_HAVE_EIGEN3
00417 
00418     WeightMatrix m_weightMatrix;
00419     WeightRowVector m_rowVector;
00420     WeightColVector m_colVector;
00421 
00422 #endif
00423 
00424     /// <summary>
00425     ///     The reference to the moab::Core object that contains source/target and overlap sets.
00426     /// </summary>
00427     moab::Interface* m_interface;
00428 
00429 #ifdef MOAB_HAVE_MPI
00430     /// <summary>
00431     ///     The reference to the parallel communicator object used by the Core object.
00432     /// </summary>
00433     moab::ParallelComm* m_pcomm;
00434 #endif
00435 
00436     /// <summary>
00437     ///     The original tag data and local to global DoF mapping to associate matrix values to
00438     /// solution    <summary>
00439     moab::Tag m_dofTagSrc, m_dofTagDest;
00440     std::vector< unsigned > row_gdofmap, col_gdofmap, srccol_gdofmap;
00441 
00442     // make it int, because it can be -1 in new logic
00443     std::vector< int > row_dtoc_dofmap, col_dtoc_dofmap, srccol_dtoc_dofmap;
00444 
00445     std::map< int, int > rowMap, colMap;
00446 
00447     DataArray3D< int > dataGLLNodesSrc, dataGLLNodesSrcCov, dataGLLNodesDest;
00448     DiscretizationType m_srcDiscType, m_destDiscType;
00449     int m_nTotDofs_Src, m_nTotDofs_SrcCov, m_nTotDofs_Dest;
00450 
00451     // Key details about the current map
00452     int m_nDofsPEl_Src, m_nDofsPEl_Dest;
00453     DiscretizationType m_eInputType, m_eOutputType;
00454     bool m_bConserved;
00455     int m_iMonotonicity;
00456 
00457     Mesh* m_meshInput;
00458     Mesh* m_meshInputCov;
00459     Mesh* m_meshOutput;
00460     Mesh* m_meshOverlap;
00461 
00462     bool is_parallel, is_root;
00463     int rank, size, root_proc;
00464 };
00465 
00466 ///////////////////////////////////////////////////////////////////////////////
00467 
00468 inline int moab::TempestOnlineMap::GetRowGlobalDoF( int localRowID ) const
00469 {
00470     return row_gdofmap[localRowID];
00471 }
00472 
00473 inline int moab::TempestOnlineMap::GetIndexOfRowGlobalDoF( int globalRowDoF ) const /* 0 based */
00474 {
00475     return globalRowDoF + 1;
00476 }
00477 ///////////////////////////////////////////////////////////////////////////////
00478 
00479 inline int moab::TempestOnlineMap::GetColGlobalDoF( int localColID ) const
00480 {
00481     return col_gdofmap[localColID];
00482 }
00483 
00484 inline int moab::TempestOnlineMap::GetIndexOfColGlobalDoF( int globalColDoF ) const /* 0 based */
00485 {
00486     return globalColDoF + 1;  // temporary
00487 }
00488 ///////////////////////////////////////////////////////////////////////////////
00489 
00490 inline int moab::TempestOnlineMap::GetSourceNDofsPerElement()
00491 {
00492     return m_nDofsPEl_Src;
00493 }
00494 
00495 ///////////////////////////////////////////////////////////////////////////////
00496 
00497 inline int moab::TempestOnlineMap::GetDestinationNDofsPerElement()
00498 {
00499     return m_nDofsPEl_Dest;
00500 }
00501 
00502 // setters for m_nDofsPEl_Src, m_nDofsPEl_Dest
00503 inline void moab::TempestOnlineMap::SetSourceNDofsPerElement( int ns )
00504 {
00505     m_nDofsPEl_Src = ns;
00506 }
00507 
00508 inline void moab::TempestOnlineMap::SetDestinationNDofsPerElement( int nt )
00509 {
00510     m_nDofsPEl_Dest = nt;
00511 }
00512 
00513 ///////////////////////////////////////////////////////////////////////////////
00514 #ifdef MOAB_HAVE_EIGEN3
00515 
00516 inline int moab::TempestOnlineMap::GetSourceGlobalNDofs()
00517 {
00518     return m_weightMatrix.cols();  // return the global number of rows from the weight matrix
00519 }
00520 
00521 // ///////////////////////////////////////////////////////////////////////////////
00522 
00523 inline int moab::TempestOnlineMap::GetDestinationGlobalNDofs()
00524 {
00525     return m_weightMatrix.rows();  // return the global number of columns from the weight matrix
00526 }
00527 
00528 ///////////////////////////////////////////////////////////////////////////////
00529 
00530 inline int moab::TempestOnlineMap::GetSourceLocalNDofs()
00531 {
00532     return m_weightMatrix.cols();  // return the local number of rows from the weight matrix
00533 }
00534 
00535 ///////////////////////////////////////////////////////////////////////////////
00536 
00537 inline int moab::TempestOnlineMap::GetDestinationLocalNDofs()
00538 {
00539     return m_weightMatrix.rows();  // return the local number of columns from the weight matrix
00540 }
00541 
00542 ///////////////////////////////////////////////////////////////////////////////
00543 
00544 inline moab::TempestOnlineMap::WeightMatrix& moab::TempestOnlineMap::GetWeightMatrix()
00545 {
00546     assert( m_weightMatrix.rows() != 0 && m_weightMatrix.cols() != 0 );
00547     return m_weightMatrix;
00548 }
00549 
00550 ///////////////////////////////////////////////////////////////////////////////
00551 
00552 inline moab::TempestOnlineMap::WeightRowVector& moab::TempestOnlineMap::GetRowVector()
00553 {
00554     assert( m_rowVector.size() != 0 );
00555     return m_rowVector;
00556 }
00557 
00558 ///////////////////////////////////////////////////////////////////////////////
00559 
00560 inline moab::TempestOnlineMap::WeightColVector& moab::TempestOnlineMap::GetColVector()
00561 {
00562     assert( m_colVector.size() != 0 );
00563     return m_colVector;
00564 }
00565 
00566 ///////////////////////////////////////////////////////////////////////////////
00567 
00568 #endif
00569 
00570 }  // namespace moab
00571 
00572 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines