 |
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 * Usage: MOAB-Tempest tool
3 *
4 * Generate a Cubed-Sphere mesh: ./mbtempest -t 0 -res 25 -f cubed_sphere_mesh.exo
5 * Generate a RLL mesh: ./mbtempest -t 1 -res 25 -f rll_mesh.exo
6 * Generate a Icosahedral-Sphere mesh: ./mbtempest -t 2 -res 25 <-dual> -f icosa_mesh.exo
7 *
8 * Now you can compute the intersections between the meshes too!
9 *
10 * Generate the overlap mesh: ./mbtempest -t 3 -l cubed_sphere_mesh.exo -l rll_mesh.exo -f
11 * overlap_mesh.exo
12 *
13 */
14
15 // standard C++ includes
16 #include <iostream>
17 #include <iomanip>
18 #include <cstdlib>
19 #include <vector>
20 #include <string>
21 #include <sstream>
22 #include <cassert>
23
24 // MOAB includes
25 #include "moab/Core.hpp"
26 #include "moab/IntxMesh/IntxUtils.hpp"
27 #include "moab/Remapping/TempestRemapper.hpp"
28 #include "moab/Remapping/TempestOnlineMap.hpp"
29 #include "moab/ProgOptions.hpp"
30 #include "moab/CpuTimer.hpp"
31 #include "DebugOutput.hpp"
32
33 //#ifndef MOAB_HAVE_MPI
34 // #error mbtempest tool requires MPI configuration
35 //#endif
36
37 #ifdef MOAB_HAVE_MPI
38 // MPI includes
39 #include "moab_mpi.h"
40 #include "moab/ParallelComm.hpp"
41 #include "MBParallelConventions.h"
42 #endif
43
44 struct ToolContext
45 {
46 moab::Core* mbcore;
47 #ifdef MOAB_HAVE_MPI
48 moab::ParallelComm* pcomm;
49 #endif
50 const int proc_id, n_procs; // MPI process id and number of processes
51 moab::DebugOutput outputFormatter; // Output formatter
52 int blockSize; // Resolution of the mesh
53 int nlayers; // Number of ghost layers in the source mesh
54 std::vector< std::string > inFilenames; // Input mesh filenames
55 std::vector< Mesh* > meshes; // Mesh references in TempestRemap format
56 std::vector< moab::EntityHandle > meshsets; // Meshsets in MOAB format
57 std::vector< int > disc_orders; // Discretization orders
58 std::vector< std::string > disc_methods; // Discretization method names (fv, cgll, dgll, etc.)
59 std::vector< std::string > doftag_names; // DoF tag names (GLOBAL_ID, etc.)
60 std::string fvMethod; // FV method name (invdist, delaunay, bilin, intbilin, intbilingb, none)
61 std::string outFilename; // Output filename
62 std::string intxFilename; // Output intersection mesh filename
63 std::string baselineFile; // Output baseline file (for verification)
64 std::string variableToVerify; // Variable name for verification
65 moab::TempestRemapper::TempestMeshType
66 meshType; // Mesh type (CS, RLL, ICO, OVERLAP_FILES, OVERLAP_MEMORY, OVERLAP_MOAB)
67 bool skip_io; // Skip I/O operations; strictly to measure the performance of the intersection/map/SpMV algorithms
68 bool computeDual; // Compute the dual of the mesh
69 bool computeWeights; // Compute the projection weights
70 bool verifyWeights; // Verify the projection weights
71 bool enforceConvexity; // Enforce convexity of the input meshes
72 int ensureMonotonicity; // Ensure monotonicity in the weight generation process
73 bool rrmGrids; // Flag specifying that we are dealing with regionally refined grids (possibly nonoverlapping)
74 bool kdtreeSearch; // Use Kd-tree based search for computing mesh intersections instead of advancing front
75 bool fCheck; // Check the generated map for conservation and consistency
76 bool fVolumetric; // Apply a volumetric projection to compute the weights
77 bool useGnomonicProjection; // Use Gnomonic plane projections to compute coverage mesh
78 moab::TempestOnlineMap::CAASType cassType; // CAAS filter type
79 GenerateOfflineMapAlgorithmOptions mapOptions; // Offline map options
80 bool print_diagnostics; // Print diagnostics
81 double boxeps; // Box error tolerance
82 double epsrel; // Relative error tolerance
83
84 #ifdef MOAB_HAVE_MPI
85 ToolContext( moab::Core* icore, moab::ParallelComm* p_pcomm )
86 : mbcore( icore ), pcomm( p_pcomm ), proc_id( pcomm->rank() ), n_procs( pcomm->size() ),
87 outputFormatter( std::cout, pcomm->rank(), 0 ),
88 #else
89 ToolContext( moab::Core* icore )
90 : mbcore( icore ), proc_id( 0 ), n_procs( 1 ), outputFormatter( std::cout, 0, 0 ),
91 #endif
92 blockSize( 5 ), nlayers( 0 ), fvMethod( "none" ), outFilename( "outputFile.nc" ), intxFilename( "" ),
93 baselineFile( "" ), variableToVerify( "" ), meshType( moab::TempestRemapper::DEFAULT ), skip_io( false ),
94 computeDual( false ), computeWeights( false ), verifyWeights( false ), enforceConvexity( false ),
95 ensureMonotonicity( 0 ), rrmGrids( false ), kdtreeSearch( true ), fCheck( false ), fVolumetric( false ),
96 useGnomonicProjection( false ), cassType( moab::TempestOnlineMap::CAAS_NONE ), print_diagnostics( false ),
97 boxeps( 1e-7 ), // Box error tolerance default value
98 epsrel( ReferenceTolerance ) // ReferenceTolerance is defined in Defines.h in TempestRemap
99 {
100 inFilenames.resize( 2 );
101 doftag_names.resize( 2 );
102 timer = new moab::CpuTimer();
103
104 outputFormatter.set_prefix( "[MBTempest]: " );
105 }
106
107 ~ToolContext()
108 {
109 // for (unsigned i=0; i < meshes.size(); ++i) delete meshes[i];
110 meshes.clear();
111 inFilenames.clear();
112 disc_orders.clear();
113 disc_methods.clear();
114 doftag_names.clear();
115 outFilename.clear();
116 intxFilename.clear();
117 baselineFile.clear();
118 variableToVerify.clear();
119 meshsets.clear();
120 delete timer;
121 }
122
123 void timer_push( std::string operation )
124 {
125 timer_ops = timer->time_since_birth();
126 opName = operation;
127 }
128
129 void timer_pop()
130 {
131 double locElapsed = timer->time_since_birth() - timer_ops, avgElapsed = 0, maxElapsed = 0;
132 #ifdef MOAB_HAVE_MPI
133 MPI_Reduce( &locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm() );
134 MPI_Reduce( &locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, pcomm->comm() );
135 #else
136 maxElapsed = locElapsed;
137 avgElapsed = locElapsed;
138 #endif
139 if( !proc_id )
140 {
141 avgElapsed /= n_procs;
142 std::cout << "[LOG] Time taken to " << opName.c_str() << ": max = " << maxElapsed
143 << ", avg = " << avgElapsed << "\n";
144 }
145 // std::cout << "\n[LOG" << proc_id << "] Time taken to " << opName << " = " <<
146 // timer->time_since_birth() - timer_ops << std::endl;
147 opName.clear();
148 }
149
150 void ParseCLOptions( int argc, char* argv[] )
151 {
152 ProgOptions opts;
153 int imeshType = 0;
154 std::string expectedFName = "output.exo";
155 std::string expectedMethod = "fv";
156 std::string expectedFVMethod = "none";
157 std::string expectedDofTagName = "GLOBAL_ID";
158 int expectedOrder = 1;
159 int useCAAS = 0;
160 int nlayer_input = 0;
161
162 if( !proc_id )
163 {
164 std::cout << "Command line options provided to mbtempest:\n ";
165 for( int iarg = 0; iarg < argc; ++iarg )
166 std::cout << argv[iarg] << " ";
167 std::cout << std::endl << std::endl;
168 }
169
170 opts.addOpt< int >( "type,t",
171 "Type of mesh (default=CS; Choose from [CS=0, RLL=1, ICO=2, OVERLAP_FILES=3, "
172 "OVERLAP_MEMORY=4, OVERLAP_MOAB=5])",
173 &imeshType );
174
175 opts.addOpt< int >( "res,r", "Resolution of the mesh (default=5)", &blockSize );
176
177 opts.addOpt< void >( "dual,d", "Output the dual of the mesh (relevant only for ICO mesh type)", &computeDual );
178
179 opts.addOpt< std::string >( "file,f", "Output computed mesh or remapping weights to specified filename",
180 &outFilename );
181 opts.addOpt< std::string >(
182 "load,l", "Input mesh filenames for source and target meshes. (relevant only when computing weights)",
183 &expectedFName );
184
185 opts.addOpt< void >( "advfront,a",
186 "Use the advancing front intersection instead of the Kd-tree based algorithm "
187 "to compute mesh intersections. (relevant only when computing weights)" );
188
189 opts.addOpt< std::string >( "intx,i", "Output TempestRemap intersection mesh filename", &intxFilename );
190
191 opts.addOpt< void >( "weights,w",
192 "Compute and output the weights using the overlap mesh (generally "
193 "relevant only for OVERLAP mesh)",
194 &computeWeights );
195
196 opts.addOpt< void >(
197 "verbose,v", "Print verbose diagnostic messages during intersection and map computation (default=false)",
198 &print_diagnostics );
199
200 opts.addOpt< std::string >( "method,m", "Discretization method for the source and target solution fields",
201 &expectedMethod );
202
203 opts.addOpt< int >( "order,o", "Discretization orders for the source and target solution fields",
204 &expectedOrder );
205
206 opts.addOpt< std::string >( "global_id,g",
207 "Tag name that contains the global DoF IDs for source and target solution fields",
208 &expectedDofTagName );
209
210 opts.addOpt< std::string >( "fvmethod",
211 "Sub-type method for FV-FV projections (invdist, delaunay, bilin, "
212 "intbilin, intbilingb, none. Default: none)",
213 &expectedFVMethod );
214
215 opts.addOpt< void >( "noconserve",
216 "Do not apply conservation to the resultant weights (relevant only "
217 "when computing weights)",
218 &mapOptions.fNoConservation );
219
220 opts.addOpt< void >( "volumetric",
221 "Apply a volumetric projection to compute the weights (relevant only "
222 "when computing weights)",
223 &fVolumetric );
224
225 opts.addOpt< void >( "skip_output", "For performance studies, skip all I/O operations.", &skip_io );
226
227 opts.addOpt< void >( "gnomonic", "Use Gnomonic plane projections to compute coverage mesh.",
228 &useGnomonicProjection );
229
230 opts.addOpt< void >( "enforce_convexity", "check convexity of input meshes to compute mesh intersections",
231 &enforceConvexity );
232
233 opts.addOpt< void >( "nobubble", "do not use bubble on interior of spectral element nodes",
234 &mapOptions.fNoBubble );
235
236 opts.addOpt< void >(
237 "sparseconstraints",
238 "Use sparse solver for constraints when we have high-valence (typical with high-res RLL mesh)",
239 &mapOptions.fSparseConstraints );
240
241 opts.addOpt< void >( "rrmgrids",
242 "At least one of the meshes is a regionally refined grid (relevant to "
243 "accelerate intersection computation)",
244 &rrmGrids );
245
246 opts.addOpt< void >( "checkmap", "Check the generated map for conservation and consistency", &fCheck );
247
248 opts.addOpt< void >( "verify",
249 "Verify the accuracy of the maps by projecting analytical functions "
250 "from source to target "
251 "grid by applying the maps",
252 &verifyWeights );
253
254 opts.addOpt< std::string >( "var",
255 "Tag name of the variable to use in the verification study "
256 "(error metrics for user defined variables may not be available)",
257 &variableToVerify );
258
259 opts.addOpt< int >( "monotonicity", "Ensure monotonicity in the weight generation. Options=[0,1,2,3]",
260 &ensureMonotonicity );
261
262 opts.addOpt< int >( "ghost", "Number of ghost layers in coverage mesh (default=1)", &nlayer_input );
263
264 opts.addOpt< double >( "boxeps", "The tolerance for boxes (default=1e-7)", &boxeps );
265
266 opts.addOpt< int >( "caas", "apply CAAS nonlinear filter after linear map application", &useCAAS );
267
268 opts.addOpt< std::string >( "baseline", "Output baseline file", &baselineFile );
269
270 opts.parseCommandLine( argc, argv );
271
272 // By default - use Kd-tree based search; if user asks for advancing front, disable Kd-tree
273 // algorithm
274 kdtreeSearch = opts.numOptSet( "advfront,a" ) == 0;
275
276 switch( imeshType )
277 {
278 case 0:
279 meshType = moab::TempestRemapper::CS;
280 break;
281
282 case 1:
283 meshType = moab::TempestRemapper::RLL;
284 break;
285
286 case 2:
287 meshType = moab::TempestRemapper::ICO;
288 break;
289
290 case 3:
291 meshType = moab::TempestRemapper::OVERLAP_FILES;
292 break;
293
294 case 4:
295 meshType = moab::TempestRemapper::OVERLAP_MEMORY;
296 break;
297
298 case 5:
299 meshType = moab::TempestRemapper::OVERLAP_MOAB;
300 break;
301
302 default:
303 meshType = moab::TempestRemapper::DEFAULT;
304 break;
305 }
306
307 // decipher whether we want to use CAAS filter when applying the projection
308 switch( useCAAS )
309 {
310 case moab::TempestOnlineMap::CAAS_GLOBAL:
311 cassType = moab::TempestOnlineMap::CAAS_GLOBAL;
312 break;
313 case moab::TempestOnlineMap::CAAS_LOCAL:
314 cassType = moab::TempestOnlineMap::CAAS_LOCAL;
315 break;
316 case moab::TempestOnlineMap::CAAS_LOCAL_ADJACENT:
317 cassType = moab::TempestOnlineMap::CAAS_LOCAL_ADJACENT;
318 break;
319 case moab::TempestOnlineMap::CAAS_QLT:
320 cassType = moab::TempestOnlineMap::CAAS_QLT;
321 break;
322 default:
323 cassType = moab::TempestOnlineMap::CAAS_NONE;
324 break;
325 }
326
327 if( meshType > moab::TempestRemapper::ICO ) // compute overlap mesh and maps possibly
328 {
329 opts.getOptAllArgs( "load,l", inFilenames );
330 opts.getOptAllArgs( "order,o", disc_orders );
331 opts.getOptAllArgs( "method,m", disc_methods );
332 opts.getOptAllArgs( "global_id,i", doftag_names );
333
334 assert( inFilenames.size() == 2 );
335 assert( ensureMonotonicity >= 0 && ensureMonotonicity <= 3 );
336
337 // get discretization order parameters
338 if( disc_orders.size() == 0 ) disc_orders.resize( 2, 1 );
339 if( disc_orders.size() == 1 ) disc_orders.push_back( disc_orders[0] );
340
341 // get discretization method parameters
342 if( disc_methods.size() == 0 ) disc_methods.resize( 2, "fv" );
343 if( disc_methods.size() == 1 ) disc_methods.push_back( disc_methods[0] );
344
345 // get DoF tagname parameters
346 if( doftag_names.size() == 0 ) doftag_names.resize( 2, "GLOBAL_ID" );
347 if( doftag_names.size() == 1 ) doftag_names.push_back( doftag_names[0] );
348
349 assert( disc_orders.size() == 2 );
350 assert( disc_methods.size() == 2 );
351 assert( doftag_names.size() == 2 );
352
353 // for computing maps and overlaps, set discretization orders
354 mapOptions.nPin = disc_orders[0];
355 mapOptions.nPout = disc_orders[1];
356 mapOptions.fSourceConcave = false;
357 mapOptions.fTargetConcave = false;
358
359 mapOptions.strMethod = "";
360
361 if( expectedFVMethod != "none" )
362 {
363 mapOptions.strMethod += expectedFVMethod + ";";
364 fvMethod = expectedFVMethod;
365
366 // These FV projection methods are non-conservative; specify it explicitly
367 mapOptions.fNoConservation = true;
368 }
369 switch( ensureMonotonicity )
370 {
371 case 0:
372 mapOptions.fMonotone = false;
373 break;
374 case 3:
375 mapOptions.strMethod += "mono3;";
376 break;
377 case 2:
378 mapOptions.strMethod += "mono2;";
379 break;
380 default:
381 mapOptions.fMonotone = true;
382 }
383 mapOptions.fNoCorrectAreas = false;
384 mapOptions.fNoCheck = !fCheck;
385
386 //assert( fVolumetric && fInverseDistanceMap == false ); // both options cannot be active
387 if( fVolumetric ) mapOptions.strMethod += "volumetric;";
388
389 // For global meshes, this default should work out of the box.
390 if( !fvMethod.compare( "bilin" ) )
391 nlayers = 3;
392 else
393 nlayers = ( mapOptions.nPin > 1 ? mapOptions.nPin + 1 : 0 );
394 if( nlayer_input ) nlayers = std::max( nlayer_input, nlayers );
395 }
396
397 // clear temporary string name
398 expectedFName.clear();
399
400 mapOptions.strOutputMapFile = outFilename;
401 mapOptions.strOutputFormat = "Netcdf4";
402 }
403
404 private:
405 moab::CpuTimer* timer;
406 double timer_ops;
407 std::string opName;
408 };
409
410 // Forward declare some methods
411 static moab::ErrorCode CreateTempestMesh( ToolContext&, moab::TempestRemapper& remapper, Mesh* );
412 inline double sample_slow_harmonic( double dLon, double dLat );
413 inline double sample_fast_harmonic( double dLon, double dLat );
414 inline double sample_constant( double dLon, double dLat );
415 inline double sample_stationary_vortex( double dLon, double dLat );
416
417 std::string get_file_read_options( ToolContext& ctx, std::string filename )
418 {
419 // if running in serial, blank option may suffice in general
420 std::string opts = "";
421 if( ctx.n_procs > 1 )
422 {
423 size_t lastindex = filename.find_last_of( "." );
424 std::string extension = filename.substr( lastindex + 1, filename.size() );
425 if( extension == "h5m" ) return "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
426 // "PARALLEL_GHOSTS=2.0.3;PARALLEL_THIN_GHOST_LAYER;";
427 else if( extension == "nc" )
428 {
429 // set default set of options
430 opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;";
431
432 if( ctx.proc_id == 0 )
433 {
434 {
435 NcFile ncFile( filename.c_str(), NcFile::ReadOnly );
436 if( !ncFile.is_valid() )
437 _EXCEPTION1( "Unable to open grid file \"%s\" for reading", filename.c_str() );
438
439 // Check for dimension names "grid_size", "grid_rank" and "grid_corners"
440 int iSCRIPFormat = 0, iMPASFormat = 0, iESMFFormat = 0;
441 for( int i = 0; i < ncFile.num_dims(); i++ )
442 {
443 std::string strDimName = ncFile.get_dim( i )->name();
444 if( strDimName == "grid_size" || strDimName == "grid_corners" || strDimName == "grid_rank" )
445 iSCRIPFormat++;
446 if( strDimName == "nodeCount" || strDimName == "elementCount" ||
447 strDimName == "maxNodePElement" )
448 iESMFFormat++;
449 if( strDimName == "nCells" || strDimName == "nEdges" || strDimName == "nVertices" ||
450 strDimName == "vertexDegree" )
451 iMPASFormat++;
452 }
453 ncFile.close();
454
455 if( iESMFFormat == 3 ) // Input from a NetCDF ESMF file
456 {
457 opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;PARALLEL_RESOLVE_SHARED_ENTS;VARIABLE=;";
458 }
459 if( iSCRIPFormat == 3 ) // Input from a NetCDF SCRIP file
460 {
461 opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;";
462 }
463 if( iMPASFormat == 4 ) // Input from a NetCDF SCRIP file
464 {
465 opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;"
466 "PARALLEL_RESOLVE_SHARED_ENTS;NO_EDGES;NO_MIXED_ELEMENTS;VARIABLE=;";
467 }
468 }
469 }
470
471 #ifdef MOAB_HAVE_MPI
472 int line_size = opts.size();
473 MPI_Bcast( &line_size, 1, MPI_INT, 0, MPI_COMM_WORLD );
474 if( ctx.proc_id != 0 ) opts.resize( line_size );
475 MPI_Bcast( const_cast< char* >( opts.data() ), line_size, MPI_CHAR, 0, MPI_COMM_WORLD );
476 #endif
477 }
478 else
479 return "PARALLEL=BCAST_DELETE;PARTITION=TRIVIAL;PARALLEL_RESOLVE_SHARED_ENTS;";
480 }
481 return opts;
482 }
483 //#define MOAB_DBG
484 int main( int argc, char* argv[] )
485 {
486 moab::ErrorCode rval;
487 NcError error( NcError::verbose_nonfatal );
488 std::stringstream sstr;
489 std::string historyStr;
490
491 int proc_id = 0, nprocs = 1;
492 #ifdef MOAB_HAVE_MPI
493 MPI_Init( &argc, &argv );
494 MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
495 MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
496 #endif
497
498 moab::Core* mbCore = new( std::nothrow ) moab::Core;
499
500 if( nullptr == mbCore )
501 {
502 return 1;
503 }
504
505 // Build the history string
506 for( int ia = 0; ia < argc; ++ia )
507 historyStr += std::string( argv[ia] ) + " ";
508
509 ToolContext* runCtx;
510 #ifdef MOAB_HAVE_MPI
511 moab::ParallelComm* pcomm = new moab::ParallelComm( mbCore, MPI_COMM_WORLD, 0 );
512
513 runCtx = new ToolContext( mbCore, pcomm );
514 const char* writeOptions = ( nprocs > 1 ? "PARALLEL=WRITE_PART" : "" );
515 #else
516 runCtx = new ToolContext( mbCore );
517 const char* writeOptions = "";
518 #endif
519 runCtx->ParseCLOptions( argc, argv );
520
521 const double radius_src = 1.0 /*2.0*acos(-1.0)*/;
522 const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
523
524 moab::DebugOutput& outputFormatter = runCtx->outputFormatter;
525
526 #ifdef MOAB_HAVE_MPI
527 moab::TempestRemapper remapper( mbCore, pcomm );
528 #else
529 moab::TempestRemapper remapper( mbCore );
530 #endif
531 remapper.meshValidate = true;
532 remapper.constructEdgeMap = true;
533 remapper.initialize();
534
535 // Default area_method = lHuiller; Options: Girard, lHuiller, GaussQuadrature (if TR is available)
536 moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::GaussQuadrature );
537
538 Mesh* tempest_mesh = new Mesh();
539 rval = CreateTempestMesh( *runCtx, remapper, tempest_mesh );MB_CHK_ERR( rval );
540
541 if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MEMORY )
542 {
543 // Compute intersections with MOAB
544 // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
545 assert( runCtx->meshes.size() == 3 );
546
547 #ifdef MOAB_HAVE_MPI
548 rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
549 #endif
550
551 // Load the meshes and validate
552 rval = remapper.ConvertTempestMesh( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
553 rval = remapper.ConvertTempestMesh( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
554 rval = remapper.ConvertTempestMesh( moab::Remapper::OverlapMesh );MB_CHK_ERR( rval );
555 if( !runCtx->skip_io )
556 {
557 rval = mbCore->write_mesh( "tempest_intersection.h5m", &runCtx->meshsets[2], 1 );MB_CHK_ERR( rval );
558 }
559
560 // print verbosely about the problem setting
561 size_t velist[6], gvelist[6];
562 {
563 moab::Range rintxverts, rintxelems;
564 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval );
565 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval );
566 velist[0] = rintxverts.size();
567 velist[1] = rintxelems.size();
568
569 moab::Range bintxverts, bintxelems;
570 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval );
571 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval );
572 velist[2] = bintxverts.size();
573 velist[3] = bintxelems.size();
574 }
575
576 moab::EntityHandle intxset; // == remapper.GetMeshSet(moab::Remapper::OverlapMesh);
577
578 // Compute intersections with MOAB
579 {
580 // Create the intersection on the sphere object
581 runCtx->timer_push( "setup the intersector" );
582
583 moab::Intx2MeshOnSphere* mbintx = new moab::Intx2MeshOnSphere( mbCore );
584 mbintx->set_error_tolerance( runCtx->epsrel );
585 mbintx->set_box_error( runCtx->boxeps );
586 mbintx->set_radius_source_mesh( radius_src );
587 mbintx->set_radius_destination_mesh( radius_dest );
588 #ifdef MOAB_HAVE_MPI
589 mbintx->set_parallel_comm( pcomm );
590 #endif
591 rval = mbintx->FindMaxEdges( runCtx->meshsets[0], runCtx->meshsets[1] );MB_CHK_ERR( rval );
592
593 #ifdef MOAB_HAVE_MPI
594 moab::Range local_verts;
595 rval = mbintx->build_processor_euler_boxes( runCtx->meshsets[1], local_verts );MB_CHK_ERR( rval );
596
597 runCtx->timer_pop();
598
599 moab::EntityHandle covering_set;
600 runCtx->timer_push( "communicate the mesh" );
601 // TODO:: needs clarification
602 // we compute just intersection here, no need for extra ghost layers anyway
603 // ghost layers are needed in coverage for bilinear map, which does not actually need intersection
604 // this will be fixed in the future, bilinear map needs just coverage, not intersection
605 // so I am not passing the ghost layer here, even though there is an option in runCtx for a ghost layer
606 rval = mbintx->construct_covering_set( runCtx->meshsets[0], covering_set );MB_CHK_ERR( rval ); // lots of communication if mesh is distributed very differently
607 runCtx->timer_pop();
608
609 // print verbosely about the problem setting
610 {
611 moab::Range cintxverts, cintxelems;
612 rval = mbCore->get_entities_by_dimension( covering_set, 0, cintxverts );MB_CHK_ERR( rval );
613 rval = mbCore->get_entities_by_dimension( covering_set, 2, cintxelems );MB_CHK_ERR( rval );
614 velist[4] = cintxverts.size();
615 velist[5] = cintxelems.size();
616 }
617
618 MPI_Reduce( velist, gvelist, 6, MPI_UINT64_T, MPI_SUM, 0, MPI_COMM_WORLD );
619
620 #else
621 moab::EntityHandle covering_set = runCtx->meshsets[0];
622 for( int i = 0; i < 6; i++ )
623 gvelist[i] = velist[i];
624 #endif
625
626 if( !proc_id )
627 {
628 outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", gvelist[0],
629 gvelist[0] );
630 outputFormatter.printf( 0, "The covering set contains %lu vertices and %lu elements \n", gvelist[2],
631 gvelist[2] );
632 outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", gvelist[1],
633 gvelist[1] );
634 }
635
636 // Now let's invoke the MOAB intersection algorithm in parallel with a
637 // source and target mesh set representing two different decompositions
638 runCtx->timer_push( "compute intersections with MOAB" );
639 rval = mbCore->create_meshset( moab::MESHSET_SET, intxset );MB_CHK_SET_ERR( rval, "Can't create new set" );
640 rval = mbintx->intersect_meshes( covering_set, runCtx->meshsets[1], intxset );MB_CHK_SET_ERR( rval, "Can't compute the intersection of meshes on the sphere" );
641 runCtx->timer_pop();
642
643 // free the memory
644 delete mbintx;
645 }
646
647 {
648 moab::Range intxelems, intxverts;
649 rval = mbCore->get_entities_by_dimension( intxset, 2, intxelems );MB_CHK_ERR( rval );
650 rval = mbCore->get_entities_by_dimension( intxset, 0, intxverts, true );MB_CHK_ERR( rval );
651 outputFormatter.printf( 0, "The intersection set contains %lu elements and %lu vertices \n",
652 intxelems.size(), intxverts.size() );
653
654 moab::IntxAreaUtils areaAdaptorHuiller( moab::IntxAreaUtils::lHuiller ); // lHuiller, GaussQuadrature
655 double initial_sarea =
656 areaAdaptorHuiller.area_on_sphere( mbCore, runCtx->meshsets[0],
657 radius_src ); // use the target to compute the initial area
658 double initial_tarea =
659 areaAdaptorHuiller.area_on_sphere( mbCore, runCtx->meshsets[1],
660 radius_dest ); // use the target to compute the initial area
661 double intx_area = areaAdaptorHuiller.area_on_sphere( mbCore, intxset, radius_src );
662
663 outputFormatter.printf( 0, "mesh areas: source = %12.10f, target = %12.10f, intersection = %12.10f \n",
664 initial_sarea, initial_tarea, intx_area );
665 outputFormatter.printf( 0, "relative error w.r.t source = %12.10e, target = %12.10e \n",
666 fabs( intx_area - initial_sarea ) / initial_sarea,
667 fabs( intx_area - initial_tarea ) / initial_tarea );
668 }
669
670 // Write out our computed intersection file
671 if( !runCtx->skip_io )
672 {
673 rval = mbCore->write_mesh( "moab_intersection.h5m", &intxset, 1 );MB_CHK_ERR( rval );
674 }
675
676 if( runCtx->computeWeights )
677 {
678 runCtx->timer_push( "compute weights with the Tempest meshes" );
679 // Call to generate an offline map with the tempest meshes
680 OfflineMap weightMap;
681 int err = GenerateOfflineMapWithMeshes( *runCtx->meshes[0], *runCtx->meshes[1], *runCtx->meshes[2],
682 runCtx->disc_methods[0], // std::string strInputType
683 runCtx->disc_methods[1], // std::string strOutputType,
684 runCtx->mapOptions, weightMap );
685 runCtx->timer_pop();
686
687 std::map< std::string, std::string > mapAttributes;
688 if( err )
689 {
690 rval = moab::MB_FAILURE;
691 }
692 else
693 {
694 if( !runCtx->skip_io ) weightMap.Write( "outWeights.nc", mapAttributes );
695 }
696 }
697 }
698 else if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MOAB )
699 {
700 // Usage: mpiexec -n 2 tools/mbtempest -t 5 -l mycs_2.h5m -l myico_2.h5m -f myoverlap_2.h5m
701 #ifdef MOAB_HAVE_MPI
702 rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
703 #endif
704
705 // print verbosely about the problem setting
706 size_t velist[4] = {}, gvelist[4] = {};
707 {
708 moab::Range srcverts, srcelems;
709 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, srcverts );MB_CHK_ERR( rval );
710 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, srcelems );MB_CHK_ERR( rval );
711 rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[0] );MB_CHK_ERR( rval );
712 if( runCtx->enforceConvexity )
713 {
714 rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[0], proc_id );MB_CHK_ERR( rval );
715 }
716 rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[0], radius_src );MB_CHK_ERR( rval );
717 // if( !proc_id )
718 // outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", srcverts.size(),
719 // srcelems.size() );
720 velist[0] = srcverts.size();
721 velist[1] = srcelems.size();
722
723 moab::Range tgtverts, tgtelems;
724 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, tgtverts );MB_CHK_ERR( rval );
725 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtelems );MB_CHK_ERR( rval );
726 rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[1] );MB_CHK_ERR( rval );
727 if( runCtx->enforceConvexity )
728 {
729 rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[1], proc_id );MB_CHK_ERR( rval );
730 }
731 rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[1], radius_dest );MB_CHK_ERR( rval );
732 // if( !proc_id )
733 // outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", tgtverts.size(),
734 // tgtelems.size() );
735 velist[2] = tgtverts.size();
736 velist[3] = tgtelems.size();
737 }
738 //rval = mbCore->write_file( "source_mesh.h5m", nullptr, writeOptions, &runCtx->meshsets[0], 1 );MB_CHK_ERR( rval );
739 //rval = mbCore->write_file( "target_mesh.h5m", nullptr, writeOptions, &runCtx->meshsets[1], 1 );MB_CHK_ERR( rval );
740
741 // if( runCtx->nlayers && nprocs > 1 )
742 // {
743 // remapper.ResetMeshSet( moab::Remapper::SourceMesh, runCtx->meshsets[3] );
744 // runCtx->meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh ); // ?
745 // }
746
747 // First compute the covering set such that the target elements are fully covered by the
748 // local source grid
749 runCtx->timer_push( "construct covering set for intersection" );
750 rval = remapper.ConstructCoveringSet( runCtx->epsrel, 1.0, 1.0, runCtx->boxeps, runCtx->rrmGrids,
751 runCtx->useGnomonicProjection, runCtx->nlayers );MB_CHK_ERR( rval );
752 runCtx->timer_pop();
753
754 #ifdef MOAB_HAVE_MPI
755 MPI_Reduce( velist, gvelist, 4, MPI_UINT64_T, MPI_SUM, 0, MPI_COMM_WORLD );
756 #else
757 for( int i = 0; i < 4; i++ )
758 gvelist[i] = velist[i];
759 #endif
760
761 if( !proc_id && runCtx->print_diagnostics )
762 {
763 outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", gvelist[0],
764 gvelist[1] );
765 outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", gvelist[2],
766 gvelist[3] );
767 }
768
769 // Compute intersections with MOAB with either the Kd-tree or the advancing front algorithm
770 runCtx->timer_push( "setup and compute mesh intersections" );
771 rval = remapper.ComputeOverlapMesh( runCtx->kdtreeSearch, false );MB_CHK_ERR( rval );
772 runCtx->timer_pop();
773
774 // print some diagnostic checks to see if the overlap grid resolved the input meshes
775 // correctly
776 double dTotalOverlapArea = 0.0;
777 if( runCtx->print_diagnostics )
778 {
779 moab::IntxAreaUtils areaAdaptorHuiller(
780 moab::IntxAreaUtils::GaussQuadrature ); // lHuiller, GaussQuadrature
781 double local_areas[3],
782 global_areas[3]; // Array for Initial area, and through Method 1 and Method 2
783 // local_areas[0] = area_on_sphere_lHuiller ( mbCore, runCtx->meshsets[1], radius_src );
784 local_areas[0] = areaAdaptorHuiller.area_on_sphere( mbCore, runCtx->meshsets[0], radius_src );
785 local_areas[1] = areaAdaptorHuiller.area_on_sphere( mbCore, runCtx->meshsets[1], radius_dest );
786 local_areas[2] = areaAdaptorHuiller.area_on_sphere( mbCore, runCtx->meshsets[2], radius_src );
787
788 #ifdef MOAB_HAVE_MPI
789 MPI_Allreduce( &local_areas[0], &global_areas[0], 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
790 #else
791 global_areas[0] = local_areas[0];
792 global_areas[1] = local_areas[1];
793 global_areas[2] = local_areas[2];
794 #endif
795 if( !proc_id )
796 {
797 outputFormatter.printf( 0,
798 "initial area: source mesh = %12.14f, target mesh = "
799 "%12.14f, overlap mesh = %12.14f\n",
800 global_areas[0], global_areas[1], global_areas[2] );
801 // outputFormatter.printf ( 0, " area with l'Huiller: %12.14f with Girard:
802 // %12.14f\n", global_areas[2], global_areas[3] ); outputFormatter.printf ( 0, "
803 // relative difference areas = %12.10e\n", fabs ( global_areas[2] - global_areas[3]
804 // ) / global_areas[2] );
805 outputFormatter.printf( 0, "relative error w.r.t source = %12.14e, and target = %12.14e\n",
806 fabs( global_areas[0] - global_areas[2] ) / global_areas[0],
807 fabs( global_areas[1] - global_areas[2] ) / global_areas[1] );
808 }
809 dTotalOverlapArea = global_areas[2];
810 }
811
812 if( runCtx->intxFilename.size() )
813 {
814 moab::EntityHandle writableOverlapSet;
815 rval = mbCore->create_meshset( moab::MESHSET_SET, writableOverlapSet );MB_CHK_SET_ERR( rval, "Can't create new set" );
816 moab::EntityHandle meshOverlapSet = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
817 moab::Range ovEnts;
818 rval = mbCore->get_entities_by_dimension( meshOverlapSet, 2, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
819 rval = mbCore->get_entities_by_dimension( meshOverlapSet, 0, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
820
821 #ifdef MOAB_HAVE_MPI
822 // Do not remove ghosted entities if we still haven't computed weights
823 // Remove ghosted entities from overlap set before writing the new mesh set to file
824 if( nprocs > 1 )
825 {
826 moab::Range ghostedEnts;
827 rval = remapper.GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
828 ovEnts = moab::subtract( ovEnts, ghostedEnts );
829 #ifdef MOAB_DBG
830 if( !runCtx->skip_io )
831 {
832 std::stringstream filename;
833 filename << "aug_overlap" << runCtx->pcomm->rank() << ".h5m";
834 rval = runCtx->mbcore->write_file( filename.str().c_str(), 0, 0, &meshOverlapSet, 1 );MB_CHK_ERR( rval );
835 }
836 #endif
837 }
838 #endif
839 rval = mbCore->add_entities( writableOverlapSet, ovEnts );MB_CHK_SET_ERR( rval, "adding local intx cells failed" );
840
841 #ifdef MOAB_HAVE_MPI
842 #ifdef MOAB_DBG
843 if( nprocs > 1 && !runCtx->skip_io )
844 {
845 std::stringstream filename;
846 filename << "writable_intx_" << runCtx->pcomm->rank() << ".h5m";
847 rval = runCtx->mbcore->write_file( filename.str().c_str(), 0, 0, &writableOverlapSet, 1 );MB_CHK_ERR( rval );
848 }
849 #endif
850 #endif
851
852 size_t lastindex = runCtx->intxFilename.find_last_of( "." );
853 sstr.str( "" );
854 sstr << runCtx->intxFilename.substr( 0, lastindex ) << ".h5m";
855 if( !runCtx->proc_id )
856 std::cout << "Writing out the MOAB intersection mesh file to " << sstr.str() << std::endl;
857
858 // Write out our computed intersection file
859 if( !runCtx->skip_io )
860 {
861 rval = mbCore->write_file( sstr.str().c_str(), nullptr, writeOptions, &writableOverlapSet, 1 );MB_CHK_ERR( rval );
862 }
863 }
864
865 if( runCtx->computeWeights )
866 {
867 runCtx->meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
868 if( !runCtx->proc_id ) std::cout << std::endl;
869
870 runCtx->timer_push( "setup computation of weights" );
871 // Call to generate the remapping weights with the tempest meshes
872 moab::TempestOnlineMap* weightMap = new moab::TempestOnlineMap( &remapper );
873 runCtx->timer_pop();
874
875 runCtx->timer_push( "compute weights with TempestRemap" );
876 rval = weightMap->GenerateRemappingWeights(
877 runCtx->disc_methods[0], // std::string strInputType
878 runCtx->disc_methods[1], // std::string strOutputType,
879 runCtx->mapOptions, // const GenerateOfflineMapAlgorithmOptions& options
880 runCtx->doftag_names[0], // const std::string& source_tag_name
881 runCtx->doftag_names[1] // const std::string& target_tag_name
882 );MB_CHK_ERR( rval );
883 runCtx->timer_pop();
884
885 weightMap->PrintMapStatistics();
886
887 // Invoke the CheckMap routine on the TempestRemap serial interface directly, if running
888 // on a single process
889 if( runCtx->fCheck )
890 {
891 const double dNormalTolerance = 1.0E-8;
892 const double dStrictTolerance = 1.0E-12;
893 weightMap->CheckMap( runCtx->fCheck, runCtx->fCheck, runCtx->fCheck && ( runCtx->ensureMonotonicity ),
894 dNormalTolerance, dStrictTolerance, dTotalOverlapArea );
895 }
896
897 if( runCtx->outFilename.size() )
898 {
899 std::map< std::string, std::string > attrMap;
900 attrMap["MOABversion"] = std::string( MOAB_VERSION );
901 attrMap["Title"] = "MOAB-TempestRemap (mbtempest) Offline Regridding Weight Generator";
902 attrMap["normalization"] = "ovarea";
903 attrMap["remap_options"] = runCtx->mapOptions.strMethod;
904 attrMap["domain_a"] = runCtx->inFilenames[0];
905 attrMap["domain_b"] = runCtx->inFilenames[1];
906 if( runCtx->intxFilename.size() ) attrMap["domain_aUb"] = runCtx->intxFilename;
907 attrMap["map_aPb"] = runCtx->outFilename;
908 attrMap["methodorder_a"] = runCtx->disc_methods[0] + ":" + std::to_string( runCtx->disc_orders[0] ) +
909 ":" + std::string( runCtx->doftag_names[0] );
910 attrMap["concave_a"] = runCtx->mapOptions.fSourceConcave ? "true" : "false";
911 attrMap["methodorder_b"] = runCtx->disc_methods[1] + ":" + std::to_string( runCtx->disc_orders[1] ) +
912 ":" + std::string( runCtx->doftag_names[1] );
913 attrMap["concave_b"] = runCtx->mapOptions.fTargetConcave ? "true" : "false";
914 attrMap["bubble"] = runCtx->mapOptions.fNoBubble ? "false" : "true";
915 attrMap["history"] = historyStr;
916
917 // Write the map file to disk in parallel using either HDF5 or SCRIP interface
918 // in extra case; maybe need a better solution, just create it with the right meshset
919 // from the beginning;
920 if( !runCtx->skip_io )
921 {
922 rval = weightMap->WriteParallelMap( runCtx->outFilename.c_str(), attrMap );MB_CHK_ERR( rval );
923 }
924 }
925
926 if( runCtx->verifyWeights )
927 {
928 // Let us pick a sampling test function for solution evaluation
929 // SH, SV, FH, USERVAR
930 bool userVariable = false;
931 moab::TempestOnlineMap::sample_function testFunction;
932 if( !runCtx->variableToVerify.compare( "SH" ) )
933 testFunction = &sample_slow_harmonic;
934 else if( !runCtx->variableToVerify.compare( "FH" ) )
935 testFunction = &sample_fast_harmonic;
936 else if( !runCtx->variableToVerify.compare( "SV" ) )
937 testFunction = &sample_stationary_vortex;
938 else
939 {
940 userVariable = runCtx->variableToVerify.size() ? true : false;
941 testFunction = runCtx->variableToVerify.size() ? nullptr : sample_stationary_vortex;
942 }
943
944 moab::Tag srcAnalyticalFunction;
945 moab::Tag tgtAnalyticalFunction;
946 moab::Tag tgtProjectedFunction;
947 if( testFunction )
948 {
949 runCtx->timer_push( "describe a solution on source grid" );
950 rval = weightMap->DefineAnalyticalSolution( srcAnalyticalFunction, "AnalyticalSolnSrcExact",
951 moab::Remapper::SourceMesh, testFunction );MB_CHK_ERR( rval );
952 runCtx->timer_pop();
953
954 runCtx->timer_push( "describe a solution on target grid" );
955
956 rval = weightMap->DefineAnalyticalSolution( tgtAnalyticalFunction, "AnalyticalSolnTgtExact",
957 moab::Remapper::TargetMesh, testFunction,
958 &tgtProjectedFunction, "ProjectedSolnTgt" );MB_CHK_ERR( rval );
959 // rval = mbCore->write_file ( "tgtWithSolnTag.h5m", nullptr, writeOptions,
960 // &runCtx->meshsets[1], 1 ); MB_CHK_ERR ( rval );
961 runCtx->timer_pop();
962 }
963 else
964 {
965 rval = mbCore->tag_get_handle( runCtx->variableToVerify.c_str(), srcAnalyticalFunction );MB_CHK_ERR( rval );
966
967 rval = mbCore->tag_get_handle( "ProjectedSolnTgt", 1, moab::MB_TYPE_DOUBLE, tgtProjectedFunction,
968 moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );MB_CHK_ERR( rval );
969 }
970
971 if( !runCtx->skip_io )
972 {
973 rval = mbCore->write_file( "srcWithSolnTag.h5m", nullptr, writeOptions, &runCtx->meshsets[0], 1 );MB_CHK_ERR( rval );
974 }
975
976 runCtx->timer_push( "compute solution projection on target grid" );
977 rval = weightMap->ApplyWeights( srcAnalyticalFunction, tgtProjectedFunction, false, runCtx->cassType );MB_CHK_ERR( rval );
978 runCtx->timer_pop();
979
980 if( !runCtx->skip_io )
981 {
982 rval = mbCore->write_file( "tgtWithSolnTag2.h5m", nullptr, writeOptions, &runCtx->meshsets[1], 1 );MB_CHK_ERR( rval );
983 }
984
985 if( nprocs == 1 && runCtx->baselineFile.size() )
986 {
987 // save the field from tgtWithSolnTag2 in a text file, and global ids for cells
988 moab::Range tgtCells;
989 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtCells );MB_CHK_ERR( rval );
990 std::vector< int > globIds;
991 globIds.resize( tgtCells.size() );
992 std::vector< double > vals;
993 vals.resize( tgtCells.size() );
994 moab::Tag projTag;
995 rval = mbCore->tag_get_handle( "ProjectedSolnTgt", projTag );MB_CHK_ERR( rval );
996 moab::Tag gid = mbCore->globalId_tag();
997 rval = mbCore->tag_get_data( gid, tgtCells, &globIds[0] );MB_CHK_ERR( rval );
998 rval = mbCore->tag_get_data( projTag, tgtCells, &vals[0] );MB_CHK_ERR( rval );
999 std::fstream fs;
1000 fs.open( runCtx->baselineFile.c_str(), std::fstream::out );
1001 fs << std::setprecision( 15 ); // maximum precision for doubles
1002 for( size_t i = 0; i < tgtCells.size(); i++ )
1003 fs << globIds[i] << " " << vals[i] << "\n";
1004 fs.close();
1005 // for good measure, save the source file too, with the tag AnalyticalSolnSrcExact
1006 // it will be used later to test, along with a target file
1007 if( !runCtx->skip_io )
1008 {
1009 rval =
1010 mbCore->write_file( "srcWithSolnTag.h5m", nullptr, writeOptions, &runCtx->meshsets[0], 1 );MB_CHK_ERR( rval );
1011 }
1012 }
1013
1014 // compute error metrics if it is a known analytical functional
1015 if( !userVariable )
1016 {
1017 runCtx->timer_push( "compute error metrics against analytical solution on target grid" );
1018 std::map< std::string, double > errMetrics;
1019 rval = weightMap->ComputeMetrics( moab::Remapper::TargetMesh, tgtAnalyticalFunction,
1020 tgtProjectedFunction, errMetrics, true );MB_CHK_ERR( rval );
1021 runCtx->timer_pop();
1022 }
1023 }
1024
1025 delete weightMap;
1026 }
1027 }
1028
1029 // Clean up
1030 remapper.clear();
1031 delete runCtx;
1032 delete mbCore;
1033
1034 #ifdef MOAB_HAVE_MPI
1035 MPI_Finalize();
1036 #endif
1037 exit( 0 );
1038 }
1039
1040 static moab::ErrorCode CreateTempestMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
1041 {
1042 moab::ErrorCode rval = moab::MB_SUCCESS;
1043 int err;
1044 moab::DebugOutput& outputFormatter = ctx.outputFormatter;
1045
1046 if( !ctx.proc_id )
1047 {
1048 outputFormatter.printf( 0, "Creating TempestRemap Mesh object ...\n" );
1049 }
1050
1051 if( ctx.meshType == moab::TempestRemapper::OVERLAP_FILES )
1052 {
1053 ctx.timer_push( "create Tempest OverlapMesh" );
1054 // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
1055 err = GenerateOverlapMesh( ctx.inFilenames[0], ctx.inFilenames[1], *tempest_mesh, ctx.outFilename, "NetCDF4",
1056 "exact", true );
1057 if( err )
1058 rval = moab::MB_FAILURE;
1059 else
1060 ctx.meshes.push_back( tempest_mesh );
1061 ctx.timer_pop();
1062 }
1063 else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MEMORY )
1064 {
1065 // Load the meshes and validate
1066 ctx.meshsets.resize( 3 );
1067 ctx.meshes.resize( 3 );
1068 ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
1069 ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
1070 ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
1071
1072 ctx.timer_push( "load MOAB Source mesh" );
1073 // First the source
1074 rval = remapper.LoadMesh( moab::Remapper::SourceMesh, ctx.inFilenames[0], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
1075 ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
1076 ctx.timer_pop();
1077
1078 ctx.timer_push( "load MOAB Target mesh" );
1079 // Next the target
1080 rval = remapper.LoadMesh( moab::Remapper::TargetMesh, ctx.inFilenames[1], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
1081 ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
1082 ctx.timer_pop();
1083
1084 ctx.timer_push( "generate TempestRemap OverlapMesh" );
1085 // Now let us construct the overlap mesh, by calling TempestRemap interface directly
1086 // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
1087 err = GenerateOverlapWithMeshes( *ctx.meshes[0], *ctx.meshes[1], *tempest_mesh, "" /*ctx.outFilename*/,
1088 "NetCDF4", "exact", false );
1089 ctx.timer_pop();
1090 if( err )
1091 rval = moab::MB_FAILURE;
1092 else
1093 {
1094 remapper.SetMesh( moab::Remapper::OverlapMesh, tempest_mesh );
1095 ctx.meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
1096 // ctx.meshes.push_back(*tempest_mesh);
1097 }
1098 }
1099 else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MOAB )
1100 {
1101 ctx.meshsets.resize( 3 );
1102 ctx.meshes.resize( 3 );
1103 ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
1104 ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
1105 ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
1106
1107 const double radius_src = 1.0 /*2.0*acos(-1.0)*/;
1108 const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
1109
1110 std::vector< int > smetadata, tmetadata;
1111 //const char* additional_read_opts = ( ctx.n_procs > 1 ? "NO_SET_CONTAINING_PARENTS;" : "" );
1112 std::string additional_read_opts_src = get_file_read_options( ctx, ctx.inFilenames[0] );
1113
1114 ctx.timer_push( "load MOAB Source mesh" );
1115 // Load the source mesh and validate
1116 rval =
1117 remapper.LoadNativeMesh( ctx.inFilenames[0], ctx.meshsets[0], smetadata, additional_read_opts_src.c_str() );MB_CHK_ERR( rval );
1118 if( smetadata.size() ) remapper.SetMeshType( moab::Remapper::SourceMesh, smetadata );
1119 ctx.timer_pop();
1120
1121 ctx.timer_push( "preprocess MOAB Source mesh" );
1122 // Rescale the radius of both to compute the intersection
1123 rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[0], radius_src );MB_CHK_ERR( rval );
1124
1125 ctx.timer_pop();
1126
1127 // Load the target mesh and validate
1128 std::string addititional_read_opts_tgt = get_file_read_options( ctx, ctx.inFilenames[1] );
1129 // addititional_read_opts_tgt += "PARALLEL_GHOSTS=2.0.3;PARALLEL_THIN_GHOST_LAYER;";
1130
1131 ctx.timer_push( "load MOAB Target mesh" );
1132 rval = remapper.LoadNativeMesh( ctx.inFilenames[1], ctx.meshsets[1], tmetadata,
1133 addititional_read_opts_tgt.c_str() );MB_CHK_ERR( rval );
1134 if( tmetadata.size() ) remapper.SetMeshType( moab::Remapper::TargetMesh, tmetadata );
1135 ctx.timer_pop();
1136
1137 ctx.timer_push( "preprocess MOAB Target mesh" );
1138 rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[1], radius_dest );MB_CHK_ERR( rval );
1139 ctx.timer_pop();
1140
1141 if( ctx.computeWeights )
1142 {
1143 ctx.timer_push( "convert MOAB meshes to TempestRemap meshes in memory" );
1144 // convert MOAB representation to TempestRemap's Mesh for source
1145 rval = remapper.ConvertMeshToTempest( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
1146 ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
1147
1148 // convert MOAB representation to TempestRemap's Mesh for target
1149 rval = remapper.ConvertMeshToTempest( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
1150 ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
1151 ctx.timer_pop();
1152 }
1153 }
1154 else if( ctx.meshType == moab::TempestRemapper::ICO )
1155 {
1156 ctx.timer_push( "generate ICO mesh with TempestRemap" );
1157 err = GenerateICOMesh( *tempest_mesh, ctx.blockSize, ctx.computeDual, ctx.outFilename, "NetCDF4" );
1158 ctx.timer_pop();
1159
1160 if( err )
1161 rval = moab::MB_FAILURE;
1162 else
1163 ctx.meshes.push_back( tempest_mesh );
1164 }
1165 else if( ctx.meshType == moab::TempestRemapper::RLL )
1166 {
1167 ctx.timer_push( "generate RLL mesh with TempestRemap" );
1168 err = GenerateRLLMesh( *tempest_mesh, // Mesh& meshOut,
1169 ctx.blockSize * 2, ctx.blockSize, // int nLongitudes, int nLatitudes,
1170 0.0, 360.0, // double dLonBegin, double dLonEnd,
1171 -90.0, 90.0, // double dLatBegin, double dLatEnd,
1172 false, false, false, // bool fGlobalCap, bool fFlipLatLon, bool fForceGlobal,
1173 "" /*ctx.inFilename*/, "",
1174 "", // std::string strInputFile, std::string strInputFileLonName, std::string
1175 // strInputFileLatName,
1176 ctx.outFilename, "NetCDF4", // std::string strOutputFile, std::string strOutputFormat
1177 true // bool fVerbose
1178 );
1179 ctx.timer_pop();
1180
1181 if( err )
1182 rval = moab::MB_FAILURE;
1183 else
1184 ctx.meshes.push_back( tempest_mesh );
1185 }
1186 else // default
1187 {
1188 ctx.timer_push( "generate CS mesh with TempestRemap" );
1189 err = GenerateCSMesh( *tempest_mesh, ctx.blockSize, ctx.outFilename, "NetCDF4" );
1190 ctx.timer_pop();
1191 if( err )
1192 rval = moab::MB_FAILURE;
1193 else
1194 ctx.meshes.push_back( tempest_mesh );
1195 }
1196
1197 if( ctx.meshType != moab::TempestRemapper::OVERLAP_MOAB && !tempest_mesh )
1198 {
1199 std::cout << "Tempest Mesh is not a complete object; Quitting...";
1200 exit( -1 );
1201 }
1202
1203 return rval;
1204 }
1205
1206 #undef MOAB_DBG
1207 ///////////////////////////////////////////////
1208 // Test functions
1209
1210 double sample_slow_harmonic( double dLon, double dLat )
1211 {
1212 return ( 2.0 + cos( dLat ) * cos( dLat ) * cos( 2.0 * dLon ) );
1213 }
1214
1215 double sample_fast_harmonic( double dLon, double dLat )
1216 {
1217 return ( 2.0 + pow( sin( 2.0 * dLat ), 16.0 ) * cos( 16.0 * dLon ) );
1218 // return (2.0 + pow(cos(2.0 * dLat), 16.0) * cos(16.0 * dLon));
1219 }
1220
1221 double sample_constant( double /*dLon*/, double /*dLat*/ )
1222 {
1223 return 1.0;
1224 }
1225
1226 double sample_stationary_vortex( double dLon, double dLat )
1227 {
1228 const double dLon0 = 0.0;
1229 const double dLat0 = 0.6;
1230 const double dR0 = 3.0;
1231 const double dD = 5.0;
1232 const double dT = 6.0;
1233
1234 /// Find the rotated longitude and latitude of a point on a sphere
1235 /// with pole at (dLonC, dLatC).
1236 {
1237 double dSinC = sin( dLat0 );
1238 double dCosC = cos( dLat0 );
1239 double dCosT = cos( dLat );
1240 double dSinT = sin( dLat );
1241
1242 double dTrm = dCosT * cos( dLon - dLon0 );
1243 double dX = dSinC * dTrm - dCosC * dSinT;
1244 double dY = dCosT * sin( dLon - dLon0 );
1245 double dZ = dSinC * dSinT + dCosC * dTrm;
1246
1247 dLon = atan2( dY, dX );
1248 if( dLon < 0.0 )
1249 {
1250 dLon += 2.0 * M_PI;
1251 }
1252 dLat = asin( dZ );
1253 }
1254
1255 double dRho = dR0 * cos( dLat );
1256 double dVt = 3.0 * sqrt( 3.0 ) / 2.0 / cosh( dRho ) / cosh( dRho ) * tanh( dRho );
1257
1258 double dOmega;
1259 if( dRho == 0.0 )
1260 {
1261 dOmega = 0.0;
1262 }
1263 else
1264 {
1265 dOmega = dVt / dRho;
1266 }
1267
1268 return ( 1.0 - tanh( dRho / dD * sin( dLon - dOmega * dT ) ) );
1269 }
1270
1271 ///////////////////////////////////////////////
1.9.1.