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