Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ExchangeHalos.cpp
Go to the documentation of this file.
1 /** @example ExchangeHalos Driver
2  * \brief Example program that shows the use case for performing tag data exchange
3  * between parallel processors in order to sync data on shared entities.
4  *
5  * <b>This example </b>:
6  * -# Initialize MPI and instantiates MOAB
7  * -# Gets user options: Input mesh file name, vector tag length, ghost layer size etc
8  * -# Create the root and partition sets
9  * -# Instantiate ParallelComm and read the mesh file in parallel using appropriate options
10  * -# Create the required number of ghost layers as requested by the user (default = 3)
11  * -# Get 2D MPAS polygonal entities in the mesh and filter to get only the "owned" entities
12  * -# Create two tags: scalar_variable (single data/cell) and vector_variable (multiple data/cell)
13  * -# Set tag data using analytical functions for both scalar and vector fields on owned entities
14  * -# Exchange shared entity information and tags between processors
15  * -# If debugging is turned on, store mesh file and tag on root process (will not contain data on shared entities)
16  * -# Perform exchange of scalar tag data on shared entities
17  * -# Perform exchange of vector tag data on shared entities
18  * -# If debugging is turned on, store mesh file and tag on root process (will now contain data on *all* entities)
19  * -# Destroy the MOAB instance and finalize MPI
20  *
21  * <b>To run: </b>
22  * mpiexec -n np ./ExchangeHalos --input <mpas_mesh_file> --nghosts <ghostlayers> --vtaglength <vector component
23  * size> \
24  * --nexchanges <number of exchange runs>
25  * <b>Example:</b>
26  * mpiexec -n 4 ./ExchangeHalos --input $MOAB_DIR/MeshFiles/unittest/io/mpasx1.642.t.2.nc --nghosts 3 --vtaglength 100
27  *
28  * NOTE: --debug option can be added to write out extra files in h5m format to visualize some output (written from root
29  * task only)
30  *
31  */
32 
33 // MOAB includes
34 #include "moab/Core.hpp"
35 #include "moab/CpuTimer.hpp"
36 #include "moab/ProgOptions.hpp"
37 
38 #ifndef MOAB_HAVE_MPI
39 #error "Please build MOAB with MPI..."
40 #endif
41 
42 #include "moab/ParallelComm.hpp"
43 #include "MBParallelConventions.h"
44 
45 // C++ includes
46 #include <iostream>
47 #include <string>
48 
49 using namespace moab;
50 using namespace std;
51 #define dbgprint( MSG ) \
52  do \
53  { \
54  if( context.proc_id == 0 ) std::cout << MSG << std::endl; \
55  } while( false )
56 
57 #define runchk( CODE, MSG ) \
58  do \
59  { \
60  moab::ErrorCode err = CODE; \
61  MB_CHK_SET_ERR( err, MSG ); \
62  } while( false )
63 
64 #define runchk_cont( CODE, MSG ) \
65  do \
66  { \
67  moab::ErrorCode err = CODE; \
68  MB_CHK_ERR_CONT( err ); \
69  if( err ) std::cout << "Error:: " << MSG << std::endl; \
70  } while( false )
71 
72 // Run time context structure
73 
74 /// @brief The RunttimeContext is an example specific class to store
75 /// the run specific input data, MOAB datastructures used during the run
76 /// and provides other utility functions to profile operations etc
78 {
79  public:
80  int dimension{ 2 }; /// dimension of the problem
81  std::string input_filename; /// input file name (nc format)
82  std::string output_filename; /// output file name (h5m format)
83  int ghost_layers{ 3 }; /// number of ghost layers
84  std::string scalar_tagname; /// scalar tag name
85  std::string vector_tagname; /// vector tag name
86  int vector_length{ 3 }; /// length of the vector tag components
87  int num_max_exchange{ 10 }; /// total number of exchange iterations
88  bool debug_output{ false }; /// write debug output information?
89  int proc_id{ 1 }; /// process identifier
90  int num_procs{ 1 }; /// total number of processes
91  double last_counter{ 0.0 }; /// last time counter between push/pop timer
92 
93  // MOAB objects
96  moab::EntityHandle fileset{ 0 }, partnset{ 0 };
97 
98  /// @brief Constructor: allocate MOAB interface and communicator, and initialize
99  /// other data members with some default values
100  explicit RuntimeContext( MPI_Comm comm = MPI_COMM_WORLD );
101 
102  /// @brief Destructor: deallocate MOAB interface and communicator
103  ~RuntimeContext();
104 
105  /// @brief Parse the runtime command line options
106  /// @param argc - number of command line arguments
107  /// @param argv - command line arguments as string list
108  void ParseCLOptions( int argc, char* argv[] );
109 
110  /// @brief Measure and start the timer to profile a task
111  /// @param operation String name of the task being measured
112  inline void timer_push( const std::string& operation );
113 
114  /// @brief Stop the timer and store the elapsed duration
115  /// @param nruns Optional argument used to average the measured time
116  void timer_pop( const int nruns = 1 );
117 
118  /// @brief Return the last elapsed time
119  /// @return last_counter from timer_pop was called
120  inline double last_elapsed() const;
121 
122  /// @brief Load a MOAB supported file (h5m or nc format) from disk
123  /// representing an MPAS mesh
124  /// @param load_ghosts Optional boolean to specify whether to load ghosts
125  /// when reading the file (only relevant for h5m)
126  /// @return Error code if any (else MB_SUCCESS)
127  moab::ErrorCode load_file( bool load_ghosts = false );
128 
129  /// @brief Create scalar and vector tags in the MOAB mesh instance
130  /// @param tagScalar Tag reference to the scalar field
131  /// @param tagVector Tag reference to the vector field
132  /// @param entities Entities on which both the scalar and vector fields are defined
133  /// @return Error code if any (else MB_SUCCESS)
134  moab::ErrorCode create_sv_tags( moab::Tag& tagScalar, moab::Tag& tagVector, moab::Range& entities );
135 
136  /// @brief Evaluate some closed-form Spherical Harmonic functions with an optional multiplier term
137  /// @param lon Longitude in lat-lon space
138  /// @param lat Latitude in lat-lon space
139  /// @param type Function type
140  /// @param multiplier Optional multiplier to scale value (default=1.0)
141  /// @return value of the evaluated function
142  inline double evaluate_function( double lon, double lat, int type = 1, double multiplier = 1.0 ) const
143  {
144  switch( type )
145  {
146  case 1:
147  return ( 2.0 + std::pow( sin( 2.0 * lat ), 16.0 ) * cos( 16.0 * lon ) ) * multiplier;
148  default:
149  return ( 2.0 + cos( lon ) * cos( lon ) * cos( 2.0 * lat ) ) * multiplier;
150  }
151  }
152 
153  private:
154  /// @brief Compute the centroids of elements in 2D lat/lon space
155  /// @param entities Entities to compute centroids
156  /// @return Vector of centroids (as lat/lon)
157  std::vector< double > compute_centroids( const moab::Range& entities ) const;
158 
160  double mTimerOps{ 0.0 };
161  std::string mOpName;
162 };
163 
164 //
165 // Start of main test program
166 //
167 int main( int argc, char** argv )
168 {
169  // Initialize MPI first
170  MPI_Init( &argc, &argv );
171 
172  {
173  // Create our context for this example run
175  dbgprint( "********** Exchange halos example **********\n" );
176 
177  // Get the input options
178  context.ParseCLOptions( argc, argv );
179 
180  /////////////////////////////////////////////////////////////////////////
181  // Print out the input parameters in use
182  dbgprint( " -- Input Parameters -- " );
183  dbgprint( " Number of Processes = " << context.num_procs );
184  dbgprint( " Input mesh = " << context.input_filename );
185  dbgprint( " Ghost Layers = " << context.ghost_layers );
186  dbgprint( " Scalar Tag name = " << context.scalar_tagname );
187  dbgprint( " Vector Tag name = " << context.vector_tagname );
188  dbgprint( " Vector Tag length = " << context.vector_length << endl );
189  /////////////////////////////////////////////////////////////////////////
190 
191  // Timer storage for all phases
192  double elapsed_times[4];
193 
194  // Read the input file specified by user, in parallel, using appropriate options
195  // Supports reading partitioned h5m files and MPAS nc files directly with online Zoltan partitioning
196  context.timer_push( "Read input file" );
197  {
198  // Load the file from disk with given options
199  runchk( context.load_file( false ), "MOAB::load_file failed for filename: " << context.input_filename );
200  }
201  context.timer_pop();
202  elapsed_times[0] = context.last_elapsed();
203 
204  // Let the actual measurements begin...
205  dbgprint( "\n- Starting execution -\n" );
206 
207  // We need to set up the ghost layers requested by the user. First correct for thin layers and then
208  // call `exchange_ghost_cells` to prepare the mesh for use with halo regions
209  context.timer_push( "Setup ghost layers" );
210  {
211  // Loop over the number of ghost layers needed and ask MOAB for layers 1 at a time
212  for( int ighost = 0; ighost < context.ghost_layers; ++ighost )
213  {
214  // Exchange ghost cells
215  int ghost_dimension = context.dimension;
216  int bridge_dimension = context.dimension - 1;
217  // Let us now get all ghost layers from adjacent parts
218  runchk( context.parallel_communicator->exchange_ghost_cells(
219  ghost_dimension, bridge_dimension, ( ighost + 1 ), 0, true /* store_remote_handles */,
220  true /* wait_all */, &context.fileset ),
221  "Exchange ghost cells failed" ); // true to store remote handles
222 
223  // Ensure that all processes understand about multi-shared vertices and entities
224  // in case some adjacent parts are only m layers thick (where m < context.ghost_layers)
225  if( ighost < context.ghost_layers - 1 )
226  runchk( context.parallel_communicator->correct_thin_ghost_layers(),
227  "Thin layer correction failed" );
228  }
229  }
230  context.timer_pop();
231  elapsed_times[1] = context.last_elapsed();
232 
233  // Get the 2D MPAS elements and filter it so that we have only owned elements
234  Range dimEnts;
235  {
236  // Get all entities of dimension = dim
237  runchk( context.moab_interface.get_entities_by_dimension( context.fileset, context.dimension, dimEnts ),
238  "Getting 2D entities failed" );
239  // Get only owned entities! The ghosted/shared entities will get their data when we exchange
240  // So let us filter entities based on the status: NOT x NOT_OWNED = OWNED status :-)
241  runchk( context.parallel_communicator->filter_pstatus( dimEnts, PSTATUS_NOT_OWNED, PSTATUS_NOT ),
242  "Filtering pstatus failed" );
243 
244  // Aggregate the total number of elements in the mesh
245  auto numEntities = dimEnts.size();
246  int numTotalEntities = 0;
247  MPI_Reduce( &numEntities, &numTotalEntities, 1, MPI_INT, MPI_SUM, 0,
248  context.parallel_communicator->proc_config().proc_comm() );
249 
250  // We expect the total number of elements to be constant, immaterial of number of processes.
251  // If not, we have a bug!
252  dbgprint( "Total number of " << context.dimension << "D elements in the mesh = " << numTotalEntities );
253  }
254 
255  Tag tagScalar = nullptr;
256  Tag tagVector = nullptr;
257  // Create two tag handles: scalar_variable and vector_variable
258  // Set these tags with appropriate closed form functional data
259  // based on element centroid information
260  runchk( context.create_sv_tags( tagScalar, tagVector, dimEnts ), "Unable to create scalar and vector tags" );
261 
262  // let us write out the local mesh before tag_exchange is called
263  // we expect to see data only on the owned entities - and ghosted entities should have default values
264  if( context.debug_output && ( context.proc_id == 0 ) ) // only on root process, for debugging
265  {
266  dbgprint( "> Writing to file *before* ghost exchange " );
267  runchk( context.moab_interface.write_file( "exchangeHalos_output_rank0_pre.h5m", "H5M", "" ),
268  "Writing to disk failed" );
269  }
270 
271  // Perform exchange of tag data between neighboring tasks
272  dbgprint( "> Exchanging tags between processors " );
273  context.timer_push( "Exchange scalar tag data" );
274  for( auto irun = 0; irun < context.num_max_exchange; ++irun )
275  {
276  // Exchange scalar tags between processors
277  runchk( context.parallel_communicator->exchange_tags( tagScalar, dimEnts ),
278  "Exchanging scalar tag between processors failed" );
279  }
280  context.timer_pop( context.num_max_exchange );
281  elapsed_times[2] = context.last_elapsed();
282 
283  context.timer_push( "Exchange vector tag data" );
284  for( auto irun = 0; irun < context.num_max_exchange; ++irun )
285  {
286  // Exchange vector tags between processors
287  runchk( context.parallel_communicator->exchange_tags( tagVector, dimEnts ),
288  "Exchanging vector tag between processors failed" );
289  }
290  context.timer_pop( context.num_max_exchange );
291  elapsed_times[3] = context.last_elapsed();
292 
293  // let us write out the local mesh after tag_exchange is called
294  // we expect to see real data on both owned and ghost entities in halo regions (non-default values)
295  if( context.debug_output && ( context.proc_id == 0 ) ) // only on root process, for debugging
296  {
297  dbgprint( "> Writing to file *after* ghost exchange " );
298  runchk( context.moab_interface.write_file( "exchangeHalos_output_rank0_post.h5m", "H5M", "" ),
299  "Writing to disk failed" );
300  }
301 
302  // Write out the final mesh with the tag data and mesh -- just for verification
303  if( context.debug_output )
304  {
305  dbgprint( "> Writing out the final mesh and data in MOAB h5m format. File = " << context.output_filename );
306  string write_options = ( context.num_procs > 1 ? "PARALLEL=WRITE_PART;DEBUG_IO=0;" : "" );
307  // Write out to output file to visualize reduction/exchange of tag data
308  runchk( context.moab_interface.write_file( context.output_filename.c_str(), "H5M", write_options.c_str() ),
309  "File write failed" );
310  }
311 
312  // Consolidated timing results: the data is listed as follows
313  // [ntasks, nghosts, load_mesh(I/O), exchange_ghost_cells(setup), exchange_tags(scalar),
314  // exchange_tags(vector)]
315  dbgprint( "\n> Consolidated: [" << context.num_procs << ", " << context.ghost_layers << ", " << elapsed_times[0]
316  << ", " << elapsed_times[1] << ", " << elapsed_times[2] << ", "
317  << elapsed_times[3] << "]," );
318 
319  // execution finished
320  dbgprint( "\n********** ExchangeHalos Example DONE! **********" );
321  }
322  // Done, cleanup
323  MPI_Finalize();
324 
325  return 0;
326 }
327 
328 /// Implementation details ///
329 
330 ///
332  : input_filename( std::string( MESH_DIR ) + std::string( "/io/mpasx1.642.t.2.nc" ) ),
333  output_filename( "exchangeHalos_output.h5m" ), scalar_tagname( "scalar_variable" ),
334  vector_tagname( "vector_variable" )
335 {
336  // Create sets for the mesh and partition. Then pass these to the load_file functions to populate the mesh.
337  runchk_cont( moab_interface.create_meshset( moab::MESHSET_SET, fileset ), "Creating root set failed" );
338  runchk_cont( moab_interface.create_meshset( moab::MESHSET_SET, partnset ), "Creating partition set failed" );
339 
340  // Create the parallel communicator object with the partition handle associated with MOAB
342 
345 }
346 
348 {
349  delete parallel_communicator;
350 }
351 
352 void RuntimeContext::ParseCLOptions( int argc, char* argv[] )
353 {
354  ProgOptions opts;
355  // Input mesh
356  opts.addOpt< std::string >( "input", "Input mesh filename to load in parallel. Default=data/default_mesh_holes.h5m",
357  &input_filename );
358  // Output mesh
359  opts.addOpt< void >( "debug", "Should we write output file? Default=false", &debug_output );
360  opts.addOpt< std::string >( "output",
361  "Output mesh filename for verification (use --debug). Default=exchangeHalos_output.h5m",
362  &output_filename );
363  // Dimension of the input mesh
364  // Vector tag length
365  opts.addOpt< int >( "vtaglength", "Size of vector components per each entity. Default=3", &vector_length );
366  // Number of halo (ghost) regions
367  opts.addOpt< int >( "nghosts", "Number of ghost layers (halos) to exchange. Default=3", &ghost_layers );
368  // Number of times to perform the halo exchange for timing
369  opts.addOpt< int >( "nexchanges", "Number of ghost-halo exchange iterations to perform. Default=10",
370  &num_max_exchange );
371 
372  opts.parseCommandLine( argc, argv );
373 }
374 
375 void RuntimeContext::timer_push( const std::string& operation )
376 {
378  mOpName = operation;
379 }
380 
381 void RuntimeContext::timer_pop( const int nruns )
382 {
383  double locElapsed = mTimer.time_since_birth() - mTimerOps;
384  double avgElapsed = 0;
385  double maxElapsed = 0;
386  MPI_Reduce( &locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, parallel_communicator->comm() );
387  MPI_Reduce( &locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, parallel_communicator->comm() );
388  if( proc_id == 0 )
389  {
390  avgElapsed /= num_procs;
391  if( nruns > 1 )
392  std::cout << "[LOG] Time taken to " << mOpName.c_str() << ", averaged over " << nruns
393  << " runs : max = " << maxElapsed / nruns << ", avg = " << avgElapsed / nruns << "\n";
394  else
395  std::cout << "[LOG] Time taken to " << mOpName.c_str() << " : max = " << maxElapsed
396  << ", avg = " << avgElapsed << "\n";
397 
398  last_counter = maxElapsed / nruns;
399  }
400  mOpName.clear();
401 }
402 
404 {
405  return last_counter;
406 }
407 
409 {
410  // Get element (centroid) coordinates so that we can evaluate some arbitrary data
411  std::vector< double > entCoords = compute_centroids( entities ); // [entities * [lon, lat]]
412 
413  if( proc_id == 0 ) std::cout << "> Getting scalar tag handle " << scalar_tagname << "..." << std::endl;
414  double defSTagValue = -1.0;
415  bool createdTScalar = false;
416  // Get or create the scalar tag: default name = "scalar_variable"
417  // Type: double, Components: 1, Layout: Dense (all entities potentially), Default: -1.0
419  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &defSTagValue, &createdTScalar ),
420  "Retrieving scalar tag handle failed" );
421 
422  // we expect to create a new tag -- fail if Tag already exists since we do not want to overwrite data
423  assert( createdTScalar );
424  // set the data for scalar tag with an analytical Spherical Harmonic function
425  {
426  std::vector< double > tagValues( entities.size(), -1.0 );
427  std::generate( tagValues.begin(), tagValues.end(), [=, &entCoords]() {
428  static int index = 0;
429  const int offset = index++ * 2;
430  return evaluate_function( entCoords[offset], entCoords[offset + 1] );
431  } );
432  // Set local scalar tag data for exchange
433  runchk( moab_interface.tag_set_data( tagScalar, entities, tagValues.data() ),
434  "Setting scalar tag data failed" );
435  }
436 
437  if( proc_id == 0 ) std::cout << "> Getting vector tag handle " << vector_tagname << "..." << std::endl;
438  std::vector< double > defVTagValue( vector_length, -1.0 );
439  bool createdTVector = false;
440  // Get or create the scalar tag: default name = "vector_variable"
441  // Type: double, Components: vector_length, Layout: Dense (all entities potentially), Default: [-1.0,..]
443  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, defVTagValue.data(),
444  &createdTVector ),
445  "Retrieving vector tag handle failed" );
446 
447  // we expect to create a new tag -- fail if Tag already exists since we do not want to overwrite data
448  assert( createdTVector );
449  // set the data for vector tag with an analytical Spherical Harmonic function
450  // with an optional scaling for each component; just to make it look different :-)
451  {
452  const int veclength = vector_length;
453  std::vector< double > tagValues( entities.size() * veclength, -1.0 );
454  std::generate( tagValues.begin(), tagValues.end(), [=, &entCoords]() {
455  static int index = 0;
456  const int offset = ( index++ / veclength ) * 2;
457  return this->evaluate_function( entCoords[offset], entCoords[offset + 1], 2, ( index % veclength + 1.0 ) );
458  } );
459  // Set local tag data for exchange
460  runchk( moab_interface.tag_set_data( tagVector, entities, tagValues.data() ),
461  "Setting vector tag data failed" );
462  }
463 
464  return moab::MB_SUCCESS;
465 }
466 
468 {
469  /// Parallel Read options:
470  /// PARALLEL = type {READ_PART} : Read on all tasks
471  /// PARTITION_METHOD = RCBZOLTAN : Use Zoltan partitioner to compute an online partition and redistribute on the
472  /// fly PARTITION = PARALLEL_PARTITION : Partition as you read based on part information stored in h5m file
473  /// PARALLEL_RESOLVE_SHARED_ENTS : Communicate to all processors to get the shared adjacencies
474  /// consistently in parallel
475  /// PARALLEL_GHOSTS : a.b.c
476  /// : a = 2 - highest dimension of entities (2D in this case)
477  /// : b = 1 - dimension of entities to calculate adjacencies (vertex=0, edges=1)
478  /// : c = 3 - number of ghost layers needed (3 in this case)
479  std::string read_options = "DEBUG_IO=0;";
480  std::string::size_type idx = input_filename.rfind( '.' );
481  if( num_procs > 1 && idx != std::string::npos )
482  {
483  std::string extension = input_filename.substr( idx + 1 );
484  if( !extension.compare( "nc" ) )
485  // PARTITION_METHOD= [RCBZOLTAN, TRIVIAL]
486  read_options += "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;"
487  "PARALLEL_RESOLVE_SHARED_ENTS;NO_EDGES;NO_MIXED_ELEMENTS;VARIABLE=;";
488  else if( !extension.compare( "h5m" ) )
489  read_options +=
490  "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;"
491  "PARALLEL_RESOLVE_SHARED_ENTS;" +
492  ( load_ghosts ? "PARALLEL_THIN_GHOST_LAYER;PARALLEL_GHOSTS=2.1." + std::to_string( ghost_layers ) + ";"
493  : "" );
494  else
495  {
496  std::cout << "Error unsupported file type (only h5m and nc) for this example: " << input_filename
497  << std::endl;
499  }
500  }
501 
502  // Load the file from disk with given read options in parallel and associate all entities to fileset
503  return moab_interface.load_file( input_filename.c_str(), &fileset, read_options.c_str() );
504 }
505 
506 std::vector< double > RuntimeContext::compute_centroids( const moab::Range& entities ) const
507 {
508  double node[3];
509  std::vector< double > eCentroids( entities.size() * 2 ); // [lon, lat]
510  size_t offset = 0;
511  for( auto entity : entities )
512  {
513  // Get the element coordinates (centroid) on the real mesh
514  runchk_cont( moab_interface.get_coords( &entity, 1, node ), "Getting entity coordinates failed" );
515 
516  // scale by magnitude so that element is on unit sphere
517  double magnitude = std::sqrt( node[0] * node[0] + node[1] * node[1] + node[2] * node[2] );
518  node[0] /= magnitude;
519  node[1] /= magnitude;
520  node[2] /= magnitude;
521 
522  // compute the spherical transformation onto unit sphere
523  eCentroids[offset] = atan2( node[1], node[0] );
524  if( eCentroids[offset] < 0.0 ) eCentroids[offset] += 2.0 * M_PI;
525  eCentroids[offset + 1] = asin( node[2] );
526 
527  offset += 2; // increment the offset
528  }
529  // return centroid list for elements
530  return eCentroids;
531 }