MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* 00002 * Usage: MOAB-Tempest tool 00003 * 00004 * Generate a Cubed-Sphere mesh: ./mbtempest -t 0 -res 25 -f cubed_sphere_mesh.exo 00005 * Generate a RLL mesh: ./mbtempest -t 1 -res 25 -f rll_mesh.exo 00006 * Generate a Icosahedral-Sphere mesh: ./mbtempest -t 2 -res 25 <-dual> -f icosa_mesh.exo 00007 * 00008 * Now you can compute the intersections between the meshes too! 00009 * 00010 * Generate the overlap mesh: ./mbtempest -t 3 -l cubed_sphere_mesh.exo -l rll_mesh.exo -f 00011 * overlap_mesh.exo 00012 * 00013 */ 00014 00015 #include <iostream> 00016 #include <cstdlib> 00017 #include <vector> 00018 #include <string> 00019 #include <sstream> 00020 #include <cassert> 00021 00022 #include "moab/Core.hpp" 00023 #include "moab/IntxMesh/IntxUtils.hpp" 00024 #include "moab/Remapping/TempestRemapper.hpp" 00025 #include "moab/Remapping/TempestOnlineMap.hpp" 00026 #include "moab/ProgOptions.hpp" 00027 #include "moab/CpuTimer.hpp" 00028 #include "DebugOutput.hpp" 00029 00030 //#ifndef MOAB_HAVE_MPI 00031 // #error mbtempest tool requires MPI configuration 00032 //#endif 00033 00034 #ifdef MOAB_HAVE_MPI 00035 // MPI includes 00036 #include "moab_mpi.h" 00037 #include "moab/ParallelComm.hpp" 00038 #include "MBParallelConventions.h" 00039 #endif 00040 00041 struct ToolContext 00042 { 00043 moab::Interface* mbcore; 00044 #ifdef MOAB_HAVE_MPI 00045 moab::ParallelComm* pcomm; 00046 #endif 00047 const int proc_id, n_procs; 00048 moab::DebugOutput outputFormatter; 00049 int blockSize; 00050 std::vector< std::string > inFilenames; 00051 std::vector< Mesh* > meshes; 00052 std::vector< moab::EntityHandle > meshsets; 00053 std::vector< int > disc_orders; 00054 std::vector< std::string > disc_methods; 00055 std::vector< std::string > doftag_names; 00056 std::string outFilename; 00057 std::string intxFilename; 00058 moab::TempestRemapper::TempestMeshType meshType; 00059 bool computeDual; 00060 bool computeWeights; 00061 bool verifyWeights; 00062 bool enforceConvexity; 00063 int ensureMonotonicity; 00064 bool fNoConservation; 00065 bool fVolumetric; 00066 bool rrmGrids; 00067 bool kdtreeSearch; 00068 bool fNoBubble, fInputConcave, fOutputConcave, fNoCheck; 00069 00070 #ifdef MOAB_HAVE_MPI 00071 ToolContext( moab::Interface* icore, moab::ParallelComm* p_pcomm ) 00072 : mbcore( icore ), pcomm( p_pcomm ), proc_id( pcomm->rank() ), n_procs( pcomm->size() ), 00073 outputFormatter( std::cout, pcomm->rank(), 0 ), 00074 #else 00075 ToolContext ( moab::Interface* icore ) : 00076 mbcore(icore), 00077 proc_id ( 0 ), n_procs ( 1 ), 00078 outputFormatter ( std::cout, 0, 0 ), 00079 #endif 00080 blockSize( 5 ), outFilename( "output.exo" ), intxFilename( "" ), meshType( moab::TempestRemapper::DEFAULT ), 00081 computeDual( false ), computeWeights( false ), verifyWeights( false ), enforceConvexity(false), ensureMonotonicity( 0 ), 00082 fNoConservation( false ), fVolumetric( false ), rrmGrids( false ), kdtreeSearch( true ), fNoBubble( false ), 00083 fInputConcave( false ), fOutputConcave( false ), fNoCheck( false ) 00084 { 00085 inFilenames.resize( 2 ); 00086 doftag_names.resize( 2 ); 00087 timer = new moab::CpuTimer(); 00088 00089 outputFormatter.set_prefix( "[MBTempest]: " ); 00090 } 00091 00092 ~ToolContext() 00093 { 00094 // for (unsigned i=0; i < meshes.size(); ++i) delete meshes[i]; 00095 meshes.clear(); 00096 inFilenames.clear(); 00097 disc_orders.clear(); 00098 disc_methods.clear(); 00099 doftag_names.clear(); 00100 outFilename.clear(); 00101 intxFilename.clear(); 00102 meshsets.clear(); 00103 delete timer; 00104 } 00105 00106 void timer_push( std::string operation ) 00107 { 00108 timer_ops = timer->time_since_birth(); 00109 opName = operation; 00110 } 00111 00112 void timer_pop() 00113 { 00114 double locElapsed = timer->time_since_birth() - timer_ops, avgElapsed = 0, maxElapsed = 0; 00115 #ifdef MOAB_HAVE_MPI 00116 MPI_Reduce( &locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm() ); 00117 MPI_Reduce( &locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, pcomm->comm() ); 00118 #else 00119 maxElapsed = locElapsed; 00120 avgElapsed = locElapsed; 00121 #endif 00122 if( !proc_id ) 00123 { 00124 avgElapsed /= n_procs; 00125 std::cout << "[LOG] Time taken to " << opName.c_str() << ": max = " << maxElapsed 00126 << ", avg = " << avgElapsed << "\n"; 00127 } 00128 // std::cout << "\n[LOG" << proc_id << "] Time taken to " << opName << " = " << 00129 // timer->time_since_birth() - timer_ops << std::endl; 00130 opName.clear(); 00131 } 00132 00133 void ParseCLOptions( int argc, char* argv[] ) 00134 { 00135 ProgOptions opts; 00136 int imeshType = 0; 00137 std::string expectedFName = "output.exo"; 00138 std::string expectedMethod = "fv"; 00139 std::string expectedDofTagName = "GLOBAL_ID"; 00140 int expectedOrder = 1; 00141 00142 if( !proc_id ) 00143 { 00144 std::cout << "Command line options provided to mbtempest:\n "; 00145 for( int iarg = 0; iarg < argc; ++iarg ) 00146 std::cout << argv[iarg] << " "; 00147 std::cout << std::endl << std::endl; 00148 } 00149 00150 opts.addOpt< int >( "type,t", 00151 "Type of mesh (default=CS; Choose from [CS=0, RLL=1, ICO=2, OVERLAP_FILES=3, " 00152 "OVERLAP_MEMORY=4, OVERLAP_MOAB=5])", 00153 &imeshType ); 00154 opts.addOpt< int >( "res,r", "Resolution of the mesh (default=5)", &blockSize ); 00155 opts.addOpt< void >( "dual,d", "Output the dual of the mesh (generally relevant only for ICO mesh)", 00156 &computeDual ); 00157 opts.addOpt< void >( "weights,w", 00158 "Compute and output the weights using the overlap mesh (generally " 00159 "relevant only for OVERLAP mesh)", 00160 &computeWeights ); 00161 opts.addOpt< void >( "noconserve,c", 00162 "Do not apply conservation to the resultant weights (relevant only " 00163 "when computing weights)", 00164 &fNoConservation ); 00165 opts.addOpt< void >( "volumetric,v", 00166 "Apply a volumetric projection to compute the weights (relevant only " 00167 "when computing weights)", 00168 &fVolumetric ); 00169 opts.addOpt< void >( "rrmgrids", 00170 "At least one of the meshes is a regionally refined grid (relevant to " 00171 "accelerate intersection computation)", 00172 &rrmGrids ); 00173 opts.addOpt< void >( "nocheck", "Do not check the generated map for conservation and consistency", &fNoCheck ); 00174 opts.addOpt< void >( "advfront,a", 00175 "Use the advancing front intersection instead of the Kd-tree based algorithm " 00176 "to compute mesh intersections" ); 00177 opts.addOpt< void >( "verify", 00178 "Verify the accuracy of the maps by projecting analytical functions " 00179 "from source to target " 00180 "grid by applying the maps", 00181 &verifyWeights ); 00182 opts.addOpt< int >( "monotonic,n", "Ensure monotonicity in the weight generation", &ensureMonotonicity ); 00183 opts.addOpt< std::string >( "load,l", "Input mesh filenames (a source and target mesh)", &expectedFName ); 00184 opts.addOpt< int >( "order,o", "Discretization orders for the source and target solution fields", 00185 &expectedOrder ); 00186 opts.addOpt< std::string >( "method,m", "Discretization method for the source and target solution fields", 00187 &expectedMethod ); 00188 opts.addOpt< std::string >( "global_id,g", 00189 "Tag name that contains the global DoF IDs for source and target solution fields", 00190 &expectedDofTagName ); 00191 opts.addOpt< std::string >( "file,f", "Output remapping weights filename", &outFilename ); 00192 opts.addOpt< std::string >( "intx,i", "Output TempestRemap intersection mesh filename", &intxFilename ); 00193 00194 opts.addOpt< void >( "enforce_convexity", 00195 "check convexity of input meshes to compute mesh intersections" , &enforceConvexity); 00196 00197 opts.parseCommandLine( argc, argv ); 00198 00199 // By default - use Kd-tree based search; if user asks for advancing front, disable Kd-tree 00200 // algorithm 00201 kdtreeSearch = opts.numOptSet( "advfront,a" ) == 0; 00202 00203 switch( imeshType ) 00204 { 00205 case 0: 00206 meshType = moab::TempestRemapper::CS; 00207 break; 00208 00209 case 1: 00210 meshType = moab::TempestRemapper::RLL; 00211 break; 00212 00213 case 2: 00214 meshType = moab::TempestRemapper::ICO; 00215 break; 00216 00217 case 3: 00218 meshType = moab::TempestRemapper::OVERLAP_FILES; 00219 break; 00220 00221 case 4: 00222 meshType = moab::TempestRemapper::OVERLAP_MEMORY; 00223 break; 00224 00225 case 5: 00226 meshType = moab::TempestRemapper::OVERLAP_MOAB; 00227 break; 00228 00229 default: 00230 meshType = moab::TempestRemapper::DEFAULT; 00231 break; 00232 } 00233 00234 if( meshType > moab::TempestRemapper::ICO ) 00235 { 00236 opts.getOptAllArgs( "load,l", inFilenames ); 00237 opts.getOptAllArgs( "order,o", disc_orders ); 00238 opts.getOptAllArgs( "method,m", disc_methods ); 00239 opts.getOptAllArgs( "global_id,i", doftag_names ); 00240 00241 if( disc_orders.size() == 0 ) { disc_orders.resize( 2, 1 ); } 00242 00243 if( disc_orders.size() == 1 ) { disc_orders.push_back( 1 ); } 00244 00245 if( disc_methods.size() == 0 ) { disc_methods.resize( 2, "fv" ); } 00246 00247 if( disc_methods.size() == 1 ) { disc_methods.push_back( "fv" ); } 00248 00249 if( doftag_names.size() == 0 ) { doftag_names.resize( 2, "GLOBAL_ID" ); } 00250 00251 if( doftag_names.size() == 1 ) { doftag_names.push_back( "GLOBAL_ID" ); } 00252 00253 assert( inFilenames.size() == 2 ); 00254 assert( disc_orders.size() == 2 ); 00255 assert( disc_methods.size() == 2 ); 00256 } 00257 00258 expectedFName.clear(); 00259 } 00260 00261 private: 00262 moab::CpuTimer* timer; 00263 double timer_ops; 00264 std::string opName; 00265 }; 00266 00267 // Forward declare some methods 00268 static moab::ErrorCode CreateTempestMesh( ToolContext&, moab::TempestRemapper& remapper, Mesh* ); 00269 inline double sample_slow_harmonic( double dLon, double dLat ); 00270 inline double sample_fast_harmonic( double dLon, double dLat ); 00271 inline double sample_constant( double dLon, double dLat ); 00272 inline double sample_stationary_vortex( double dLon, double dLat ); 00273 00274 int main( int argc, char* argv[] ) 00275 { 00276 moab::ErrorCode rval; 00277 NcError error( NcError::verbose_nonfatal ); 00278 std::stringstream sstr; 00279 00280 int proc_id = 0, nprocs = 1; 00281 #ifdef MOAB_HAVE_MPI 00282 MPI_Init( &argc, &argv ); 00283 MPI_Comm_rank( MPI_COMM_WORLD, &proc_id ); 00284 MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); 00285 #endif 00286 00287 moab::Interface* mbCore = new( std::nothrow ) moab::Core; 00288 00289 if( NULL == mbCore ) { return 1; } 00290 00291 ToolContext* runCtx; 00292 #ifdef MOAB_HAVE_MPI 00293 moab::ParallelComm* pcomm = new moab::ParallelComm( mbCore, MPI_COMM_WORLD, 0 ); 00294 00295 runCtx = new ToolContext( mbCore, pcomm ); 00296 const char* writeOptions = ( nprocs > 1 ? "PARALLEL=WRITE_PART" : "" ); 00297 #else 00298 runCtx = new ToolContext( mbCore ); 00299 const char* writeOptions = ""; 00300 #endif 00301 runCtx->ParseCLOptions( argc, argv ); 00302 00303 const double radius_src = 1.0 /*2.0*acos(-1.0)*/; 00304 const double radius_dest = 1.0 /*2.0*acos(-1.0)*/; 00305 00306 moab::DebugOutput& outputFormatter = runCtx->outputFormatter; 00307 00308 #ifdef MOAB_HAVE_MPI 00309 moab::TempestRemapper remapper( mbCore, pcomm ); 00310 #else 00311 moab::TempestRemapper remapper( mbCore ); 00312 #endif 00313 remapper.meshValidate = true; 00314 remapper.constructEdgeMap = false; 00315 remapper.initialize(); 00316 00317 // Default area_method = lHuiller; Options: Girard, GaussQuadrature (if TR is available) 00318 #ifdef MOAB_HAVE_TEMPESTREMAP 00319 moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::GaussQuadrature ); 00320 #else 00321 moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::lHuiller ); 00322 #endif 00323 00324 Mesh* tempest_mesh = new Mesh(); 00325 runCtx->timer_push( "create Tempest mesh" ); 00326 rval = CreateTempestMesh( *runCtx, remapper, tempest_mesh );MB_CHK_ERR( rval ); 00327 runCtx->timer_pop(); 00328 00329 // Some constant parameters 00330 const double epsrel = 1.e-15; 00331 const double boxeps = 1e-6; 00332 00333 if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MEMORY ) 00334 { 00335 // Compute intersections with MOAB 00336 // For the overlap method, choose between: "fuzzy", "exact" or "mixed" 00337 assert( runCtx->meshes.size() == 3 ); 00338 00339 #ifdef MOAB_HAVE_MPI 00340 rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval ); 00341 #endif 00342 00343 // Load the meshes and validate 00344 rval = remapper.ConvertTempestMesh( moab::Remapper::SourceMesh );MB_CHK_ERR( rval ); 00345 rval = remapper.ConvertTempestMesh( moab::Remapper::TargetMesh );MB_CHK_ERR( rval ); 00346 remapper.SetMeshType( moab::Remapper::OverlapMesh, moab::TempestRemapper::OVERLAP_FILES ); 00347 rval = remapper.ConvertTempestMesh( moab::Remapper::OverlapMesh );MB_CHK_ERR( rval ); 00348 rval = mbCore->write_mesh( "tempest_intersection.h5m", &runCtx->meshsets[2], 1 );MB_CHK_ERR( rval ); 00349 00350 // print verbosely about the problem setting 00351 { 00352 moab::Range rintxverts, rintxelems; 00353 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval ); 00354 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval ); 00355 outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", rintxverts.size(), 00356 rintxelems.size() ); 00357 00358 moab::Range bintxverts, bintxelems; 00359 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval ); 00360 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval ); 00361 outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", bintxverts.size(), 00362 bintxelems.size() ); 00363 } 00364 00365 moab::EntityHandle intxset; // == remapper.GetMeshSet(moab::Remapper::OverlapMesh); 00366 00367 // Compute intersections with MOAB 00368 { 00369 // Create the intersection on the sphere object 00370 runCtx->timer_push( "setup the intersector" ); 00371 00372 moab::Intx2MeshOnSphere* mbintx = new moab::Intx2MeshOnSphere( mbCore ); 00373 mbintx->set_error_tolerance( epsrel ); 00374 mbintx->set_box_error( boxeps ); 00375 mbintx->set_radius_source_mesh( radius_src ); 00376 mbintx->set_radius_destination_mesh( radius_dest ); 00377 #ifdef MOAB_HAVE_MPI 00378 mbintx->set_parallel_comm( pcomm ); 00379 #endif 00380 rval = mbintx->FindMaxEdges( runCtx->meshsets[0], runCtx->meshsets[1] );MB_CHK_ERR( rval ); 00381 00382 #ifdef MOAB_HAVE_MPI 00383 moab::Range local_verts; 00384 rval = mbintx->build_processor_euler_boxes( runCtx->meshsets[1], local_verts );MB_CHK_ERR( rval ); 00385 00386 runCtx->timer_pop(); 00387 00388 moab::EntityHandle covering_set; 00389 runCtx->timer_push( "communicate the mesh" ); 00390 rval = mbintx->construct_covering_set( runCtx->meshsets[0], covering_set );MB_CHK_ERR( rval ); // lots of communication if mesh is distributed very differently 00391 runCtx->timer_pop(); 00392 #else 00393 moab::EntityHandle covering_set = runCtx->meshsets[0]; 00394 #endif 00395 // Now let's invoke the MOAB intersection algorithm in parallel with a 00396 // source and target mesh set representing two different decompositions 00397 runCtx->timer_push( "compute intersections with MOAB" ); 00398 rval = mbCore->create_meshset( moab::MESHSET_SET, intxset );MB_CHK_SET_ERR( rval, "Can't create new set" ); 00399 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" ); 00400 runCtx->timer_pop(); 00401 00402 // free the memory 00403 delete mbintx; 00404 } 00405 00406 { 00407 moab::Range intxelems, intxverts; 00408 rval = mbCore->get_entities_by_dimension( intxset, 2, intxelems );MB_CHK_ERR( rval ); 00409 rval = mbCore->get_entities_by_dimension( intxset, 0, intxverts, true );MB_CHK_ERR( rval ); 00410 outputFormatter.printf( 0, "The intersection set contains %lu elements and %lu vertices \n", 00411 intxelems.size(), intxverts.size() ); 00412 00413 double initial_sarea = 00414 areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0], 00415 radius_src ); // use the target to compute the initial area 00416 double initial_tarea = 00417 areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1], 00418 radius_dest ); // use the target to compute the initial area 00419 double intx_area = areaAdaptor.area_on_sphere( mbCore, intxset, radius_src ); 00420 00421 outputFormatter.printf( 0, "mesh areas: source = %12.10f, target = %12.10f, intersection = %12.10f \n", 00422 initial_sarea, initial_tarea, intx_area ); 00423 outputFormatter.printf( 0, "relative error w.r.t source = %12.10e, target = %12.10e \n", 00424 fabs( intx_area - initial_sarea ) / initial_sarea, 00425 fabs( intx_area - initial_tarea ) / initial_tarea ); 00426 } 00427 00428 // Write out our computed intersection file 00429 rval = mbCore->write_mesh( "moab_intersection.h5m", &intxset, 1 );MB_CHK_ERR( rval ); 00430 00431 if( runCtx->computeWeights ) 00432 { 00433 runCtx->timer_push( "compute weights with the Tempest meshes" ); 00434 // Call to generate an offline map with the tempest meshes 00435 OfflineMap weightMap; 00436 int err = GenerateOfflineMapWithMeshes( 00437 weightMap, *runCtx->meshes[0], *runCtx->meshes[1], *runCtx->meshes[2], "", 00438 "", // std::string strInputMeta, std::string strOutputMeta, 00439 runCtx->disc_methods[0], 00440 runCtx->disc_methods[1], // std::string strInputType, std::string strOutputType, 00441 runCtx->disc_orders[0], runCtx->disc_orders[1], // int nPin=4, int nPout=4, 00442 runCtx->fNoBubble, true, 00443 runCtx->ensureMonotonicity // bool fNoBubble = false, bool fCorrectAreas = false, 00444 // int fMonotoneTypeID = 0 00445 ); 00446 runCtx->timer_pop(); 00447 00448 std::map< std::string, std::string > mapAttributes; 00449 if( err ) { rval = moab::MB_FAILURE; } 00450 else 00451 { 00452 weightMap.Write( "outWeights.nc", mapAttributes ); 00453 } 00454 } 00455 } 00456 else if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MOAB ) 00457 { 00458 // Usage: mpiexec -n 2 tools/mbtempest -t 5 -l mycs_2.h5m -l myico_2.h5m -f myoverlap_2.h5m 00459 #ifdef MOAB_HAVE_MPI 00460 rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval ); 00461 #endif 00462 // print verbosely about the problem setting 00463 { 00464 moab::Range rintxverts, rintxelems; 00465 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval ); 00466 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval ); 00467 rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[0] );MB_CHK_ERR( rval ); 00468 if (runCtx->enforceConvexity) 00469 { 00470 rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[0], proc_id );MB_CHK_ERR( rval ); 00471 } 00472 rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[0], radius_src );MB_CHK_ERR( rval ); 00473 if( !proc_id ) 00474 outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", 00475 rintxverts.size(), rintxelems.size() ); 00476 00477 moab::Range bintxverts, bintxelems; 00478 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval ); 00479 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval ); 00480 rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[1] );MB_CHK_ERR( rval ); 00481 if (runCtx->enforceConvexity) 00482 { 00483 rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[1], proc_id );MB_CHK_ERR( rval ); 00484 } 00485 rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[1], radius_dest );MB_CHK_ERR( rval ); 00486 if( !proc_id ) 00487 outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", 00488 bintxverts.size(), bintxelems.size() ); 00489 } 00490 00491 // First compute the covering set such that the target elements are fully covered by the 00492 // lcoal source grid 00493 runCtx->timer_push( "construct covering set for intersection" ); 00494 rval = remapper.ConstructCoveringSet( epsrel, 1.0, 1.0, boxeps, runCtx->rrmGrids );MB_CHK_ERR( rval ); 00495 runCtx->timer_pop(); 00496 00497 // Compute intersections with MOAB with either the Kd-tree or the advancing front algorithm 00498 runCtx->timer_push( "setup and compute mesh intersections" ); 00499 rval = remapper.ComputeOverlapMesh( runCtx->kdtreeSearch, false );MB_CHK_ERR( rval ); 00500 runCtx->timer_pop(); 00501 00502 // print some diagnostic checks to see if the overlap grid resolved the input meshes 00503 // correctly 00504 { 00505 double local_areas[3], 00506 global_areas[3]; // Array for Initial area, and through Method 1 and Method 2 00507 // local_areas[0] = area_on_sphere_lHuiller ( mbCore, runCtx->meshsets[1], radius_src ); 00508 local_areas[0] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0], radius_src ); 00509 local_areas[1] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1], radius_dest ); 00510 local_areas[2] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[2], radius_src ); 00511 00512 #ifdef MOAB_HAVE_MPI 00513 MPI_Allreduce( &local_areas[0], &global_areas[0], 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); 00514 #else 00515 global_areas[0] = local_areas[0]; 00516 global_areas[1] = local_areas[1]; 00517 global_areas[2] = local_areas[2]; 00518 #endif 00519 if( !proc_id ) 00520 { 00521 outputFormatter.printf( 0, 00522 "initial area: source mesh = %12.14f, target mesh = " 00523 "%12.14f, overlap mesh = %12.14f\n", 00524 global_areas[0], global_areas[1], global_areas[2] ); 00525 // outputFormatter.printf ( 0, " area with l'Huiller: %12.14f with Girard: 00526 // %12.14f\n", global_areas[2], global_areas[3] ); outputFormatter.printf ( 0, " 00527 // relative difference areas = %12.10e\n", fabs ( global_areas[2] - global_areas[3] 00528 // ) / global_areas[2] ); 00529 outputFormatter.printf( 0, "relative error w.r.t source = %12.14e, and target = %12.14e\n", 00530 fabs( global_areas[0] - global_areas[2] ) / global_areas[0], 00531 fabs( global_areas[1] - global_areas[2] ) / global_areas[1] ); 00532 } 00533 } 00534 00535 if( runCtx->intxFilename.size() ) 00536 { 00537 moab::EntityHandle writableOverlapSet; 00538 rval = mbCore->create_meshset( moab::MESHSET_SET, writableOverlapSet );MB_CHK_SET_ERR( rval, "Can't create new set" ); 00539 moab::EntityHandle meshOverlapSet = remapper.GetMeshSet( moab::Remapper::OverlapMesh ); 00540 moab::Range ovEnts; 00541 rval = mbCore->get_entities_by_dimension( meshOverlapSet, 2, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" ); 00542 rval = mbCore->get_entities_by_dimension( meshOverlapSet, 0, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" ); 00543 00544 #ifdef MOAB_HAVE_MPI 00545 // Do not remove ghosted entities if we still haven't computed weights 00546 // Remove ghosted entities from overlap set before writing the new mesh set to file 00547 if( nprocs > 1 ) 00548 { 00549 moab::Range ghostedEnts; 00550 rval = remapper.GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval ); 00551 ovEnts = moab::subtract( ovEnts, ghostedEnts ); 00552 } 00553 #endif 00554 rval = mbCore->add_entities( writableOverlapSet, ovEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" ); 00555 00556 size_t lastindex = runCtx->intxFilename.find_last_of( "." ); 00557 sstr.str( "" ); 00558 sstr << runCtx->intxFilename.substr( 0, lastindex ) << ".h5m"; 00559 if( !runCtx->proc_id ) 00560 std::cout << "Writing out the MOAB intersection mesh file to " << sstr.str() << std::endl; 00561 00562 // Write out our computed intersection file 00563 rval = mbCore->write_file( sstr.str().c_str(), NULL, writeOptions, &writableOverlapSet, 1 );MB_CHK_ERR( rval ); 00564 } 00565 00566 if( runCtx->computeWeights ) 00567 { 00568 runCtx->meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh ); 00569 if( !runCtx->proc_id ) std::cout << std::endl; 00570 00571 runCtx->timer_push( "setup computation of weights" ); 00572 // Call to generate the remapping weights with the tempest meshes 00573 moab::TempestOnlineMap* weightMap = new moab::TempestOnlineMap( &remapper ); 00574 runCtx->timer_pop(); 00575 00576 runCtx->timer_push( "compute weights with TempestRemap" ); 00577 00578 rval = weightMap->GenerateRemappingWeights( 00579 runCtx->disc_methods[0], 00580 runCtx->disc_methods[1], // std::string strInputType, std::string strOutputType, 00581 runCtx->disc_orders[0], runCtx->disc_orders[1], // int nPin=4, int nPout=4, 00582 runCtx->fNoBubble, 00583 runCtx->ensureMonotonicity, // bool fNoBubble=true, int fMonotoneTypeID=0, 00584 runCtx->fVolumetric, runCtx->fNoConservation, 00585 runCtx->fNoCheck, // bool fVolumetric=false, bool fNoConservation=false, bool 00586 // fNoCheck=false, 00587 runCtx->doftag_names[0], 00588 runCtx->doftag_names[1], // std::string source_tag_name, std::string 00589 // target_tag_name, 00590 "", //"", // std::string 00591 // strVariables="", 00592 "", "", // std::string strInputData="", std::string strOutputData="", 00593 "", true, // std::string strNColName="", bool fOutputDouble=true, 00594 "", false, 00595 0.0, // std::string strPreserveVariables="", bool fPreserveAll=false, double 00596 // dFillValueOverride=0.0, 00597 runCtx->fInputConcave, 00598 runCtx->fOutputConcave // bool fInputConcave = false, bool fOutputConcave = false 00599 );MB_CHK_ERR( rval ); 00600 runCtx->timer_pop(); 00601 00602 // Invoke the CheckMap routine on the TempestRemap serial interface directly, if running 00603 // on a single process 00604 if( nprocs == 1 ) 00605 { 00606 const double dNormalTolerance = 1.0E-8; 00607 const double dStrictTolerance = 1.0E-12; 00608 weightMap->CheckMap( !runCtx->fNoCheck, !runCtx->fNoCheck, 00609 !runCtx->fNoCheck && ( runCtx->ensureMonotonicity ), dNormalTolerance, 00610 dStrictTolerance ); 00611 } 00612 00613 if( runCtx->outFilename.size() ) 00614 { 00615 // Write the map file to disk in parallel using either HDF5 or SCRIP interface 00616 rval = weightMap->WriteParallelMap( runCtx->outFilename.c_str() );MB_CHK_ERR( rval ); 00617 00618 // Write out the metadata information for the map file 00619 if( proc_id == 0 ) 00620 { 00621 size_t lastindex = runCtx->outFilename.find_last_of( "." ); 00622 sstr.str( "" ); 00623 sstr << runCtx->outFilename.substr( 0, lastindex ) << ".meta"; 00624 00625 std::ofstream metafile( sstr.str() ); 00626 metafile << "Generator = MOAB-TempestRemap (mbtempest) Offline Regridding " 00627 "Weight Generator" 00628 << std::endl; 00629 metafile << "domain_a = " << runCtx->inFilenames[0] << std::endl; 00630 metafile << "domain_b = " << runCtx->inFilenames[1] << std::endl; 00631 metafile << "domain_aNb = " 00632 << ( runCtx->intxFilename.size() ? runCtx->intxFilename : "outOverlap.h5m" ) << std::endl; 00633 metafile << "map_aPb = " << runCtx->outFilename << std::endl; 00634 metafile << "type_src = " << runCtx->disc_methods[0] << std::endl; 00635 metafile << "np_src = " << runCtx->disc_orders[0] << std::endl; 00636 metafile << "concave_src = " << ( runCtx->fInputConcave ? "true" : "false" ) << std::endl; 00637 metafile << "type_dst = " << runCtx->disc_methods[1] << std::endl; 00638 metafile << "np_dst = " << runCtx->disc_orders[1] << std::endl; 00639 metafile << "concave_dst = " << ( runCtx->fOutputConcave ? "true" : "false" ) << std::endl; 00640 metafile << "mono_type = " << runCtx->ensureMonotonicity << std::endl; 00641 metafile << "bubble = " << ( runCtx->fNoBubble ? "false" : "true" ) << std::endl; 00642 metafile << "version = " 00643 << "MOAB v5.1.0+" << std::endl; 00644 metafile.close(); 00645 } 00646 } 00647 00648 if( runCtx->verifyWeights ) 00649 { 00650 // Let us pick a sampling test function for solution evaluation 00651 moab::TempestOnlineMap::sample_function testFunction = 00652 &sample_stationary_vortex; // &sample_slow_harmonic; 00653 00654 runCtx->timer_push( "describe a solution on source grid" ); 00655 moab::Tag srcAnalyticalFunction; 00656 rval = weightMap->DefineAnalyticalSolution( srcAnalyticalFunction, "AnalyticalSolnSrcExact", 00657 moab::Remapper::SourceMesh, testFunction );MB_CHK_ERR( rval ); 00658 runCtx->timer_pop(); 00659 // rval = mbCore->write_file ( "srcWithSolnTag.h5m", NULL, writeOptions, 00660 // &runCtx->meshsets[0], 1 ); MB_CHK_ERR ( rval ); 00661 00662 runCtx->timer_push( "describe a solution on target grid" ); 00663 moab::Tag tgtAnalyticalFunction; 00664 moab::Tag tgtProjectedFunction; 00665 rval = weightMap->DefineAnalyticalSolution( tgtAnalyticalFunction, "AnalyticalSolnTgtExact", 00666 moab::Remapper::TargetMesh, testFunction, 00667 &tgtProjectedFunction, "ProjectedSolnTgt" );MB_CHK_ERR( rval ); 00668 // rval = mbCore->write_file ( "tgtWithSolnTag.h5m", NULL, writeOptions, 00669 // &runCtx->meshsets[1], 1 ); MB_CHK_ERR ( rval ); 00670 runCtx->timer_pop(); 00671 00672 runCtx->timer_push( "compute solution projection on target grid" ); 00673 rval = weightMap->ApplyWeights( srcAnalyticalFunction, tgtProjectedFunction );MB_CHK_ERR( rval ); 00674 runCtx->timer_pop(); 00675 rval = mbCore->write_file( "tgtWithSolnTag2.h5m", NULL, writeOptions, &runCtx->meshsets[1], 1 );MB_CHK_ERR( rval ); 00676 00677 runCtx->timer_push( "compute error metrics against analytical solution on target grid" ); 00678 std::map< std::string, double > errMetrics; 00679 rval = weightMap->ComputeMetrics( moab::Remapper::TargetMesh, tgtAnalyticalFunction, 00680 tgtProjectedFunction, errMetrics, true );MB_CHK_ERR( rval ); 00681 runCtx->timer_pop(); 00682 } 00683 00684 delete weightMap; 00685 } 00686 } 00687 00688 // Clean up 00689 remapper.clear(); 00690 delete runCtx; 00691 delete mbCore; 00692 00693 #ifdef MOAB_HAVE_MPI 00694 MPI_Finalize(); 00695 #endif 00696 exit( 0 ); 00697 } 00698 00699 static moab::ErrorCode CreateTempestMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh ) 00700 { 00701 moab::ErrorCode rval = moab::MB_SUCCESS; 00702 int err; 00703 moab::DebugOutput& outputFormatter = ctx.outputFormatter; 00704 00705 if( !ctx.proc_id ) { outputFormatter.printf( 0, "Creating TempestRemap Mesh object ...\n" ); } 00706 00707 if( ctx.meshType == moab::TempestRemapper::OVERLAP_FILES ) 00708 { 00709 // For the overlap method, choose between: "fuzzy", "exact" or "mixed" 00710 err = GenerateOverlapMesh( ctx.inFilenames[0], ctx.inFilenames[1], *tempest_mesh, ctx.outFilename, "NetCDF4", 00711 "exact", true ); 00712 00713 if( err ) { rval = moab::MB_FAILURE; } 00714 else 00715 { 00716 ctx.meshes.push_back( tempest_mesh ); 00717 } 00718 } 00719 else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MEMORY ) 00720 { 00721 // Load the meshes and validate 00722 ctx.meshsets.resize( 3 ); 00723 ctx.meshes.resize( 3 ); 00724 ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh ); 00725 ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh ); 00726 ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh ); 00727 00728 // First the source 00729 rval = remapper.LoadMesh( moab::Remapper::SourceMesh, ctx.inFilenames[0], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval ); 00730 ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh ); 00731 00732 // Next the target 00733 rval = remapper.LoadMesh( moab::Remapper::TargetMesh, ctx.inFilenames[1], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval ); 00734 ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh ); 00735 00736 // Now let us construct the overlap mesh, by calling TempestRemap interface directly 00737 // For the overlap method, choose between: "fuzzy", "exact" or "mixed" 00738 err = GenerateOverlapWithMeshes( *ctx.meshes[0], *ctx.meshes[1], *tempest_mesh, "" /*ctx.outFilename*/, 00739 "NetCDF4", "exact", false ); 00740 00741 if( err ) { rval = moab::MB_FAILURE; } 00742 else 00743 { 00744 remapper.SetMesh( moab::Remapper::OverlapMesh, tempest_mesh ); 00745 ctx.meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh ); 00746 // ctx.meshes.push_back(*tempest_mesh); 00747 } 00748 } 00749 else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MOAB ) 00750 { 00751 ctx.meshsets.resize( 3 ); 00752 ctx.meshes.resize( 3 ); 00753 ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh ); 00754 ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh ); 00755 ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh ); 00756 00757 const double radius_src = 1.0 /*2.0*acos(-1.0)*/; 00758 const double radius_dest = 1.0 /*2.0*acos(-1.0)*/; 00759 00760 const char* additional_read_opts = ( ctx.n_procs > 1 ? "NO_SET_CONTAINING_PARENTS;" : "" ); 00761 // Load the source mesh and validate 00762 rval = remapper.LoadNativeMesh( ctx.inFilenames[0], ctx.meshsets[0], additional_read_opts );MB_CHK_ERR( rval ); 00763 // Rescale the radius of both to compute the intersection 00764 rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[0], radius_src );MB_CHK_ERR( rval ); 00765 rval = remapper.ConvertMeshToTempest( moab::Remapper::SourceMesh );MB_CHK_ERR( rval ); 00766 ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh ); 00767 00768 // Load the target mesh and validate 00769 rval = remapper.LoadNativeMesh( ctx.inFilenames[1], ctx.meshsets[1], additional_read_opts );MB_CHK_ERR( rval ); 00770 rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[1], radius_dest );MB_CHK_ERR( rval ); 00771 rval = remapper.ConvertMeshToTempest( moab::Remapper::TargetMesh );MB_CHK_ERR( rval ); 00772 ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh ); 00773 } 00774 else if( ctx.meshType == moab::TempestRemapper::ICO ) 00775 { 00776 err = GenerateICOMesh( *tempest_mesh, ctx.blockSize, ctx.computeDual, ctx.outFilename, "NetCDF4" ); 00777 00778 if( err ) { rval = moab::MB_FAILURE; } 00779 else 00780 { 00781 ctx.meshes.push_back( tempest_mesh ); 00782 } 00783 } 00784 else if( ctx.meshType == moab::TempestRemapper::RLL ) 00785 { 00786 err = GenerateRLLMesh( *tempest_mesh, // Mesh& meshOut, 00787 ctx.blockSize * 2, ctx.blockSize, // int nLongitudes, int nLatitudes, 00788 0.0, 360.0, // double dLonBegin, double dLonEnd, 00789 -90.0, 90.0, // double dLatBegin, double dLatEnd, 00790 false, false, false, // bool fGlobalCap, bool fFlipLatLon, bool fForceGlobal, 00791 "" /*ctx.inFilename*/, "", 00792 "", // std::string strInputFile, std::string strInputFileLonName, std::string 00793 // strInputFileLatName, 00794 ctx.outFilename, "NetCDF4", // std::string strOutputFile, std::string strOutputFormat 00795 true // bool fVerbose 00796 ); 00797 00798 if( err ) { rval = moab::MB_FAILURE; } 00799 else 00800 { 00801 ctx.meshes.push_back( tempest_mesh ); 00802 } 00803 } 00804 else // default 00805 { 00806 err = GenerateCSMesh( *tempest_mesh, ctx.blockSize, ctx.outFilename, "NetCDF4" ); 00807 00808 if( err ) { rval = moab::MB_FAILURE; } 00809 else 00810 { 00811 ctx.meshes.push_back( tempest_mesh ); 00812 } 00813 } 00814 00815 if( ctx.meshType != moab::TempestRemapper::OVERLAP_MOAB && !tempest_mesh ) 00816 { 00817 std::cout << "Tempest Mesh is not a complete object; Quitting..."; 00818 exit( -1 ); 00819 } 00820 00821 return rval; 00822 } 00823 00824 /////////////////////////////////////////////// 00825 // Test functions 00826 00827 double sample_slow_harmonic( double dLon, double dLat ) 00828 { 00829 return ( 2.0 + cos( dLat ) * cos( dLat ) * cos( 2.0 * dLon ) ); 00830 } 00831 00832 double sample_fast_harmonic( double dLon, double dLat ) 00833 { 00834 return ( 2.0 + pow( sin( 2.0 * dLat ), 16.0 ) * cos( 16.0 * dLon ) ); 00835 // return (2.0 + pow(cos(2.0 * dLat), 16.0) * cos(16.0 * dLon)); 00836 } 00837 00838 double sample_constant( double /*dLon*/, double /*dLat*/ ) 00839 { 00840 return 1.0; 00841 } 00842 00843 double sample_stationary_vortex( double dLon, double dLat ) 00844 { 00845 const double dLon0 = 0.0; 00846 const double dLat0 = 0.6; 00847 const double dR0 = 3.0; 00848 const double dD = 5.0; 00849 const double dT = 6.0; 00850 00851 /// Find the rotated longitude and latitude of a point on a sphere 00852 /// with pole at (dLonC, dLatC). 00853 { 00854 double dSinC = sin( dLat0 ); 00855 double dCosC = cos( dLat0 ); 00856 double dCosT = cos( dLat ); 00857 double dSinT = sin( dLat ); 00858 00859 double dTrm = dCosT * cos( dLon - dLon0 ); 00860 double dX = dSinC * dTrm - dCosC * dSinT; 00861 double dY = dCosT * sin( dLon - dLon0 ); 00862 double dZ = dSinC * dSinT + dCosC * dTrm; 00863 00864 dLon = atan2( dY, dX ); 00865 if( dLon < 0.0 ) { dLon += 2.0 * M_PI; } 00866 dLat = asin( dZ ); 00867 } 00868 00869 double dRho = dR0 * cos( dLat ); 00870 double dVt = 3.0 * sqrt( 3.0 ) / 2.0 / cosh( dRho ) / cosh( dRho ) * tanh( dRho ); 00871 00872 double dOmega; 00873 if( dRho == 0.0 ) { dOmega = 0.0; } 00874 else 00875 { 00876 dOmega = dVt / dRho; 00877 } 00878 00879 return ( 1.0 - tanh( dRho / dD * sin( dLon - dOmega * dT ) ) ); 00880 } 00881 00882 ///////////////////////////////////////////////