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 TempestOfflineMap;
35 
36 class TempestRemapper : public Remapper
37 {
38 public:
39 #ifdef MOAB_HAVE_MPI
TempestRemapper(moab::Interface * mbInt,moab::ParallelComm * pcomm=NULL)40     TempestRemapper(moab::Interface* mbInt, moab::ParallelComm* pcomm = NULL) :
41         Remapper(mbInt, pcomm),
42 #else
43     TempestRemapper(moab::Interface* mbInt) :
44         Remapper(mbInt),
45 #endif
46         meshValidate(false), constructEdgeMap(false), m_source_type(DEFAULT), m_target_type(DEFAULT)
47     {
48     }
49 
50     virtual ~TempestRemapper();
51 
52     virtual ErrorCode initialize(bool initialize_fsets=true);
53 
54     // Mesh type with a correspondence to Tempest/Climate formats
55     enum TempestMeshType {
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 TempestOfflineMap;
67 
68     moab::ErrorCode GenerateMesh(Remapper::IntersectionContext ctx, TempestMeshType type);
69 
70     moab::ErrorCode LoadMesh(Remapper::IntersectionContext ctx, std::string inputFilename, TempestMeshType type);
71 
72     moab::ErrorCode ComputeOverlapMesh(double tolerance=1e-8, double radius_src=1.0, double radius_tgt=1.0, double boxeps=0.1, bool use_tempest=false);
73 
74     // Converters between MOAB and Tempest representations
75     moab::ErrorCode ConvertTempestMesh(Remapper::IntersectionContext ctx);
76 
77     moab::ErrorCode ConvertMeshToTempest(Remapper::IntersectionContext ctx);
78 
79     Mesh* GetMesh(Remapper::IntersectionContext ctx);
80 
81     void SetMesh(Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite=true);
82 
83     Mesh* GetCoveringMesh();
84 
85     moab::EntityHandle& GetMeshSet(Remapper::IntersectionContext ctx);
86 
87     moab::EntityHandle GetMeshSet(Remapper::IntersectionContext ctx) const;
88 
89     moab::Range& GetMeshEntities(Remapper::IntersectionContext ctx);
90 
91     const moab::Range& GetMeshEntities(Remapper::IntersectionContext ctx) const;
92 
93     moab::EntityHandle& GetCoveringSet();
94 
95     void SetMeshType(Remapper::IntersectionContext ctx, TempestMeshType type);
96 
97     TempestMeshType GetMeshType(Remapper::IntersectionContext ctx) const;
98 
99     int GetGlobalID(Remapper::IntersectionContext ctx, int localID);
100 
101     int GetLocalID(Remapper::IntersectionContext ctx, int globalID);
102 
103     // public members
104     bool meshValidate;  // Validate the mesh after loading from file
105 
106     bool constructEdgeMap;  //  Construct the edge map within the TempestRemap datastructures
107 
108     static const bool verbose = true;
109 
110 private:
111 
112     moab::ErrorCode AssociateSrcTargetInOverlap();
113 
114     moab::ErrorCode ConvertMOABMesh_WithSortedEntitiesBySource();
115 
116     // private methods
117     moab::ErrorCode LoadTempestMesh_Private(std::string inputFilename, Mesh** tempest_mesh);
118 
119     moab::ErrorCode ConvertMOABMeshToTempest_Private(Mesh* mesh, moab::EntityHandle meshset, moab::Range& entities);
120 
121     moab::ErrorCode ConvertTempestMeshToMOAB_Private(TempestMeshType type, Mesh* mesh, moab::EntityHandle& meshset);
122 
123     moab::ErrorCode augment_overlap_set();
124 
125     // Source, Target amd Overlap meshes
126     Mesh* m_source;
127     TempestMeshType m_source_type;
128     moab::Range m_source_entities;
129     moab::EntityHandle m_source_set;
130 
131     Mesh* m_target;
132     TempestMeshType m_target_type;
133     moab::Range m_target_entities;
134     moab::EntityHandle m_target_set;
135 
136     // Overlap meshes
137     Mesh* m_overlap;
138     TempestMeshType m_overlap_type;
139     moab::Range m_overlap_entities;
140     moab::EntityHandle m_overlap_set;
141     std::vector<std::pair<int,int> > m_sorted_overlap_order;
142     // Mesh* m_sorted_overlap;
143 
144     // Parallel - migrated mesh that is in the local view
145     Mesh* m_covering_source;
146     moab::EntityHandle m_covering_source_set;
147     moab::Range m_covering_source_entities;
148 
149     std::map<int,int> gid_to_lid_src, gid_to_lid_covsrc, gid_to_lid_tgt;
150     std::map<int,int> lid_to_gid_src, lid_to_gid_covsrc, lid_to_gid_tgt;
151 
152     bool is_parallel, is_root;
153     int rank, size;
154 };
155 
156 // Inline functions
157 inline
GetMesh(Remapper::IntersectionContext ctx)158 Mesh* TempestRemapper::GetMesh(Remapper::IntersectionContext ctx)
159 {
160     switch(ctx)
161     {
162         case Remapper::SourceMesh:
163             return m_source;
164         case Remapper::TargetMesh:
165             return m_target;
166         case Remapper::IntersectedMesh:
167             return m_overlap;
168         case Remapper::CoveringMesh:
169             return m_covering_source;
170         case Remapper::DEFAULT:
171         default:
172             return NULL;
173     }
174 }
175 
176 inline
SetMesh(Remapper::IntersectionContext ctx,Mesh * mesh,bool overwrite)177 void TempestRemapper::SetMesh(Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite)
178 {
179     switch(ctx)
180     {
181         case Remapper::SourceMesh:
182             if (!overwrite && m_source) return;
183             if (overwrite && m_source) delete m_source;
184             m_source = mesh;
185             break;
186         case Remapper::TargetMesh:
187             if (!overwrite && m_target) return;
188             if (overwrite && m_target) delete m_target;
189             m_target = mesh;
190             break;
191         case Remapper::IntersectedMesh:
192             if (!overwrite && m_overlap) return;
193             if (overwrite && m_overlap) delete m_overlap;
194             m_overlap = mesh;
195             break;
196         case Remapper::CoveringMesh:
197             if (!overwrite && m_covering_source) return;
198             if (overwrite && m_covering_source) delete m_covering_source;
199             m_covering_source = mesh;
200             break;
201         case Remapper::DEFAULT:
202         default:
203             break;
204     }
205 }
206 
207 inline
GetMeshSet(Remapper::IntersectionContext ctx)208 moab::EntityHandle& TempestRemapper::GetMeshSet(Remapper::IntersectionContext ctx)
209 {
210     switch(ctx)
211     {
212         case Remapper::SourceMesh:
213             return m_source_set;
214         case Remapper::TargetMesh:
215             return m_target_set;
216         case Remapper::IntersectedMesh:
217             return m_overlap_set;
218         case Remapper::CoveringMesh:
219             return m_covering_source_set;
220         case Remapper::DEFAULT:
221         default:
222             MB_SET_ERR_RET_VAL("Invalid context passed to GetMeshSet", m_overlap_set);
223     }
224 }
225 
226 inline
GetMeshSet(Remapper::IntersectionContext ctx) const227 moab::EntityHandle TempestRemapper::GetMeshSet(Remapper::IntersectionContext ctx) const
228 {
229     switch(ctx)
230     {
231         case Remapper::SourceMesh:
232             return m_source_set;
233         case Remapper::TargetMesh:
234             return m_target_set;
235         case Remapper::IntersectedMesh:
236             return m_overlap_set;
237         case Remapper::CoveringMesh:
238             return m_covering_source_set;
239         case Remapper::DEFAULT:
240         default:
241             MB_SET_ERR_RET_VAL("Invalid context passed to GetMeshSet", m_overlap_set);
242     }
243 }
244 
245 
246 inline
GetMeshEntities(Remapper::IntersectionContext ctx)247 moab::Range& TempestRemapper::GetMeshEntities(Remapper::IntersectionContext ctx)
248 {
249     switch(ctx)
250     {
251         case Remapper::SourceMesh:
252             return m_source_entities;
253         case Remapper::TargetMesh:
254             return m_target_entities;
255         case Remapper::IntersectedMesh:
256             return m_overlap_entities;
257         case Remapper::CoveringMesh:
258             return m_covering_source_entities;
259         case Remapper::DEFAULT:
260         default:
261             MB_SET_ERR_RET_VAL("Invalid context passed to GetMeshSet", m_overlap_entities);
262     }
263 }
264 
265 inline
GetMeshEntities(Remapper::IntersectionContext ctx) const266 const moab::Range& TempestRemapper::GetMeshEntities(Remapper::IntersectionContext ctx) const
267 {
268     switch(ctx)
269     {
270         case Remapper::SourceMesh:
271             return m_source_entities;
272         case Remapper::TargetMesh:
273             return m_target_entities;
274         case Remapper::IntersectedMesh:
275             return m_overlap_entities;
276         case Remapper::CoveringMesh:
277             return m_covering_source_entities;
278         case Remapper::DEFAULT:
279         default:
280             MB_SET_ERR_RET_VAL("Invalid context passed to GetMeshSet", m_overlap_entities);
281     }
282 }
283 
284 inline
SetMeshType(Remapper::IntersectionContext ctx,TempestRemapper::TempestMeshType type)285 void TempestRemapper::SetMeshType(Remapper::IntersectionContext ctx, TempestRemapper::TempestMeshType type)
286 {
287     switch(ctx)
288     {
289         case Remapper::SourceMesh:
290             m_source_type = type;
291             break;
292         case Remapper::TargetMesh:
293             m_target_type = type;
294             break;
295         case Remapper::IntersectedMesh:
296             m_overlap_type = type;
297             break;
298         case Remapper::DEFAULT:
299         default:
300             break;
301     }
302 }
303 
304 inline
GetMeshType(Remapper::IntersectionContext ctx) const305 TempestRemapper::TempestMeshType TempestRemapper::GetMeshType(Remapper::IntersectionContext ctx) const
306 {
307     switch(ctx)
308     {
309         case Remapper::SourceMesh:
310             return m_source_type;
311         case Remapper::TargetMesh:
312             return m_target_type;
313         case Remapper::IntersectedMesh:
314             return m_overlap_type;
315         case Remapper::DEFAULT:
316         default:
317             return TempestRemapper::DEFAULT;
318     }
319 }
320 
321 inline
GetCoveringMesh()322 Mesh* TempestRemapper::GetCoveringMesh() {
323     return m_covering_source;
324 }
325 
326 inline
GetCoveringSet()327 moab::EntityHandle& TempestRemapper::GetCoveringSet() {
328     return m_covering_source_set;
329 }
330 
331 inline
GetGlobalID(Remapper::IntersectionContext ctx,int localID)332 int TempestRemapper::GetGlobalID(Remapper::IntersectionContext ctx, int localID)
333 {
334     switch(ctx)
335     {
336         case Remapper::SourceMesh:
337             return lid_to_gid_src[localID];
338         case Remapper::TargetMesh:
339             return lid_to_gid_tgt[localID];
340         case Remapper::CoveringMesh:
341             return lid_to_gid_covsrc[localID];
342         case Remapper::IntersectedMesh:
343         case Remapper::DEFAULT:
344         default:
345             return -1;
346     }
347 }
348 
349 inline
GetLocalID(Remapper::IntersectionContext ctx,int globalID)350 int TempestRemapper::GetLocalID(Remapper::IntersectionContext ctx, int globalID)
351 {
352     switch(ctx)
353     {
354         case Remapper::SourceMesh:
355             return gid_to_lid_src[globalID];
356         case Remapper::TargetMesh:
357             return gid_to_lid_tgt[globalID];
358         case Remapper::CoveringMesh:
359             return gid_to_lid_covsrc[globalID];
360         case Remapper::DEFAULT:
361         case Remapper::IntersectedMesh:
362         default:
363             return -1;
364     }
365 }
366 
367 
368 }
369 
370 #endif // MB_TEMPESTREMAPPER_HPP
371