Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
TempestOnlineMap.hpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 ///
3 /// \file TempestOnlineMap.hpp
4 /// \author Vijay Mahadevan
5 /// \version November 20, 2017
6 ///
7 
8 #ifndef _TEMPESTONLINEMAP_H_
9 #define _TEMPESTONLINEMAP_H_
10 
11 #include "moab/MOABConfig.h"
12 #ifndef MOAB_HAVE_TEMPESTREMAP
13 #error Re-configure with TempestRemap
14 #endif
15 
17 
18 // Tempest includes
19 #pragma GCC diagnostic push
20 #pragma GCC diagnostic ignored "-Wunused-but-set-variable"
21 #pragma GCC diagnostic ignored "-Wunused-variable"
22 #include "OfflineMap.h"
23 
24 #ifdef MOAB_HAVE_EIGEN3
25 // #ifdef MOAB_HAVE_BLAS
26 // #define EIGEN_USE_BLAS
27 // #endif
28 // #define EIGEN_DONT_PARALLELIZE
29 // #define EIGEN_DONT_VECTORIZE
30 // #define EIGEN_STRONG_INLINE
31 #define EIGEN_DEFAULT_TO_ROW_MAJOR
32 #define EIGEN_RUNTIME_NO_MALLOC
33 // #define EIGEN_UNROLLING_LIMIT 0
34 #include <Eigen/Sparse>
35 #endif
36 
37 #include <unordered_set>
38 
39 #pragma GCC diagnostic pop
40 
41 ///////////////////////////////////////////////////////////////////////////////
42 
43 #if !defined( RECTANGULAR_TRUNCATION ) && !defined( TRIANGULAR_TRUNCATION )
44 #define RECTANGULAR_TRUNCATION
45 // #define TRIANGULAR_TRUNCATION
46 #endif
47 
48 ///////////////////////////////////////////////////////////////////////////////
49 
50 // Forward declarations
51 class Mesh;
52 
53 ///////////////////////////////////////////////////////////////////////////////
54 
55 namespace moab
56 {
57 
58 /// <summary>
59 /// An offline map between two Meshes.
60 /// </summary>
61 // The base class OfflineMap declares private multi-arg overloads alongside public
62 // single-arg virtual functions (IsConsistent, IsConservative, IsMonotone).
63 // Overriding only the public single-arg versions is intentional; suppress the
64 // partial_override warning that NVHPC emits for this pattern.
65 #ifdef __NVCOMPILER
66 #pragma diag_suppress partial_override
67 #endif
68 class TempestOnlineMap : public OfflineMap
69 {
70 
71  public:
72  /// <summary>
73  /// Generate the metadata associated with the offline map.
74  /// </summary>
76 
77  /// <summary>
78  /// Define a virtual destructor.
79  /// </summary>
80  virtual ~TempestOnlineMap();
81 
82  public:
83  // Input / Output types
85  {
90  };
91 
92  // Type of limiter
93  enum CAASType
94  {
95  CAAS_NONE = 0,
99  CAAS_QLT = 4
100  };
101 
102  /// <summary>
103  /// Generate the offline map, given the source and target mesh and discretization details.
104  /// This method generates the mapping between the two meshes based on the overlap and stores
105  /// the result in the SparseMatrix.
106  /// </summary>
107  moab::ErrorCode GenerateRemappingWeights( std::string strInputType,
108  std::string strOutputType,
109  const GenerateOfflineMapAlgorithmOptions& mapOptions,
110  const std::string& srcDofTagName = "GLOBAL_ID",
111  const std::string& tgtDofTagName = "GLOBAL_ID" );
112 
113  /// <summary>
114  /// Generate the metadata associated with the offline map.
115  /// </summary>
116  // moab::ErrorCode GenerateMetaData();
117 
118  /// <summary>
119  /// Read the OfflineMap from a NetCDF file.
120  /// </summary>
121  moab::ErrorCode ReadParallelMap( const char* strSource,
122  const std::vector< int >& tgt_dof_ids,
123  int arearead,
124  std::vector< double >& areaA,
125  int& nA,
126  std::vector< double >& areaB,
127  int& nB );
128 
129  /// <summary>
130  /// Write the TempestOnlineMap to a parallel NetCDF file.
131  /// </summary>
132  moab::ErrorCode WriteParallelMap( const std::string& strTarget,
133  const std::map< std::string, std::string >& attrMap );
134 
135  /// <summary>
136  /// Determine if the map is first-order accurate.
137  /// </summary>
138  virtual int IsConsistent( double dTolerance );
139 
140  /// <summary>
141  /// Determine if the map is conservative.
142  /// </summary>
143  virtual int IsConservative( double dTolerance );
144 
145  /// <summary>
146  /// Determine if the map is monotone.
147  /// </summary>
148  virtual int IsMonotone( double dTolerance );
149 
150  /// <summary>
151  /// If we computed the reduction, get the vector representing the source areas for all entities
152  /// in the mesh
153  /// </summary>
154  const DataArray1D< double >& GetGlobalSourceAreas() const;
155 
156  /// <summary>
157  /// If we computed the reduction, get the vector representing the target areas for all entities
158  /// in the mesh
159  /// </summary>
160  const DataArray1D< double >& GetGlobalTargetAreas() const;
161 
162  /// <summary>
163  /// Print information and metadata about the remapping weights.
164  /// </summary>
165  void PrintMapStatistics();
166 
167  private:
168  /// <summary>
169  /// Compute the remapping weights as a permutation matrix that relates DoFs on the source
170  /// mesh
171  /// to DoFs on the target mesh.
172  /// </summary>
173  moab::ErrorCode LinearRemapNN_MOAB( bool use_GID_matching = false, bool strict_check = false );
174 
175  /// <summary>
176  /// Compute the remapping weights for a FV field defined on the source to a
177  /// FV field defined on the target mesh.
178  /// </summary>
179  void LinearRemapFVtoFV_Tempest_MOAB( int nOrder );
180 
181  /// <summary>
182  /// Generate the OfflineMap for linear conserative element-average
183  /// spectral element to element average remapping.
184  /// </summary>
185  void LinearRemapSE0_Tempest_MOAB( const DataArray3D< int >& dataGLLNodes,
186  const DataArray3D< double >& dataGLLJacobian );
187 
188  /// <summary>
189  /// Generate the OfflineMap for cubic conserative element-average
190  /// spectral element to element average remapping.
191  /// </summary>
192  void LinearRemapSE4_Tempest_MOAB( const DataArray3D< int >& dataGLLNodes,
193  const DataArray3D< double >& dataGLLJacobian,
194  int nMonotoneType,
195  bool fContinuousIn,
196  bool fNoConservation );
197 
198  /// <summary>
199  /// Generate the OfflineMap for remapping from finite volumes to finite
200  /// elements.
201  /// </summary>
202  void LinearRemapFVtoGLL_MOAB( const DataArray3D< int >& dataGLLNodes,
203  const DataArray3D< double >& dataGLLJacobian,
204  const DataArray1D< double >& dataGLLNodalArea,
205  int nOrder,
206  int nMonotoneType,
207  bool fContinuous,
208  bool fNoConservation );
209 
210  /// <summary>
211  /// Generate the OfflineMap for remapping from finite elements to finite
212  /// elements.
213  /// </summary>
214  void LinearRemapGLLtoGLL2_MOAB( const DataArray3D< int >& dataGLLNodesIn,
215  const DataArray3D< double >& dataGLLJacobianIn,
216  const DataArray3D< int >& dataGLLNodesOut,
217  const DataArray3D< double >& dataGLLJacobianOut,
218  const DataArray1D< double >& dataNodalAreaOut,
219  int nPin,
220  int nPout,
221  int nMonotoneType,
222  bool fContinuousIn,
223  bool fContinuousOut,
224  bool fNoConservation );
225 
226  /// <summary>
227  /// Generate the OfflineMap for remapping from finite elements to finite
228  /// elements (pointwise interpolation).
229  /// </summary>
230  void LinearRemapGLLtoGLL2_Pointwise_MOAB( const DataArray3D< int >& dataGLLNodesIn,
231  const DataArray3D< double >& dataGLLJacobianIn,
232  const DataArray3D< int >& dataGLLNodesOut,
233  const DataArray3D< double >& dataGLLJacobianOut,
234  const DataArray1D< double >& dataNodalAreaOut,
235  int nPin,
236  int nPout,
237  int nMonotoneType,
238  bool fContinuousIn,
239  bool fContinuousOut );
240 
241  /// <summary>
242  /// Copy the local matrix from Tempest SparseMatrix representation (ELL)
243  /// to the parallel CSR Eigen Matrix for scalable application of matvec
244  /// needed for projections.
245  /// </summary>
246 #ifdef MOAB_HAVE_EIGEN3
247  void copy_tempest_sparsemat_to_eigen3();
248 #endif
249 
250  /// <summary>
251  /// Parallel I/O with HDF5 to write out the remapping weights from multiple processors.
252  /// </summary>
253  moab::ErrorCode WriteSCRIPMapFile( const std::string& strOutputFile,
254  const std::map< std::string, std::string >& attrMap );
255 
256  /// <summary>
257  /// Parallel I/O with NetCDF to write out the SCRIP file from multiple processors.
258  /// </summary>
259  moab::ErrorCode WriteHDF5MapFile( const std::string& filename );
260 
261  public:
262  /// <summary>
263  /// Store the tag names associated with global DoF ids for source and target meshes to be used for mapping.
264  /// <param name="srcDofTagName">The tag name associated with global DoF ids for the source mesh</param>
265  /// <param name="tgtDofTagName">The tag name associated with global DoF ids for the target mesh</param>
266  /// </summary>
267  moab::ErrorCode SetDOFmapTags( const std::string srcDofTagName, const std::string tgtDofTagName );
268 
269  /// <summary>
270  /// @brief Compute the association between the solution tag global DoF numbering and
271  /// the local matrix numbering so that matvec operations can be performed
272  /// consistently.
273  /// <param name="srcType">The discretization type of the source mesh</param>
274  /// <param name="srcOrder">The order of the discretization on the source mesh</param>
275  /// <param name="isSrcContinuous">The continuity of the discretization on the source mesh</param>
276  /// <param name="srcdataGLLNodes">The GLL nodes on the source mesh</param>
277  /// <param name="srcdataGLLNodesSrc">The GLL nodes on the source mesh</param>
278  /// <param name="destType">The discretization type of the destination mesh</param>
279  /// <param name="destOrder">The order of the discretization on the destination mesh</param>
280  /// <param name="isTgtContinuous">The continuity of the discretization on the destination mesh</param>
281  /// <param name="tgtdataGLLNodes">The GLL nodes on the destination mesh</param>
282  /// </summary>
284  int srcOrder,
285  bool isSrcContinuous,
286  DataArray3D< int >* srcdataGLLNodes,
287  DataArray3D< int >* srcdataGLLNodesSrc,
288  DiscretizationType destType,
289  int destOrder,
290  bool isTgtContinuous,
291  DataArray3D< int >* tgtdataGLLNodes );
292 
293  /// <summary>
294  /// @brief ApplyBoundsLimiting - Apply bounds limiting to the data field
295  /// @param dataInDouble - input data field
296  /// @param dataOutDouble - output data field
297  /// @param caasType - type of limiter
298  /// @param caasIteration - iteration number of limiter
299  /// @return - pair of mass defect pre and post limiter application
300  /// </summary>
301  std::pair< double, double > ApplyBoundsLimiting( std::vector< double >& dataInDouble,
302  std::vector< double >& dataOutDouble,
303  CAASType caasType = CAAS_GLOBAL,
304  int caasIteration = 0,
305  double mismatch = 0.0 );
306 
307  /// @brief
308  /// @param vecAdjFaces
309  /// @param nrings
310  /// @param entities
311  /// @param useMOABAdjacencies
312  /// @param trMesh
313  /// @return
314  void ComputeAdjacencyRelations( std::vector< std::unordered_set< int > >& vecAdjFaces,
315  int nrings,
316  const Range& entities,
317  bool useMOABAdjacencies = true,
318  Mesh* trMesh = nullptr );
319 
320 #ifdef MOAB_HAVE_EIGEN3
321 
322  typedef Eigen::Matrix< double, 1, Eigen::Dynamic > WeightDRowVector;
323  typedef Eigen::Matrix< double, Eigen::Dynamic, 1 > WeightDColVector;
324  typedef Eigen::SparseVector< double > WeightSVector;
325  typedef Eigen::SparseMatrix< double, Eigen::RowMajor > WeightRMatrix;
326  typedef Eigen::SparseMatrix< double, Eigen::ColMajor > WeightCMatrix;
327 
328  typedef WeightDRowVector WeightRowVector;
329  typedef WeightDColVector WeightColVector;
330  typedef WeightRMatrix WeightMatrix;
331 
332  /// <summary>
333  /// Get the raw reference to the Eigen weight matrix representing the projection from source to
334  /// destination mesh.
335  /// </summary>
336  WeightMatrix& GetWeightMatrix();
337 
338  /// <summary>
339  /// Get the row vector that is amenable for application of A*x operation.
340  /// </summary>
341  WeightRowVector& GetRowVector();
342 
343  /// <summary>
344  /// Get the column vector that is amenable for application of A^T*x operation.
345  /// </summary>
346  WeightColVector& GetColVector();
347 
348 #endif
349 
350  /// <summary>
351  /// Get the number of total Degrees-Of-Freedom defined on the source mesh.
352  /// </summary>
354 
355  /// <summary>
356  /// Get the number of total Degrees-Of-Freedom defined on the destination mesh.
357  /// </summary>
359 
360  /// <summary>
361  /// Get the number of local Degrees-Of-Freedom defined on the source mesh.
362  /// </summary>
364 
365  /// <summary>
366  /// Get the number of local Degrees-Of-Freedom defined on the destination mesh.
367  /// </summary>
369 
370  /// <summary>
371  /// Get the number of Degrees-Of-Freedom per element on the source mesh.
372  /// </summary>
374 
375  /// <summary>
376  /// Get the number of Degrees-Of-Freedom per element on the destination mesh.
377  /// </summary>
379 
380  /// <summary>
381  /// Set the number of Degrees-Of-Freedom per element on the source mesh.
382  /// </summary>
383  void SetSourceNDofsPerElement( int ns );
384 
385  /// <summary>
386  /// Get the number of Degrees-Of-Freedom per element on the destination mesh.
387  /// </summary>
388  void SetDestinationNDofsPerElement( int nt );
389 
390  /// <summary>
391  /// Get the global Degrees-Of-Freedom ID on the destination mesh.
392  /// </summary>
393  int GetRowGlobalDoF( int localID ) const;
394 
395  /// <summary>
396  /// Get the index of globaRowDoF.
397  /// </summary>
398  inline int GetIndexOfRowGlobalDoF( int globalRowDoF ) const;
399 
400  /// <summary>
401  /// Get the global Degrees-Of-Freedom ID on the source mesh.
402  /// </summary>
403  int GetColGlobalDoF( int localID ) const;
404 
405  /// <summary>
406  /// Get the index of globaColDoF.
407  /// </summary>
408  inline int GetIndexOfColGlobalDoF( int globalColDoF ) const;
409 
410  /// <summary>
411  /// Apply the weight matrix onto the source vector (tag) provided as input, and return the
412  /// column vector (solution projection) in a tag, after the map application
413  /// Compute: \p tgtVals = A(S->T) * \note Source values, or
414  /// if (transpose) \p tgtVals = [A(T->S)]^T * \note Source values
415  /// </summary>
416  moab::ErrorCode ApplyWeights( moab::Tag srcSolutionTag,
417  moab::Tag tgtSolutionTag,
418  bool transpose = false,
419  CAASType caasType = CAAS_NONE,
420  double default_projection = 0.0 );
421 
422  /// <summary>
423  /// Apply the high-order weight matrix onto the source vector (tag), and then enforce
424  /// bounds computed from the stencil of a separate low-order weight map using the
425  /// Clip-And-Assert-Sum (CAAS) algorithm. This implements the dual-map nonlinear
426  /// remapping pattern used in E3SM coupling.
427  /// \p loWeightMap provides the low-order (monotone) map whose per-row stencil defines
428  /// the min/max bounds. The CAAS filter clips the high-order result to those bounds
429  /// and redistributes mass proportionally to maintain conservation.
430  /// </summary>
432  moab::Tag tgtSolutionTag,
433  TempestOnlineMap* loWeightMap,
434  CAASType caasType = CAAS_LOCAL );
435 
436  typedef double ( *sample_function )( double, double );
437 
438  /// <summary>
439  /// Define an analytical solution over the given (source or target) mesh, as specificed in
440  /// the context. This routine will define a tag that is compatible with the specified
441  /// discretization method type and order and sample the solution exactly using the
442  /// analytical function provided by the user.
443  /// </summary>
445  const std::string& solnName,
447  sample_function testFunction,
448  moab::Tag* clonedSolnTag = NULL,
449  std::string cloneSolnName = "" );
450 
451  /// <summary>
452  /// Compute the error between a sampled (exact) solution and a projected solution in various
453  /// error norms.
454  /// </summary>
456  moab::Tag& exactTag,
457  moab::Tag& approxTag,
458  std::map< std::string, double >& metrics,
459  bool verbose = true );
460 
461  moab::ErrorCode fill_col_ids( std::vector< int >& ids_of_interest )
462  {
463  ids_of_interest.reserve( col_gdofmap.size() );
464  // need to add 1
465  for( auto it = col_gdofmap.begin(); it != col_gdofmap.end(); it++ )
466  ids_of_interest.push_back( *it + 1 );
467  return moab::MB_SUCCESS;
468  }
469 
470  moab::ErrorCode set_col_dc_dofs( std::vector< int >& values_entities );
471 
472  moab::ErrorCode set_row_dc_dofs( std::vector< int >& values_entities );
473 
474  // hack
475  void SetMeshInput( Mesh* imesh )
476  {
477  m_meshInput = imesh;
478  };
479 
480  /// Read-only access to the matrix-row -> matrix-col DOF index maps. Used
481  /// by callers (e.g. iMOAB diagnostic helpers) that need to translate
482  /// matrix indices back into source/target tag-vector indices.
483  const std::vector< int >& GetRowDofMap() const { return row_dtoc_dofmap; }
484  const std::vector< int >& GetColDofMap() const { return col_dtoc_dofmap; }
485 
486  private:
487  template < typename SparseMatrixType >
488  void serializeSparseMatrix( const SparseMatrixType& mat, const std::string& filename );
489 
490  void setup_sizes_dimensions();
491 
492  void CAASLimiter( std::vector< double >& dataCorrectedField,
493  std::vector< double >& dataLowerBound,
494  std::vector< double >& dataUpperBound,
495  double& dMass );
496  double QLTLimiter( int caasIteration,
497  std::vector< double >& dataCorrectedField,
498  std::vector< double >& dataLowerBound,
499  std::vector< double >& dataUpperBound,
500  std::vector< double >& dMassDefect );
501 
502  /// <summary>
503  /// Apply the weight matrix onto the source vector provided as input, and return the column
504  /// vector (solution projection) after the map application
505  /// Compute: \p tgtVals = A(S->T) * \note Source values, or
506  /// if (transpose) \p tgtVals = [A(T->S)]^T * \note Source values
507  /// </summary>
508  moab::ErrorCode ApplyWeights( std::vector< double >& srcVals,
509  std::vector< double >& tgtVals,
510  bool transpose = false );
511 
512 #ifdef MOAB_HAVE_MPI
513  int rearrange_arrays_by_dofs( const std::vector< unsigned int >& gdofmap,
514  DataArray1D< double >& vecFaceArea,
515  DataArray1D< double >& dCenterLon,
516  DataArray1D< double >& dCenterLat,
517  DataArray2D< double >& dVertexLat,
518  DataArray2D< double >& dVertexLon,
519  std::vector< int >& masks,
520  unsigned& N, // this will be output too now
521  int nv,
522  int& maxdof );
523 #endif
524 
525  /// <summary>
526  /// The fundamental remapping operator object.
527  /// </summary>
529 
530 #ifdef MOAB_HAVE_EIGEN3
531 
532  int num_rows, num_cols;
533  WeightMatrix m_weightMatrix;
534  WeightRowVector m_rowVector;
535  WeightColVector m_colVector;
536 
537 #endif
538 
539  /// <summary>
540  /// The reference to the moab::Core object that contains source/target and overlap sets.
541  /// </summary>
543 
544 #ifdef MOAB_HAVE_MPI
545  /// <summary>
546  /// The reference to the parallel communicator object used by the Core object.
547  /// </summary>
548  moab::ParallelComm* m_pcomm;
549 #endif
550 
551  /// <summary>
552  /// The original tag data and local to global DoF mapping to associate matrix values to
553  /// solution <summary>
555  std::vector< unsigned > row_gdofmap, col_gdofmap, srccol_gdofmap;
556 
557  // make it int, because it can be -1 in new logic
559 
560  std::map< int, int > rowMap, colMap;
562 
566 
567  // Key details about the current map
572 
573  Mesh* m_meshInput;
577 
579  int rank, size;
580 };
581 #ifdef __NVCOMPILER
582 #pragma diag_default partial_override
583 #endif
584 
585 ///////////////////////////////////////////////////////////////////////////////
586 
587 inline int moab::TempestOnlineMap::GetRowGlobalDoF( int localRowID ) const
588 {
589  return row_gdofmap[localRowID];
590 }
591 
592 inline int moab::TempestOnlineMap::GetIndexOfRowGlobalDoF( int globalRowDoF ) const /* 0 based */
593 {
594  return globalRowDoF + 1;
595 }
596 ///////////////////////////////////////////////////////////////////////////////
597 
598 inline int moab::TempestOnlineMap::GetColGlobalDoF( int localColID ) const
599 {
600  return col_gdofmap[localColID];
601 }
602 
603 inline int moab::TempestOnlineMap::GetIndexOfColGlobalDoF( int globalColDoF ) const /* 0 based */
604 {
605  return globalColDoF + 1; // temporary
606 }
607 ///////////////////////////////////////////////////////////////////////////////
608 
610 {
611  return m_nDofsPEl_Src;
612 }
613 
614 ///////////////////////////////////////////////////////////////////////////////
615 
617 {
618  return m_nDofsPEl_Dest;
619 }
620 
621 // setters for m_nDofsPEl_Src, m_nDofsPEl_Dest
623 {
624  m_nDofsPEl_Src = ns;
625 }
626 
628 {
629  m_nDofsPEl_Dest = nt;
630 }
631 
632 ///////////////////////////////////////////////////////////////////////////////
633 #ifdef MOAB_HAVE_EIGEN3
634 
636 {
637  return m_weightMatrix.cols(); // return the global number of rows from the weight matrix
638 }
639 
640 // ///////////////////////////////////////////////////////////////////////////////
641 
643 {
644  return m_weightMatrix.rows(); // return the global number of columns from the weight matrix
645 }
646 
647 ///////////////////////////////////////////////////////////////////////////////
648 
650 {
651  return m_weightMatrix.cols(); // return the local number of rows from the weight matrix
652 }
653 
654 ///////////////////////////////////////////////////////////////////////////////
655 
657 {
658  return m_weightMatrix.rows(); // return the local number of columns from the weight matrix
659 }
660 
661 ///////////////////////////////////////////////////////////////////////////////
662 
663 inline moab::TempestOnlineMap::WeightMatrix& moab::TempestOnlineMap::GetWeightMatrix()
664 {
665  assert( m_weightMatrix.rows() != 0 && m_weightMatrix.cols() != 0 );
666  return m_weightMatrix;
667 }
668 
669 ///////////////////////////////////////////////////////////////////////////////
670 
671 inline moab::TempestOnlineMap::WeightRowVector& moab::TempestOnlineMap::GetRowVector()
672 {
673  assert( m_rowVector.size() != 0 );
674  return m_rowVector;
675 }
676 
677 ///////////////////////////////////////////////////////////////////////////////
678 
679 inline moab::TempestOnlineMap::WeightColVector& moab::TempestOnlineMap::GetColVector()
680 {
681  assert( m_colVector.size() != 0 );
682  return m_colVector;
683 }
684 
685 ///////////////////////////////////////////////////////////////////////////////
686 
687 #endif
688 
689 } // namespace moab
690 
691 #endif