MOAB: Mesh Oriented datABase  (version 5.2.1)
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(
00078         std::string strInputType = "fv", std::string strOutputType = "fv", const int nPin = 1, const int nPout = 1,
00079         bool fBubble = false, int fMonotoneTypeID = 0, bool fVolumetric = false, bool fNoConservation = false,
00080         bool fNoCheck = false, const std::string srcDofTagName = "GLOBAL_ID",
00081         const std::string tgtDofTagName = "GLOBAL_ID", const std::string strVariables = "",
00082         const std::string strInputData = "", const std::string strOutputData = "", const std::string strNColName = "",
00083         const bool fOutputDouble = false, const std::string strPreserveVariables = "", const bool fPreserveAll = false,
00084         const double dFillValueOverride = 0.0, const bool fInputConcave = false, const bool fOutputConcave = false );
00085 
00086     /// <summary>
00087     ///     Generate the metadata associated with the offline map.
00088     /// </summary>
00089     // moab::ErrorCode GenerateMetaData();
00090 
00091     /// <summary>
00092     ///     Read the OfflineMap from a NetCDF file.
00093     /// </summary>
00094     moab::ErrorCode ReadParallelMap( const char* strSource, const std::vector< int >& owned_dof_ids,
00095                                      bool row_major_ownership = true );
00096 
00097     /// <summary>
00098     ///     Write the TempestOnlineMap to a parallel NetCDF file.
00099     /// </summary>
00100     moab::ErrorCode WriteParallelMap( const std::string& strTarget );
00101 
00102     /// <summary>
00103     ///     Determine if the map is first-order accurate.
00104     /// </summary>
00105     virtual int IsConsistent( double dTolerance );
00106 
00107     /// <summary>
00108     ///     Determine if the map is conservative.
00109     /// </summary>
00110     virtual int IsConservative( double dTolerance );
00111 
00112     /// <summary>
00113     ///     Determine if the map is monotone.
00114     /// </summary>
00115     virtual int IsMonotone( double dTolerance );
00116 
00117     /// <summary>
00118     ///     If we computed the reduction, get the vector representing the source areas for all entities
00119     /// in the mesh
00120     /// </summary>
00121     const DataArray1D< double >& GetGlobalSourceAreas() const;
00122 
00123     /// <summary>
00124     ///     If we computed the reduction, get the vector representing the target areas for all entities
00125     /// in the mesh
00126     /// </summary>
00127     const DataArray1D< double >& GetGlobalTargetAreas() const;
00128 
00129   private:
00130     /// <summary>
00131     ///     Compute the remapping weights as a permutation matrix that relates DoFs on the source
00132     /// mesh
00133     ///     to DoFs on the target mesh.
00134     /// </summary>
00135     moab::ErrorCode LinearRemapNN_MOAB( bool use_GID_matching = false, bool strict_check = false );
00136 
00137     /// <summary>
00138     ///     Compute the remapping weights for a FV field defined on the source to a
00139     ///     FV field defined on the target mesh.
00140     /// </summary>
00141     void LinearRemapFVtoFV_Tempest_MOAB( int nOrder );
00142 
00143     /// <summary>
00144     ///     Generate the OfflineMap for linear conserative element-average
00145     ///     spectral element to element average remapping.
00146     /// </summary>
00147     void LinearRemapSE0_Tempest_MOAB( const DataArray3D< int >& dataGLLNodes,
00148                                       const DataArray3D< double >& dataGLLJacobian );
00149 
00150     /// <summary>
00151     ///     Generate the OfflineMap for cubic conserative element-average
00152     ///     spectral element to element average remapping.
00153     /// </summary>
00154     void LinearRemapSE4_Tempest_MOAB( const DataArray3D< int >& dataGLLNodes,
00155                                       const DataArray3D< double >& dataGLLJacobian, int nMonotoneType,
00156                                       bool fContinuousIn, 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, const DataArray3D< double >& dataGLLJacobian,
00163                                   const DataArray1D< double >& dataGLLNodalArea, int nOrder, int nMonotoneType,
00164                                   bool fContinuous, bool fNoConservation );
00165 
00166     /// <summary>
00167     ///     Generate the OfflineMap for remapping from finite elements to finite
00168     ///     elements.
00169     /// </summary>
00170     void LinearRemapGLLtoGLL2_MOAB( const DataArray3D< int >& dataGLLNodesIn,
00171                                     const DataArray3D< double >& dataGLLJacobianIn,
00172                                     const DataArray3D< int >& dataGLLNodesOut,
00173                                     const DataArray3D< double >& dataGLLJacobianOut,
00174                                     const DataArray1D< double >& dataNodalAreaOut, int nPin, int nPout,
00175                                     int nMonotoneType, bool fContinuousIn, bool fContinuousOut, bool fNoConservation );
00176 
00177     /// <summary>
00178     ///     Generate the OfflineMap for remapping from finite elements to finite
00179     ///     elements (pointwise interpolation).
00180     /// </summary>
00181     void LinearRemapGLLtoGLL2_Pointwise_MOAB( const DataArray3D< int >& dataGLLNodesIn,
00182                                               const DataArray3D< double >& dataGLLJacobianIn,
00183                                               const DataArray3D< int >& dataGLLNodesOut,
00184                                               const DataArray3D< double >& dataGLLJacobianOut,
00185                                               const DataArray1D< double >& dataNodalAreaOut, int nPin, int nPout,
00186                                               int nMonotoneType, bool fContinuousIn, bool fContinuousOut );
00187 
00188     /// <summary>
00189     ///     Copy the local matrix from Tempest SparseMatrix representation (ELL)
00190     ///     to the parallel CSR Eigen Matrix for scalable application of matvec
00191     ///     needed for projections.
00192     /// </summary>
00193 #ifdef MOAB_HAVE_EIGEN3
00194     void copy_tempest_sparsemat_to_eigen3();
00195 #endif
00196 
00197     /// <summary>
00198     ///     Parallel I/O with HDF5 to write out the remapping weights from multiple processors.
00199     /// </summary>
00200     moab::ErrorCode WriteSCRIPMapFile( const std::string& strOutputFile );
00201 
00202     /// <summary>
00203     ///     Parallel I/O with NetCDF to write out the SCRIP file from multiple processors.
00204     /// </summary>
00205     moab::ErrorCode WriteHDF5MapFile( const std::string& filename );
00206 
00207 
00208   public:
00209     /// <summary>
00210     ///     Store the tag names associated with global DoF ids for source and target meshes
00211     /// </summary>
00212     moab::ErrorCode SetDOFmapTags( const std::string srcDofTagName, const std::string tgtDofTagName );
00213 
00214     /// <summary>
00215     ///     Compute the association between the solution tag global DoF numbering and
00216     ///     the local matrix numbering so that matvec operations can be performed
00217     ///     consistently.
00218     /// </summary>
00219     moab::ErrorCode SetDOFmapAssociation( DiscretizationType srcType, bool isSrcContinuous,
00220                                           DataArray3D< int >* srcdataGLLNodes, DataArray3D< int >* srcdataGLLNodesSrc,
00221                                           DiscretizationType destType, bool isDestContinuous,
00222                                           DataArray3D< int >* tgtdataGLLNodes );
00223 
00224 #ifdef MOAB_HAVE_EIGEN3
00225 
00226     typedef Eigen::Matrix< double, 1, Eigen::Dynamic > WeightDRowVector;
00227     typedef Eigen::Matrix< double, Eigen::Dynamic, 1 > WeightDColVector;
00228     typedef Eigen::SparseVector< double > WeightSVector;
00229     typedef Eigen::SparseMatrix< double, Eigen::RowMajor > WeightRMatrix;
00230     typedef Eigen::SparseMatrix< double, Eigen::ColMajor > WeightCMatrix;
00231 
00232     typedef WeightDRowVector WeightRowVector;
00233     typedef WeightDColVector WeightColVector;
00234     typedef WeightRMatrix WeightMatrix;
00235 
00236     /// <summary>
00237     ///     Get the raw reference to the Eigen weight matrix representing the projection from source to
00238     /// destination mesh.
00239     /// </summary>
00240     WeightMatrix& GetWeightMatrix();
00241 
00242     /// <summary>
00243     ///     Get the row vector that is amenable for application of A*x operation.
00244     /// </summary>
00245     WeightRowVector& GetRowVector();
00246 
00247     /// <summary>
00248     ///     Get the column vector that is amenable for application of A^T*x operation.
00249     /// </summary>
00250     WeightColVector& GetColVector();
00251 
00252 #endif
00253 
00254     /// <summary>
00255     ///     Get the number of total Degrees-Of-Freedom defined on the source mesh.
00256     /// </summary>
00257     int GetSourceGlobalNDofs();
00258 
00259     /// <summary>
00260     ///     Get the number of total Degrees-Of-Freedom defined on the destination mesh.
00261     /// </summary>
00262     int GetDestinationGlobalNDofs();
00263 
00264     /// <summary>
00265     ///     Get the number of local Degrees-Of-Freedom defined on the source mesh.
00266     /// </summary>
00267     int GetSourceLocalNDofs();
00268 
00269     /// <summary>
00270     ///     Get the number of local Degrees-Of-Freedom defined on the destination mesh.
00271     /// </summary>
00272     int GetDestinationLocalNDofs();
00273 
00274     /// <summary>
00275     ///     Get the number of Degrees-Of-Freedom per element on the source mesh.
00276     /// </summary>
00277     int GetSourceNDofsPerElement();
00278 
00279     /// <summary>
00280     ///     Get the number of Degrees-Of-Freedom per element on the destination mesh.
00281     /// </summary>
00282     int GetDestinationNDofsPerElement();
00283 
00284     /// <summary>
00285     ///     Get the global Degrees-Of-Freedom ID on the destination mesh.
00286     /// </summary>
00287     int GetRowGlobalDoF( int localID ) const;
00288 
00289     /// <summary>
00290     ///     Get the global Degrees-Of-Freedom ID on the source mesh.
00291     /// </summary>
00292     int GetColGlobalDoF( int localID ) const;
00293 
00294     /// <summary>
00295     ///     Apply the weight matrix onto the source vector provided as input, and return the column
00296     /// vector (solution projection) after the map application
00297     ///     Compute:        \p tgtVals = A(S->T) * \srcVals, or
00298     ///     if (transpose)  \p tgtVals = [A(T->S)]^T * \srcVals
00299     /// </summary>
00300     moab::ErrorCode ApplyWeights( std::vector< double >& srcVals, std::vector< double >& tgtVals,
00301                                   bool transpose = false );
00302 
00303     /// <summary>
00304     ///     Apply the weight matrix onto the source vector (tag) provided as input, and return the
00305     /// column vector (solution projection) in a tag, after the map application
00306     ///     Compute:        \p tgtVals = A(S->T) * \srcVals, or
00307     ///     if (transpose)  \p tgtVals = [A(T->S)]^T * \srcVals
00308     /// </summary>
00309     moab::ErrorCode ApplyWeights( moab::Tag srcSolutionTag, moab::Tag tgtSolutionTag, bool transpose = false );
00310 
00311     typedef double ( *sample_function )( double, double );
00312 
00313     /// <summary>
00314     ///     Define an analytical solution over the given (source or target) mesh, as specificed in
00315     ///     the context. This routine will define a tag that is compatible with the specified
00316     ///     discretization method type and order and sample the solution exactly using the
00317     ///     analytical function provided by the user.
00318     /// </summary>
00319     moab::ErrorCode DefineAnalyticalSolution( moab::Tag& exactSolnTag, const std::string& solnName,
00320                                               Remapper::IntersectionContext ctx, sample_function testFunction,
00321                                               moab::Tag* clonedSolnTag = NULL, std::string cloneSolnName = "" );
00322 
00323     /// <summary>
00324     ///     Compute the error between a sampled (exact) solution and a projected solution in various
00325     ///     error norms.
00326     /// </summary>
00327     moab::ErrorCode ComputeMetrics( Remapper::IntersectionContext ctx, moab::Tag& exactTag, moab::Tag& approxTag,
00328                                     std::map< std::string, double >& metrics, bool verbose = true );
00329 
00330   private:
00331 
00332 #ifdef MOAB_HAVE_MPI
00333     int rearrange_arrays_by_dofs( const std::vector<unsigned int> & gdofmap,
00334                 DataArray1D< double > &  vecFaceArea,
00335                 DataArray1D< double > &  dCenterLon,
00336                 DataArray1D< double > & dCenterLat,
00337                 DataArray2D< double > & dVertexLat,
00338                 DataArray2D< double > & dVertexLon,
00339                 unsigned & N, // this will be output too now
00340                 int nv,
00341                 int & maxdof);
00342 #endif
00343 
00344     /// <summary>
00345     ///     The fundamental remapping operator object.
00346     /// </summary>
00347     moab::TempestRemapper* m_remapper;
00348 
00349 #ifdef MOAB_HAVE_EIGEN3
00350 
00351     WeightMatrix m_weightMatrix;
00352     WeightRowVector m_rowVector;
00353     WeightColVector m_colVector;
00354 
00355 #endif
00356 
00357     /// <summary>
00358     ///     The reference to the moab::Core object that contains source/target and overlap sets.
00359     /// </summary>
00360     moab::Interface* m_interface;
00361 
00362 #ifdef MOAB_HAVE_MPI
00363     /// <summary>
00364     ///     The reference to the parallel communicator object used by the Core object.
00365     /// </summary>
00366     moab::ParallelComm* m_pcomm;
00367 #endif
00368 
00369     /// <summary>
00370     ///     The original tag data and local to global DoF mapping to associate matrix values to
00371     /// solution    <summary>
00372     moab::Tag m_dofTagSrc, m_dofTagDest;
00373     std::vector< unsigned > row_gdofmap, col_gdofmap, srccol_gdofmap;
00374     std::vector< unsigned > row_dtoc_dofmap, col_dtoc_dofmap, srccol_dtoc_dofmap;
00375 
00376     DataArray3D< int > dataGLLNodesSrc, dataGLLNodesSrcCov, dataGLLNodesDest;
00377     DiscretizationType m_srcDiscType, m_destDiscType;
00378     int m_nTotDofs_Src, m_nTotDofs_SrcCov, m_nTotDofs_Dest;
00379 
00380     // Key details about the current map
00381     int m_nDofsPEl_Src, m_nDofsPEl_Dest;
00382     DiscretizationType m_eInputType, m_eOutputType;
00383     bool m_bConserved;
00384     int m_iMonotonicity;
00385 
00386     Mesh* m_meshInput;
00387     Mesh* m_meshInputCov;
00388     Mesh* m_meshOutput;
00389     Mesh* m_meshOverlap;
00390 
00391     bool is_parallel, is_root;
00392     int rank, size, root_proc;
00393 };
00394 
00395 ///////////////////////////////////////////////////////////////////////////////
00396 
00397 inline int moab::TempestOnlineMap::GetRowGlobalDoF( int localRowID ) const
00398 {
00399     return row_gdofmap[localRowID];
00400 }
00401 
00402 ///////////////////////////////////////////////////////////////////////////////
00403 
00404 inline int moab::TempestOnlineMap::GetColGlobalDoF( int localColID ) const
00405 {
00406     return col_gdofmap[localColID];
00407 }
00408 
00409 ///////////////////////////////////////////////////////////////////////////////
00410 #ifdef MOAB_HAVE_EIGEN3
00411 
00412 inline int moab::TempestOnlineMap::GetSourceGlobalNDofs()
00413 {
00414     return m_weightMatrix.cols();  // return the global number of rows from the weight matrix
00415 }
00416 
00417 // ///////////////////////////////////////////////////////////////////////////////
00418 
00419 inline int moab::TempestOnlineMap::GetDestinationGlobalNDofs()
00420 {
00421     return m_weightMatrix.rows();  // return the global number of columns from the weight matrix
00422 }
00423 
00424 ///////////////////////////////////////////////////////////////////////////////
00425 
00426 inline int moab::TempestOnlineMap::GetSourceLocalNDofs()
00427 {
00428     return m_weightMatrix.cols();  // return the local number of rows from the weight matrix
00429 }
00430 
00431 ///////////////////////////////////////////////////////////////////////////////
00432 
00433 inline int moab::TempestOnlineMap::GetDestinationLocalNDofs()
00434 {
00435     return m_weightMatrix.rows();  // return the local number of columns from the weight matrix
00436 }
00437 
00438 ///////////////////////////////////////////////////////////////////////////////
00439 
00440 inline int moab::TempestOnlineMap::GetSourceNDofsPerElement()
00441 {
00442     return m_nDofsPEl_Src;
00443 }
00444 
00445 ///////////////////////////////////////////////////////////////////////////////
00446 
00447 inline int moab::TempestOnlineMap::GetDestinationNDofsPerElement()
00448 {
00449     return m_nDofsPEl_Dest;
00450 }
00451 
00452 ///////////////////////////////////////////////////////////////////////////////
00453 
00454 inline moab::TempestOnlineMap::WeightMatrix& moab::TempestOnlineMap::GetWeightMatrix()
00455 {
00456     assert( m_weightMatrix.rows() != 0 && m_weightMatrix.cols() != 0 );
00457     return m_weightMatrix;
00458 }
00459 
00460 ///////////////////////////////////////////////////////////////////////////////
00461 
00462 inline moab::TempestOnlineMap::WeightRowVector& moab::TempestOnlineMap::GetRowVector()
00463 {
00464     assert( m_rowVector.size() != 0 );
00465     return m_rowVector;
00466 }
00467 
00468 ///////////////////////////////////////////////////////////////////////////////
00469 
00470 inline moab::TempestOnlineMap::WeightColVector& moab::TempestOnlineMap::GetColVector()
00471 {
00472     assert( m_colVector.size() != 0 );
00473     return m_colVector;
00474 }
00475 
00476 ///////////////////////////////////////////////////////////////////////////////
00477 
00478 #endif
00479 
00480 }  // namespace moab
00481 
00482 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines