 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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), mahadevan@anl.gov
11 *
12 * =====================================================================================
13 */
14
15 #ifndef MB_TEMPESTREMAPPER_HPP
16 #define MB_TEMPESTREMAPPER_HPP
17
18 #include "moab/Remapping/Remapper.hpp"
19 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
20 #include "moab/IntxMesh/IntxUtils.hpp"
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 ),
47 m_target_type( DEFAULT )
48 {
49 }
50
51 virtual ~TempestRemapper();
52
53 // Mesh type with a correspondence to Tempest/Climate formats
54 enum TempestMeshType
55 {
56 DEFAULT = -1,
57 CS = 0,
58 RLL = 1,
59 ICO = 2,
60 ICOD = 3,
61 OVERLAP_FILES = 4,
62 OVERLAP_MEMORY = 5,
63 OVERLAP_MOAB = 6
64 };
65
66 friend class TempestOnlineMap;
67
68 public:
69 /// <summary>
70 /// Initialize the TempestRemapper object internal datastructures including the mesh sets
71 /// and TempestRemap mesh references.
72 /// </summary>
73 virtual ErrorCode initialize( bool initialize_fsets = true );
74
75 /// <summary>
76 /// Deallocate and clear any memory initialized in the TempestRemapper object
77 /// </summary>
78 virtual ErrorCode clear();
79
80 /// <summary>
81 /// Generate a mesh in memory of given type (CS/RLL/ICO/MPAS(structured)) and store it
82 /// under the context specified by the user.
83 /// </summary>
84 moab::ErrorCode GenerateMesh( Remapper::IntersectionContext ctx, TempestMeshType type );
85
86 /// <summary>
87 /// Load a mesh from disk of given type and store it under the context specified by the
88 /// user.
89 /// </summary>
90 moab::ErrorCode LoadMesh( Remapper::IntersectionContext ctx, std::string inputFilename, TempestMeshType type );
91
92 /// <summary>
93 /// Construct a source covering mesh such that it completely encompasses the target grid in
94 /// parallel. This operation is critical to ensure that the parallel advancing-front
95 /// intersection algorithm can compute the intersection mesh only locally without any process
96 /// communication.
97 /// </summary>
98 moab::ErrorCode ConstructCoveringSet( double tolerance = 1e-8,
99 double radius_src = 1.0,
100 double radius_tgt = 1.0,
101 double boxeps = 0.1,
102 bool regional_mesh = false,
103 bool gnomonic = true,
104 int nb_ghost_layers = 0 );
105
106 /// <summary>
107 /// Compute the intersection mesh between the source and target grids that have been
108 /// instantiated in the Remapper. This function invokes the parallel advancing-front
109 /// intersection algorithm internally for spherical meshes and can handle arbitrary
110 /// unstructured grids (CS, RLL, ICO, MPAS) with and without holes.
111 /// </summary>
112 moab::ErrorCode ComputeOverlapMesh( bool kdtree_search = true, bool use_tempest = false);
113
114 /* Converters between MOAB and Tempest representations */
115
116 /// <summary>
117 /// Convert the TempestRemap mesh object to a corresponding MOAB mesh representation
118 /// according to the intersection context.
119 /// </summary>
120 moab::ErrorCode ConvertTempestMesh( Remapper::IntersectionContext ctx );
121
122 /// <summary>
123 /// Convert the MOAB mesh representation to a corresponding TempestRemap mesh object
124 /// according to the intersection context.
125 /// </summary>
126 moab::ErrorCode ConvertMeshToTempest( Remapper::IntersectionContext ctx );
127
128 /// <summary>
129 /// Get the TempestRemap mesh object according to the intersection context.
130 /// </summary>
131 Mesh* GetMesh( Remapper::IntersectionContext ctx );
132
133 /// <summary>
134 /// Set the TempestRemap mesh object according to the intersection context.
135 /// </summary>
136 void SetMesh( Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite = true );
137
138 void SetMeshSet( Remapper::IntersectionContext ctx /* Remapper::CoveringMesh*/,
139 moab::EntityHandle mset,
140 moab::Range* entities = nullptr );
141 /// <summary>
142 /// Get the covering mesh (TempestRemap) object.
143 /// </summary>
144 Mesh* GetCoveringMesh();
145
146 /// <summary>
147 /// Get the MOAB mesh set corresponding to the intersection context.
148 /// </summary>
149 moab::EntityHandle& GetMeshSet( Remapper::IntersectionContext ctx );
150
151 /// <summary>
152 /// Const overload. Get the MOAB mesh set corresponding to the intersection context.
153 /// </summary>
154 moab::EntityHandle GetMeshSet( Remapper::IntersectionContext ctx ) const;
155
156 /// <summary>
157 /// Get the mesh element entities corresponding to the intersection context.
158 /// </summary>
159 moab::Range& GetMeshEntities( Remapper::IntersectionContext ctx );
160
161 /// <summary>
162 /// Const overload. Get the mesh element entities corresponding to the intersection context.
163 /// </summary>
164 const moab::Range& GetMeshEntities( Remapper::IntersectionContext ctx ) const;
165
166 /// <summary>
167 /// Get the mesh vertices corresponding to the intersection context. Useful for point-cloud
168 /// meshes
169 /// </summary>
170 moab::Range& GetMeshVertices( Remapper::IntersectionContext ctx );
171
172 /// <summary>
173 /// Const overload. Get the mesh vertices corresponding to the intersection context. Useful
174 /// for point-cloud meshes.
175 /// </summary>
176 const moab::Range& GetMeshVertices( Remapper::IntersectionContext ctx ) const;
177
178 /// <summary>
179 /// Get access to the underlying source covering set if available. Else return the source
180 /// set.
181 /// </summary>
182 moab::EntityHandle& GetCoveringSet();
183
184 /// <summary>
185 /// Set the mesh type corresponding to the intersection context
186 /// </summary>
187 void SetMeshType( Remapper::IntersectionContext ctx, const std::vector< int >& metadata );
188
189 /// <summary>
190 /// reconstruct mesh, used now only for IO; need a better solution maybe
191 /// </summary>
192 void ResetMeshSet( Remapper::IntersectionContext ctx, moab::EntityHandle meshSet );
193
194 /// <summary>
195 /// Get the mesh type corresponding to the intersection context
196 /// </summary>
197 TempestMeshType GetMeshType( Remapper::IntersectionContext ctx ) const;
198
199 /// <summary>
200 /// Get the global ID corresponding to the local entity ID according to the context (source,
201 /// target, intersection)
202 /// </summary>
203 int GetGlobalID( Remapper::IntersectionContext ctx, int localID );
204
205 /// <summary>
206 /// Get the local ID corresponding to the global entity ID according to the context (source,
207 /// target, intersection)
208 /// </summary>
209 int GetLocalID( Remapper::IntersectionContext ctx, int globalID );
210
211 /// <summary>
212 /// Gather the overlap mesh and asssociated source/target data and write it out to disk
213 /// using the TempestRemap output interface. This information can then be used with the
214 /// "GenerateOfflineMap" tool in TempestRemap as needed.
215 /// </summary>
216 moab::ErrorCode WriteTempestIntersectionMesh( std::string strOutputFileName,
217 const bool fAllParallel,
218 const bool fInputConcave,
219 const bool fOutputConcave );
220
221 /// <summary>
222 /// Generate the necessary metadata and specifically the GLL node numbering for DoFs for a
223 /// CS mesh. This negates the need for running external code like HOMME to output the
224 /// numbering needed for computing maps. The functionality is used through the `mbconvert`
225 /// tool to compute processor-invariant Global DoF IDs at GLL nodes.
226 /// </summary>
227 moab::ErrorCode GenerateCSMeshMetadata( const int ntot_elements,
228 moab::Range& entities,
229 moab::Range* secondary_entities,
230 const std::string& dofTagName,
231 int nP );
232
233 /// <summary>
234 /// Generate the necessary metadata for DoF node numbering in a given mesh.
235 /// Currently, only the functionality to generate numbering on CS grids is supported.
236 /// </summary>
237 moab::ErrorCode GenerateMeshMetadata( Mesh& mesh,
238 const int ntot_elements,
239 moab::Range& entities,
240 moab::Range* secondary_entities,
241 const std::string dofTagName,
242 int nP );
243
244 /// <summary>
245 /// Compute the local and global IDs for elements in source/target/coverage meshes.
246 /// </summary>
247 moab::ErrorCode ComputeGlobalLocalMaps();
248
249 /// <summary>
250 /// Get all the ghosted overlap entities that were accumulated to enable conservation in
251 /// parallel
252 /// </summary>
253 moab::ErrorCode GetOverlapAugmentedEntities( moab::Range& sharedGhostEntities );
254
255 #ifndef MOAB_HAVE_MPI
256 /// <summary>
257 /// Internal method to assign vertex and element global IDs if one does not exist already
258 /// </summary>
259 moab::ErrorCode assign_vertex_element_IDs( Tag idtag,
260 EntityHandle this_set,
261 const int dimension = 2,
262 const int start_id = 1 );
263 #endif
264 // <summary>
265 /// Get the masks that could have been defined
266 /// </summary>
267 ErrorCode GetIMasks( Remapper::IntersectionContext ctx, std::vector< int >& masks );
268
269 public: // public members
270 const bool offlineWorkflow; // check whether we are in an offline workflow context (mbtempest)
271 bool meshValidate; // Validate the mesh after loading from file
272
273 bool constructEdgeMap; // Construct the edge map within the TempestRemap datastructures
274
275 static const bool verbose = true;
276
277 private:
278 moab::ErrorCode convert_overlap_mesh_sorted_by_source();
279
280 // private methods
281 moab::ErrorCode load_tempest_mesh_private( std::string inputFilename, Mesh** tempest_mesh );
282
283 moab::ErrorCode convert_mesh_to_tempest_private( Mesh* mesh,
284 moab::EntityHandle meshset,
285 moab::Range& entities,
286 moab::Range* pverts );
287
288 moab::ErrorCode convert_tempest_mesh_private( TempestMeshType type,
289 Mesh* mesh,
290 moab::EntityHandle& meshset,
291 moab::Range& entities,
292 moab::Range* vertices );
293
294 moab::ErrorCode augment_overlap_set();
295
296 /* Source meshset, mesh and entity references */
297 Mesh* m_source;
298 TempestMeshType m_source_type;
299 moab::Range m_source_entities;
300 moab::Range m_source_vertices;
301 moab::EntityHandle m_source_set;
302 int max_source_edges;
303 bool point_cloud_source;
304 std::vector< int > m_source_metadata;
305
306 /* Target meshset, mesh and entity references */
307 Mesh* m_target;
308 TempestMeshType m_target_type;
309 moab::Range m_target_entities;
310 moab::Range m_target_vertices;
311 moab::EntityHandle m_target_set;
312 int max_target_edges;
313 bool point_cloud_target;
314 std::vector< int > m_target_metadata;
315
316 /* Overlap meshset, mesh and entity references */
317 Mesh* m_overlap;
318 TempestMeshType m_overlap_type;
319 moab::Range m_overlap_entities;
320 moab::EntityHandle m_overlap_set;
321 std::vector< std::pair< int, int > > m_sorted_overlap_order;
322
323 /* Intersection context on a sphere */
324 moab::Intx2MeshOnSphere* mbintx;
325
326 /* Parallel - migrated mesh that is in the local view */
327 Mesh* m_covering_source;
328 moab::EntityHandle m_covering_source_set;
329 moab::Range m_covering_source_entities;
330 moab::Range m_covering_source_vertices;
331
332 /* local to glboal and global to local ID maps */
333 std::map< int, int > gid_to_lid_src, gid_to_lid_covsrc, gid_to_lid_tgt;
334 std::map< int, int > lid_to_gid_src, lid_to_gid_covsrc, lid_to_gid_tgt;
335
336 IntxAreaUtils::AreaMethod m_area_method;
337
338 bool rrmgrids;
339 bool is_parallel, is_root;
340 int rank, size;
341 };
342
343 // Inline functions
344 inline Mesh* TempestRemapper::GetMesh( Remapper::IntersectionContext ctx )
345 {
346 switch( ctx )
347 {
348 case Remapper::SourceMesh:
349 return m_source;
350 case Remapper::TargetMesh:
351 return m_target;
352 case Remapper::OverlapMesh:
353 return m_overlap;
354 case Remapper::CoveringMesh:
355 return m_covering_source;
356 case Remapper::DEFAULT:
357 default:
358 return NULL;
359 }
360 }
361
362 inline void TempestRemapper::SetMesh( Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite )
363 {
364 switch( ctx )
365 {
366 case Remapper::SourceMesh:
367 if( !overwrite && m_source ) return;
368 if( overwrite && m_source ) delete m_source;
369 m_source = mesh;
370 break;
371 case Remapper::TargetMesh:
372 if( !overwrite && m_target ) return;
373 if( overwrite && m_target ) delete m_target;
374 m_target = mesh;
375 break;
376 case Remapper::OverlapMesh:
377 if( !overwrite && m_overlap ) return;
378 if( overwrite && m_overlap ) delete m_overlap;
379 m_overlap = mesh;
380 break;
381 case Remapper::CoveringMesh:
382 if( !overwrite && m_covering_source ) return;
383 if( overwrite && m_covering_source ) delete m_covering_source;
384 m_covering_source = mesh;
385 break;
386 case Remapper::DEFAULT:
387 default:
388 break;
389 }
390 }
391
392 inline void TempestRemapper::ResetMeshSet( Remapper::IntersectionContext ctx, moab::EntityHandle meshSet )
393 {
394 switch( ctx )
395 {
396 case Remapper::SourceMesh:
397 delete m_source;
398 m_source = new Mesh;
399 m_source_set = meshSet;
400 convert_mesh_to_tempest_private( m_source, m_source_set, m_source_entities, &m_source_vertices );
401 m_source->CalculateFaceAreas( false ); // fInputConcave is false ?
402 break;
403 case Remapper::TargetMesh:
404 // not needed yet
405 break;
406 case Remapper::OverlapMesh:
407 // not needed yet
408 break;
409 case Remapper::CoveringMesh:
410 // not needed yet
411 break;
412 case Remapper::DEFAULT:
413 default:
414 break;
415 }
416 }
417
418
419 inline moab::EntityHandle& TempestRemapper::GetMeshSet( Remapper::IntersectionContext ctx )
420 {
421 switch( ctx )
422 {
423 case Remapper::SourceMesh:
424 return m_source_set;
425 case Remapper::TargetMesh:
426 return m_target_set;
427 case Remapper::OverlapMesh:
428 return m_overlap_set;
429 case Remapper::CoveringMesh:
430 return m_covering_source_set;
431 case Remapper::DEFAULT:
432 default:
433 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_set );
434 }
435 }
436
437 inline moab::EntityHandle TempestRemapper::GetMeshSet( Remapper::IntersectionContext ctx ) const
438 {
439 switch( ctx )
440 {
441 case Remapper::SourceMesh:
442 return m_source_set;
443 case Remapper::TargetMesh:
444 return m_target_set;
445 case Remapper::OverlapMesh:
446 return m_overlap_set;
447 case Remapper::CoveringMesh:
448 return m_covering_source_set;
449 case Remapper::DEFAULT:
450 default:
451 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_set );
452 }
453 }
454
455 inline moab::Range& TempestRemapper::GetMeshEntities( Remapper::IntersectionContext ctx )
456 {
457 switch( ctx )
458 {
459 case Remapper::SourceMesh:
460 return m_source_entities;
461 case Remapper::TargetMesh:
462 return m_target_entities;
463 case Remapper::OverlapMesh:
464 return m_overlap_entities;
465 case Remapper::CoveringMesh:
466 return m_covering_source_entities;
467 case Remapper::DEFAULT:
468 default:
469 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_entities );
470 }
471 }
472
473 inline const moab::Range& TempestRemapper::GetMeshEntities( Remapper::IntersectionContext ctx ) const
474 {
475 switch( ctx )
476 {
477 case Remapper::SourceMesh:
478 return m_source_entities;
479 case Remapper::TargetMesh:
480 return m_target_entities;
481 case Remapper::OverlapMesh:
482 return m_overlap_entities;
483 case Remapper::CoveringMesh:
484 return m_covering_source_entities;
485 case Remapper::DEFAULT:
486 default:
487 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_entities );
488 }
489 }
490
491 inline moab::Range& TempestRemapper::GetMeshVertices( Remapper::IntersectionContext ctx )
492 {
493 switch( ctx )
494 {
495 case Remapper::SourceMesh:
496 return m_source_vertices;
497 case Remapper::TargetMesh:
498 return m_target_vertices;
499 case Remapper::CoveringMesh:
500 return m_covering_source_vertices;
501 case Remapper::DEFAULT:
502 default:
503 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_source_vertices );
504 }
505 }
506
507 inline const moab::Range& TempestRemapper::GetMeshVertices( Remapper::IntersectionContext ctx ) const
508 {
509 switch( ctx )
510 {
511 case Remapper::SourceMesh:
512 return m_source_vertices;
513 case Remapper::TargetMesh:
514 return m_target_vertices;
515 case Remapper::CoveringMesh:
516 return m_covering_source_vertices;
517 case Remapper::DEFAULT:
518 default:
519 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_source_vertices );
520 }
521 }
522
523 inline void TempestRemapper::SetMeshType( Remapper::IntersectionContext ctx, const std::vector< int >& metadata )
524 {
525 switch( ctx )
526 {
527 case Remapper::SourceMesh:
528 m_source_type = static_cast< moab::TempestRemapper::TempestMeshType >( metadata[0] );
529 if( metadata[0] == 1 ) // RLL mesh
530 {
531 m_source_metadata.resize( 2 );
532 m_source_metadata[0] = metadata[1];
533 m_source_metadata[1] = metadata[2];
534 }
535 else
536 {
537 m_source_metadata.resize( 1 );
538 m_source_metadata[0] = metadata[1];
539 }
540 break;
541 case Remapper::TargetMesh:
542 m_target_type = static_cast< moab::TempestRemapper::TempestMeshType >( metadata[0] );
543 if( metadata[0] == 1 ) // RLL mesh
544 {
545 m_target_metadata.resize( 2 );
546 m_target_metadata[0] = metadata[1];
547 m_target_metadata[1] = metadata[2];
548 }
549 else
550 {
551 m_target_metadata.resize( 1 );
552 m_target_metadata[0] = metadata[1];
553 }
554 break;
555 case Remapper::OverlapMesh:
556 m_overlap_type = moab::TempestRemapper::OVERLAP_FILES;
557 default:
558 break;
559 }
560 }
561
562 inline TempestRemapper::TempestMeshType TempestRemapper::GetMeshType( Remapper::IntersectionContext ctx ) const
563 {
564 switch( ctx )
565 {
566 case Remapper::SourceMesh:
567 return m_source_type;
568 case Remapper::TargetMesh:
569 return m_target_type;
570 case Remapper::OverlapMesh:
571 return m_overlap_type;
572 case Remapper::DEFAULT:
573 default:
574 return TempestRemapper::DEFAULT;
575 }
576 }
577
578 inline Mesh* TempestRemapper::GetCoveringMesh()
579 {
580 return m_covering_source;
581 }
582
583 inline moab::EntityHandle& TempestRemapper::GetCoveringSet()
584 {
585 return m_covering_source_set;
586 }
587
588 inline int TempestRemapper::GetGlobalID( Remapper::IntersectionContext ctx, int localID )
589 {
590 switch( ctx )
591 {
592 case Remapper::SourceMesh:
593 return lid_to_gid_src[localID];
594 case Remapper::TargetMesh:
595 return lid_to_gid_tgt[localID];
596 case Remapper::CoveringMesh:
597 return lid_to_gid_covsrc[localID];
598 case Remapper::OverlapMesh:
599 case Remapper::DEFAULT:
600 default:
601 return -1;
602 }
603 }
604
605 inline int TempestRemapper::GetLocalID( Remapper::IntersectionContext ctx, int globalID )
606 {
607 switch( ctx )
608 {
609 case Remapper::SourceMesh:
610 return gid_to_lid_src[globalID];
611 case Remapper::TargetMesh:
612 return gid_to_lid_tgt[globalID];
613 case Remapper::CoveringMesh:
614 return gid_to_lid_covsrc[globalID];
615 case Remapper::DEFAULT:
616 case Remapper::OverlapMesh:
617 default:
618 return -1;
619 }
620 }
621
622 } // namespace moab
623
624 #endif // MB_TEMPESTREMAPPER_HPP
1.9.1.