Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
TempestRemapper.hpp
Go to the documentation of this file.
1 /*
2  * =====================================================================================
3  *
4  * Filename: TempestRemapper.hpp
5  *
6  * Description: Interface to the TempestRemap library to enable intersection and
7  * high-order conservative remapping of climate solution from
8  * arbitrary resolution of source and target grids on the sphere.
9  *
10  * Author: Vijay S. Mahadevan (vijaysm), [email protected]
11  *
12  * =====================================================================================
13  */
14 
15 #ifndef MB_TEMPESTREMAPPER_HPP
16 #define MB_TEMPESTREMAPPER_HPP
17 
21 
22 // Tempest includes
23 #ifdef MOAB_HAVE_TEMPESTREMAP
24 #include "netcdfcpp.h"
25 #ifdef MOAB_HAVE_MPI
26 #define TEMPEST_MPIOMP
27 #endif
28 #include "Announce.h"
29 #include "TempestRemapAPI.h"
30 #else
31 #error "This tool depends on TempestRemap library. Reconfigure using --with-tempestremap"
32 #endif
33 
34 namespace moab
35 {
36 
37 // Forward declare our friend, the mapper
38 class TempestOnlineMap;
39 
40 class TempestRemapper : public Remapper
41 {
42  public:
43 #ifdef MOAB_HAVE_MPI
44  TempestRemapper( moab::Interface* mbInt, moab::ParallelComm* pcomm = NULL, bool offlineMode = false )
45  : Remapper( mbInt, pcomm ),
46 #else
47  TempestRemapper( moab::Interface* mbInt, bool offlineMode = false )
48  : Remapper( mbInt ),
49 #endif
50  offlineWorkflow( offlineMode ), meshValidate( false ), constructEdgeMap( false ), m_source_type( DEFAULT ),
52  {
53 #ifdef MOAB_HAVE_MPI
54  AnnounceOnlyOutputOnRankZero();
55 #endif
56  }
57 
58  virtual ~TempestRemapper();
59 
60  // Mesh type with a correspondence to Tempest/Climate formats
62  {
63  DEFAULT = -1,
64  CS = 0,
65  RLL = 1,
66  ICO = 2,
67  ICOD = 3,
70  OVERLAP_MOAB = 6
71  };
72 
73  friend class TempestOnlineMap;
74 
75  public:
76  /**
77  * @brief Initialize the TempestRemapper object internal data structures including the mesh sets
78  * and TempestRemap mesh references.
79  *
80  * @param initialize_fsets Flag to initialize the mesh sets (default: true)
81  * @return ErrorCode indicating the status of the initialization
82  */
83  virtual ErrorCode initialize( bool initialize_fsets = true );
84 
85  /**
86  * @brief Deallocate and clear any memory initialized in the TempestRemapper object
87  *
88  * @return ErrorCode indicating the status of the clear operation
89  */
90  virtual ErrorCode clear();
91 
92  /**
93  * @brief Generate a mesh in memory of given type (CS/RLL/ICO/MPAS(structured)) and store it
94  * under the context specified by the user.
95  *
96  * @param ctx Intersection context
97  * @param type Type of mesh to generate
98  * @return ErrorCode indicating the status of the mesh generation
99  */
101 
102  /**
103  * @brief Load a mesh from disk of given type and store it under the context specified by the user.
104  *
105  * @param ctx Intersection context
106  * @param inputFilename File name of the mesh to load
107  * @param type Type of mesh to load
108  * @return ErrorCode indicating the status of the mesh loading
109  */
110  moab::ErrorCode LoadMesh( Remapper::IntersectionContext ctx, std::string inputFilename, TempestMeshType type );
111 
112  /**
113  * @brief Construct a source covering mesh such that it completely encompasses the target grid in
114  * parallel. This operation is critical to ensure that the parallel advancing-front
115  * intersection algorithm can compute the intersection mesh only locally without any process
116  * communication.
117  *
118  * @param tolerance Tolerance for the covering mesh construction (default: 1e-8)
119  * @param radius_src Radius of the source mesh (default: 1.0)
120  * @param radius_tgt Radius of the target mesh (default: 1.0)
121  * @param boxeps Box epsilon value (default: 0.1)
122  * @param regional_mesh Flag to indicate if the mesh is regional (default: false)
123  * @param gnomonic Flag to indicate if the mesh is gnomonic (default: true)
124  * @param nb_ghost_layers Number of ghost layers (default: 0)
125  * @return ErrorCode indicating the status of the covering mesh construction
126  */
128  double radius_src = 1.0,
129  double radius_tgt = 1.0,
130  double boxeps = 0.1,
131  bool regional_mesh = false,
132  bool gnomonic = true,
133  int nb_ghost_layers = 0 );
134 
135  /**
136  * @brief Compute the intersection mesh between the source and target grids that have been
137  * instantiated in the Remapper. This function invokes the parallel advancing-front
138  * intersection algorithm internally for spherical meshes and can handle arbitrary
139  * unstructured grids (CS, RLL, ICO, MPAS) with and without holes.
140  *
141  * @param kdtree_search Flag to enable k-d tree search (default: true)
142  * @param use_tempest Flag to use TempestRemap (default: false)
143  * @return ErrorCode indicating the status of the intersection mesh computation
144  */
145  moab::ErrorCode ComputeOverlapMesh( bool kdtree_search = true, bool use_tempest = false );
146 
147  /**
148  * @brief Convert the TempestRemap mesh object to a corresponding MOAB mesh representation
149  * according to the intersection context.
150  *
151  * @param ctx Intersection context
152  * @return ErrorCode indicating the status of the mesh conversion
153  */
155 
156  /**
157  * @brief Convert the MOAB mesh representation to a corresponding TempestRemap mesh object
158  * according to the intersection context.
159  *
160  * @param ctx Intersection context
161  * @return ErrorCode indicating the status of the mesh conversion
162  */
164 
165  /**
166  * @brief Get the TempestRemap mesh object according to the intersection context.
167  *
168  * @param ctx Intersection context
169  * @return Pointer to the TempestRemap mesh object
170  */
172 
173  /**
174  * @brief Set the TempestRemap mesh object according to the intersection context.
175  *
176  * @param ctx Intersection context
177  * @param mesh Pointer to the TempestRemap mesh object
178  * @param overwrite Flag to overwrite the existing mesh (default: true)
179  */
180  void SetMesh( Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite = true );
181 
182  /**
183  * @brief Set the mesh set according to the intersection context.
184  *
185  * @param ctx Intersection context
186  * @param mset MOAB mesh set handle
187  * @param entities MOAB range of entities (optional)
188  */
189  void SetMeshSet( Remapper::IntersectionContext ctx, moab::EntityHandle mset, moab::Range* entities = nullptr );
190 
191  /**
192  * @brief Get the covering mesh (TempestRemap) object.
193  *
194  * @return Pointer to the covering mesh object
195  */
196  Mesh* GetCoveringMesh();
197 
198  /**
199  * @brief Get the MOAB mesh set corresponding to the intersection context.
200  *
201  * @param ctx Intersection context
202  * @return MOAB mesh set handle
203  */
205 
206  /**
207  * @brief Const overload. Get the MOAB mesh set corresponding to the intersection context.
208  *
209  * @param ctx Intersection context
210  * @return MOAB mesh set handle
211  */
213 
214  /**
215  * @brief Get the mesh element entities corresponding to the intersection context.
216  *
217  * @param ctx Intersection context
218  * @return MOAB range of mesh element entities
219  */
221 
222  /**
223  * @brief Const overload. Get the mesh element entities corresponding to the intersection context.
224  *
225  * @param ctx Intersection context
226  * @return MOAB range of mesh element entities
227  */
229 
230  /**
231  * @brief Get the mesh vertices corresponding to the intersection context. Useful for point-cloud
232  * meshes.
233  *
234  * @param ctx Intersection context
235  * @return MOAB range of mesh vertices
236  */
238 
239  /**
240  * @brief Const overload. Get the mesh vertices corresponding to the intersection context. Useful
241  * for point-cloud meshes.
242  *
243  * @param ctx Intersection context
244  * @return MOAB range of mesh vertices
245  */
247 
248  /**
249  * @brief Get access to the underlying source covering set if available. Else return the source
250  * set.
251  *
252  * @return MOAB entity handle of the covering set
253  */
255 
256  /**
257  * @brief Set the mesh type corresponding to the intersection context
258  *
259  * @param ctx Intersection context
260  * @param metadata Vector of mesh type metadata
261  */
262  void SetMeshType( Remapper::IntersectionContext ctx, const std::vector< int >& metadata );
263 
264  /**
265  * @brief Reconstruct mesh, used now only for IO; need a better solution maybe
266  *
267  * @param ctx Intersection context
268  * @param meshSet MOAB mesh set handle
269  */
271 
272  /**
273  * @brief Get the mesh type corresponding to the intersection context
274  *
275  * @param ctx Intersection context
276  * @return Mesh type
277  */
279 
280  /**
281  * @brief Gather the overlap mesh and associated source/target data and write it out to disk
282  * using the TempestRemap output interface. This information can then be used with the
283  * "GenerateOfflineMap" tool in TempestRemap as needed.
284  *
285  * @param strOutputFileName Output file name
286  * @param fAllParallel Flag to write all parallel data (default: false)
287  * @param fInputConcave Flag to indicate if the input mesh is concave (default: false)
288  * @param fOutputConcave Flag to indicate if the output mesh is concave (default: false)
289  * @return ErrorCode indicating the status of the write operation
290  */
291  moab::ErrorCode WriteTempestIntersectionMesh( std::string strOutputFileName,
292  const bool fAllParallel,
293  const bool fInputConcave,
294  const bool fOutputConcave );
295 
296  /**
297  * @brief Generate the necessary metadata and specifically the GLL node numbering for DoFs for a
298  * CS mesh. This negates the need for running external code like HOMME to output the
299  * numbering needed for computing maps. The functionality is used through the `mbconvert`
300  * tool to compute processor-invariant Global DoF IDs at GLL nodes.
301  *
302  * @param ntot_elements Total number of elements
303  * @param entities MOAB range of entities
304  * @param secondary_entities MOAB range of secondary entities (optional)
305  * @param dofTagName Name of the DoF tag
306  * @param nP Number of points
307  * @return ErrorCode indicating the status of the metadata generation
308  */
309  moab::ErrorCode GenerateCSMeshMetadata( const int ntot_elements,
310  moab::Range& entities,
311  moab::Range* secondary_entities,
312  const std::string& dofTagName,
313  int nP );
314 
315  /**
316  * @brief Generate the necessary metadata for DoF node numbering in a given mesh.
317  * Currently, only the functionality to generate numbering on CS grids is supported.
318  *
319  * @param mesh Mesh object
320  * @param ntot_elements Total number of elements
321  * @param entities MOAB range of entities
322  * @param secondary_entities MOAB range of secondary entities (optional)
323  * @param dofTagName Name of the DoF tag
324  * @param nP Number of points
325  * @return ErrorCode indicating the status of the metadata generation
326  */
328  const int ntot_elements,
329  moab::Range& entities,
330  moab::Range* secondary_entities,
331  const std::string& dofTagName,
332  int nP );
333 
334  /**
335  * @brief Get all the ghosted overlap entities that were accumulated to enable conservation in
336  * parallel
337  *
338  * @param sharedGhostEntities MOAB range of ghosted overlap entities
339  * @return ErrorCode indicating the status of the get operation
340  */
342 
343 #ifndef MOAB_HAVE_MPI
344  /**
345  * @brief Internal method to assign vertex and element global IDs if one does not exist already
346  *
347  * @param idtag Tag for the global IDs
348  * @param this_set MOAB entity handle of the mesh set
349  * @param dimension Dimension of the mesh (default: 2)
350  * @param start_id Starting ID (default: 1)
351  * @return ErrorCode indicating the status of the assignment
352  */
354  EntityHandle this_set,
355  const int dimension = 2,
356  const int start_id = 1 );
357 #endif
358 
359  /**
360  * @brief Get the masks that could have been defined
361  *
362  * @param ctx Intersection context
363  * @param masks Vector of masks
364  * @return ErrorCode indicating the status of the get operation
365  */
366  ErrorCode GetIMasks( Remapper::IntersectionContext ctx, std::vector< int >& masks );
367 
368  public: // public members
369  /**
370  * @brief Flag indicating whether the workflow is in offline mode.
371  *
372  * This flag is used to determine the context of the workflow, specifically
373  * whether it is running in an offline mode (mbtempest).
374  */
375  const bool offlineWorkflow;
376 
377  /**
378  * @brief Flag to enable mesh validation after loading from file.
379  *
380  * If set to true, the mesh will be validated after it is loaded from a file.
381  */
383 
384  /**
385  * @brief Flag to construct the edge map within the TempestRemap data structures.
386  *
387  * If set to true, the edge map will be constructed within the TempestRemap
388  * data structures.
389  */
391 
392  /**
393  * @brief Global verbosity flag.
394  *
395  * This flag controls the verbosity of the output. If set to true, more
396  * detailed output will be generated.
397  */
398  static const bool verbose = true;
399 
400  private:
401  /**
402  * @brief Convert all MOAB meshes to TempestRemap format.
403  *
404  * Utility method primarily used in online workflows to convert all meshes
405  * to TempestRemap format.
406  *
407  * @return ErrorCode Status of the conversion
408  */
410 
411  /**
412  * @brief Convert overlap mesh to source-ordered format.
413  *
414  * Transforms the overlap mesh into a source-ordered format compatible
415  * with TempestRemap.
416  *
417  * @return moab::ErrorCode Status of the conversion
418  */
420 
421  // Private methods
422 
423  /**
424  * @brief Load a mesh from disk in TempestRemap format.
425  *
426  * @param inputFilename Path to the input mesh file
427  * @param tempest_mesh Pointer to store the loaded TempestRemap mesh
428  * @return moab::ErrorCode Status of the load operation
429  */
430  moab::ErrorCode load_tempest_mesh_private( std::string inputFilename, Mesh** tempest_mesh );
431 
432  /**
433  * @brief Convert a MOAB mesh to TempestRemap format.
434  *
435  * @param mesh Target TempestRemap mesh pointer
436  * @param meshset MOAB mesh set to convert
437  * @param entities Range of entities to convert
438  * @param pverts Optional pointer to vertex range
439  * @return moab::ErrorCode Status of the conversion
440  */
442  moab::EntityHandle meshset,
443  moab::Range& entities,
444  moab::Range* pverts );
445 
446  /**
447  * @brief Convert a TempestRemap mesh to MOAB format.
448  *
449  * @param type Type of the TempestRemap mesh
450  * @param mesh Source TempestRemap mesh
451  * @param meshset Target MOAB mesh set
452  * @param entities Range to store converted entities
453  * @param vertices Optional range to store vertices
454  * @return moab::ErrorCode Status of the conversion
455  */
457  Mesh* mesh,
458  moab::EntityHandle& meshset,
459  moab::Range& entities,
460  moab::Range* vertices );
461 
462  /**
463  * @brief Augment overlap mesh with ghosted entities.
464  *
465  * Adds ghosted entities to the overlap mesh to ensure conservation
466  * in parallel operations.
467  *
468  * @return moab::ErrorCode Status of the augmentation
469  */
471 
472  /* Source meshset, mesh and entity references */
473  Mesh* m_source = nullptr;
474  TempestMeshType m_source_type; // initialized in ctor member init list
479  bool point_cloud_source = false;
480  std::vector< int > m_source_metadata;
481 
482  /* Target meshset, mesh and entity references */
483  Mesh* m_target = nullptr;
484  TempestMeshType m_target_type; // initialized in ctor member init list
489  bool point_cloud_target = false;
490  std::vector< int > m_target_metadata;
491 
492  /* Overlap meshset, mesh and entity references */
493  Mesh* m_overlap = nullptr;
497  std::vector< std::pair< int, int > > m_sorted_overlap_order;
498 
499  /* Intersection context on a sphere */
501 
502  /* Parallel - migrated mesh that is in the local view */
503  Mesh* m_covering_source = nullptr;
507 
508  /* local to glboal and global to local ID maps */
509  // std::map< int, int > gid_to_lid_src, gid_to_lid_covsrc, gid_to_lid_tgt;
510  // std::map< int, int > lid_to_gid_src, lid_to_gid_covsrc, lid_to_gid_tgt;
511 
513 
514  bool rrmgrids = false;
515  bool is_parallel = false;
516  bool is_root = false;
517  int rank = 0;
518  int size = 0;
519 };
520 
521 // Inline functions
523 {
524  switch( ctx )
525  {
527  return m_source;
529  return m_target;
531  return m_overlap;
533  return m_covering_source;
534  case Remapper::DEFAULT:
535  default:
536  return NULL;
537  }
538 }
539 
540 inline void TempestRemapper::SetMesh( Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite )
541 {
542  switch( ctx )
543  {
545  if( !overwrite && m_source ) return;
546  if( overwrite && m_source ) delete m_source;
547  m_source = mesh;
548  break;
550  if( !overwrite && m_target ) return;
551  if( overwrite && m_target ) delete m_target;
552  m_target = mesh;
553  break;
555  if( !overwrite && m_overlap ) return;
556  if( overwrite && m_overlap ) delete m_overlap;
557  m_overlap = mesh;
558  break;
560  if( !overwrite && m_covering_source ) return;
561  if( overwrite && m_covering_source ) delete m_covering_source;
562  m_covering_source = mesh;
563  break;
564  case Remapper::DEFAULT:
565  default:
566  break;
567  }
568 }
569 
571 {
572  switch( ctx )
573  {
575  delete m_source;
576  m_source = new Mesh;
577  m_source_set = meshSet;
579  m_source->CalculateFaceAreas( false ); // fInputConcave is false ?
580  break;
582  // not needed yet
583  break;
585  // not needed yet
586  break;
588  // not needed yet
589  break;
590  case Remapper::DEFAULT:
591  default:
592  break;
593  }
594 }
595 
597 {
598  switch( ctx )
599  {
601  return m_source_set;
603  return m_target_set;
605  return m_overlap_set;
607  return m_covering_source_set;
608  case Remapper::DEFAULT:
609  default:
610  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_set );
611  }
612 }
613 
615 {
616  switch( ctx )
617  {
619  return m_source_set;
621  return m_target_set;
623  return m_overlap_set;
625  return m_covering_source_set;
626  case Remapper::DEFAULT:
627  default:
628  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_set );
629  }
630 }
631 
633 {
634  switch( ctx )
635  {
637  return m_source_entities;
639  return m_target_entities;
641  return m_overlap_entities;
644  case Remapper::DEFAULT:
645  default:
646  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_entities );
647  }
648 }
649 
651 {
652  switch( ctx )
653  {
655  return m_source_entities;
657  return m_target_entities;
659  return m_overlap_entities;
662  case Remapper::DEFAULT:
663  default:
664  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_entities );
665  }
666 }
667 
669 {
670  switch( ctx )
671  {
673  return m_source_vertices;
675  return m_target_vertices;
678  case Remapper::DEFAULT:
679  default:
680  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_source_vertices );
681  }
682 }
683 
685 {
686  switch( ctx )
687  {
689  return m_source_vertices;
691  return m_target_vertices;
694  case Remapper::DEFAULT:
695  default:
696  MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_source_vertices );
697  }
698 }
699 
700 inline void TempestRemapper::SetMeshType( Remapper::IntersectionContext ctx, const std::vector< int >& metadata )
701 {
702  switch( ctx )
703  {
705  m_source_type = static_cast< moab::TempestRemapper::TempestMeshType >( metadata[0] );
706  if( metadata[0] == 1 ) // RLL mesh
707  {
708  m_source_metadata.resize( 2 );
709  m_source_metadata[0] = metadata[1];
710  m_source_metadata[1] = metadata[2];
711  }
712  else
713  {
714  m_source_metadata.resize( 1 );
715  m_source_metadata[0] = metadata[1];
716  }
717  break;
719  m_target_type = static_cast< moab::TempestRemapper::TempestMeshType >( metadata[0] );
720  if( metadata[0] == 1 ) // RLL mesh
721  {
722  m_target_metadata.resize( 2 );
723  m_target_metadata[0] = metadata[1];
724  m_target_metadata[1] = metadata[2];
725  }
726  else
727  {
728  m_target_metadata.resize( 1 );
729  m_target_metadata[0] = metadata[1];
730  }
731  break;
734  default:
735  break;
736  }
737 }
738 
740 {
741  switch( ctx )
742  {
744  return m_source_type;
746  return m_target_type;
748  return m_overlap_type;
749  case Remapper::DEFAULT:
750  default:
752  }
753 }
754 
756 {
757  return m_covering_source;
758 }
759 
761 {
762  return m_covering_source_set;
763 }
764 
765 // inline int TempestRemapper::GetGlobalID( Remapper::IntersectionContext ctx, int localID )
766 // {
767 // switch( ctx )
768 // {
769 // case Remapper::SourceMesh:
770 // return lid_to_gid_src[localID];
771 // case Remapper::TargetMesh:
772 // return lid_to_gid_tgt[localID];
773 // case Remapper::CoveringMesh:
774 // return lid_to_gid_covsrc[localID];
775 // case Remapper::OverlapMesh:
776 // case Remapper::DEFAULT:
777 // default:
778 // return -1;
779 // }
780 // }
781 
782 // inline int TempestRemapper::GetLocalID( Remapper::IntersectionContext ctx, int globalID )
783 // {
784 // switch( ctx )
785 // {
786 // case Remapper::SourceMesh:
787 // return gid_to_lid_src[globalID];
788 // case Remapper::TargetMesh:
789 // return gid_to_lid_tgt[globalID];
790 // case Remapper::CoveringMesh:
791 // return gid_to_lid_covsrc[globalID];
792 // case Remapper::DEFAULT:
793 // case Remapper::OverlapMesh:
794 // default:
795 // return -1;
796 // }
797 // }
798 
799 } // namespace moab
800 
801 #endif // MB_TEMPESTREMAPPER_HPP