MOAB: Mesh Oriented datABase  (version 5.3.1)
parallel_write_test.cpp
Go to the documentation of this file.
00001 #include "moab/Core.hpp"
00002 #include "moab/ParallelComm.hpp"
00003 #include "MBTagConventions.hpp"
00004 #include "moab_mpi.h"
00005 #include <cstdlib>
00006 #include <iostream>
00007 #include <ctime>
00008 #include <cstring>
00009 #include <cmath>
00010 #include <cassert>
00011 #include <cstdio>
00012 #include <sstream>
00013 
00014 using namespace moab;
00015 
00016 #define TPRINT( A ) tprint( ( A ) )
00017 static void tprint( const char* A )
00018 {
00019     int rank;
00020     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00021     char buffer[128];
00022     sprintf( buffer, "%02d: %6.2f: %s\n", rank, (double)clock() / CLOCKS_PER_SEC, A );
00023     fputs( buffer, stderr );
00024 }
00025 
00026 const int DEFAULT_INTERVALS   = 2;
00027 const char* DEFAULT_FILE_NAME = "parallel_write_test.h5m";
00028 
00029 // Create mesh for each processor that is a cube of hexes
00030 // with the specified interval count along each edge.  Cubes
00031 // of mesh will be positioned and vertex IDs assigned such that
00032 // there are shared entities.  If the cubic root of the number
00033 // of processors is a whole number, then each processors mesh
00034 // will be a cube in a grid with that number of processor blocks
00035 // along each edge.  Otherwise processor blocks will be arranged
00036 // within the subset of the grid that is the ceiling of the cubic
00037 // root of the comm size such that there are no disjoint regions.
00038 ErrorCode generate_mesh( Interface& moab, int intervals );
00039 
00040 const char args[] = "[-i <intervals>] [-o <filename>] [-L <filename>] [-g <n>]";
00041 void help()
00042 {
00043     std::cout << "parallel_write_test " << args << std::endl
00044               << "  -i <N>    Each processor owns an NxNxN cube of hex elements (default: " << DEFAULT_INTERVALS << ")"
00045               << std::endl
00046               << "  -o <name> Retain output file and name it as specified." << std::endl
00047               << "  -L <name> Write local mesh to file name prefixed with MPI rank" << std::endl
00048               << "  -g <n>    Specify writer debug output level" << std::endl
00049               << "  -R        Skip resolve of shared entities (interface ents will be duplicated in file)" << std::endl
00050               << std::endl
00051               << "This program creates a (non-strict) subset of a regular hex mesh "
00052                  "such that the mesh is already partitioned, and then attempts to "
00053                  "write that mesh using MOAB's parallel HDF5 writer.  The mesh size "
00054                  "will scale with the number of processors and the number of elements "
00055                  "per processor (the latter is a function of the value specified "
00056                  "with the '-i' flag.)"
00057               << std::endl
00058               << std::endl
00059               << "Let N = ceil(cbrt(P)), where P is the number of processes.  "
00060                  "The mesh will be some subset of a cube with one corner at the "
00061                  "origin and the other at (N,N,N).  Each processor will own a "
00062                  "non-overlapping 1x1x1 unit block of mesh within that cube.  "
00063                  "If P is a power of 3, then the entire NxNxN cube will be "
00064                  "filled with hex elements.  Otherwise, some connected subset "
00065                  "of the cube will be meshed.  Each processor is assigned a "
00066                  "sub-block of the cube by rank where the blocks are enumerated "
00067                  "sequentally with x increasing most rapidly and z least rapidly."
00068               << std::endl
00069               << std::endl
00070               << "The size of the mesh owned by each processor is controlled by "
00071                  "the number of intervals along each edge of its block of mesh.  "
00072                  "If each block has N intervals, than each processor will have "
00073                  "N^3 hex elements."
00074               << std::endl
00075               << std::endl;
00076 }
00077 
00078 int main( int argc, char* argv[] )
00079 {
00080     int ierr = MPI_Init( &argc, &argv );
00081     if( ierr )
00082     {
00083         std::cerr << "MPI_Init failed with error code: " << ierr << std::endl;
00084         return ierr;
00085     }
00086     int rank;
00087     ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00088     if( ierr )
00089     {
00090         std::cerr << "MPI_Comm_rank failed with error code: " << ierr << std::endl;
00091         return ierr;
00092     }
00093     int size;
00094     ierr = MPI_Comm_size( MPI_COMM_WORLD, &size );
00095     if( ierr )
00096     {
00097         std::cerr << "MPI_Comm_size failed with error code: " << ierr << std::endl;
00098         return ierr;
00099     }
00100 
00101     // settings controlled by CL flags
00102     const char* output_file_name = 0;
00103     const char* indiv_file_name  = 0;
00104     int intervals = 0, debug_level = 0;
00105     // state for CL flag processing
00106     bool expect_intervals    = false;
00107     bool expect_file_name    = false;
00108     bool expect_indiv_file   = false;
00109     bool skip_resolve_shared = false;
00110     bool expect_debug_level  = false;
00111     // process CL args
00112     for( int i = 1; i < argc; ++i )
00113     {
00114         if( expect_intervals )
00115         {
00116             char* endptr = 0;
00117             intervals    = (int)strtol( argv[i], &endptr, 0 );
00118             if( *endptr || intervals < 1 )
00119             {
00120                 std::cerr << "Invalid block interval value: " << argv[i] << std::endl;
00121                 return 1;
00122             }
00123             expect_intervals = false;
00124         }
00125         else if( expect_indiv_file )
00126         {
00127             indiv_file_name   = argv[i];
00128             expect_indiv_file = false;
00129         }
00130         else if( expect_file_name )
00131         {
00132             output_file_name = argv[i];
00133             expect_file_name = false;
00134         }
00135         else if( expect_debug_level )
00136         {
00137             debug_level = atoi( argv[i] );
00138             if( debug_level < 1 )
00139             {
00140                 std::cerr << "Invalid argument following -g flag: \"" << argv[i] << '"' << std::endl;
00141                 return 1;
00142             }
00143             expect_debug_level = false;
00144         }
00145         else if( !strcmp( "-i", argv[i] ) )
00146             expect_intervals = true;
00147         else if( !strcmp( "-o", argv[i] ) )
00148             expect_file_name = true;
00149         else if( !strcmp( "-L", argv[i] ) )
00150             expect_indiv_file = true;
00151         else if( !strcmp( "-R", argv[i] ) )
00152             skip_resolve_shared = true;
00153         else if( !strcmp( "-g", argv[i] ) )
00154             expect_debug_level = true;
00155         else if( !strcmp( "-h", argv[i] ) )
00156         {
00157             help();
00158             return 0;
00159         }
00160         else
00161         {
00162             std::cerr << "Unexpected argument: " << argv[i] << std::endl
00163                       << "Usage: " << argv[0] << " " << args << std::endl
00164                       << "       " << argv[0] << " -h" << std::endl
00165                       << " Try '-h' for help." << std::endl;
00166             return 1;
00167         }
00168     }
00169     // Check for missing argument after last CL flag
00170     if( expect_file_name || expect_intervals || expect_indiv_file )
00171     {
00172         std::cerr << "Missing argument for '" << argv[argc - 1] << "'" << std::endl;
00173         return 1;
00174     }
00175     // If intervals weren't specified, use default
00176     if( intervals == 0 )
00177     {
00178         std::cout << "Using default interval count: " << DEFAULT_INTERVALS << std::endl;
00179         intervals = DEFAULT_INTERVALS;
00180     }
00181     // If no output file was specified, use default one and note that
00182     // we need to delete it when the test completes.
00183     bool keep_output_file = true;
00184     if( !output_file_name )
00185     {
00186         output_file_name = DEFAULT_FILE_NAME;
00187         keep_output_file = false;
00188     }
00189 
00190     // Create mesh
00191     TPRINT( "Generating mesh" );
00192     double gen_time = MPI_Wtime();
00193     Core mb;
00194     Interface& moab = mb;
00195     ErrorCode rval  = generate_mesh( moab, intervals );
00196     if( MB_SUCCESS != rval )
00197     {
00198         std::cerr << "Mesh creation failed with error code: " << rval << std::endl;
00199         return (int)rval;
00200     }
00201     gen_time = MPI_Wtime() - gen_time;
00202 
00203     // Write out local mesh on each processor if requested.
00204     if( indiv_file_name )
00205     {
00206         TPRINT( "Writing individual file" );
00207         char buffer[64];
00208         int width = (int)ceil( log10( size ) );
00209         sprintf( buffer, "%0*d-", width, rank );
00210         std::string name( buffer );
00211         name += indiv_file_name;
00212         rval = moab.write_file( name.c_str() );
00213         if( MB_SUCCESS != rval )
00214         {
00215             std::cerr << "Failed to write file: " << name << std::endl;
00216             return (int)rval;
00217         }
00218     }
00219 
00220     double res_time = MPI_Wtime();
00221     Range hexes;
00222     moab.get_entities_by_type( 0, MBHEX, hexes );
00223     if( !skip_resolve_shared )
00224     {
00225         TPRINT( "Resolving shared entities" );
00226         // Negotiate shared entities using vertex global IDs
00227         ParallelComm* pcomm = new ParallelComm( &moab, MPI_COMM_WORLD );
00228         rval                = pcomm->resolve_shared_ents( 0, hexes, 3, 0 );
00229         if( MB_SUCCESS != rval )
00230         {
00231             std::cerr << "ParallelComm::resolve_shared_ents failed" << std::endl;
00232             return rval;
00233         }
00234     }
00235     res_time = MPI_Wtime() - res_time;
00236 
00237     TPRINT( "Beginning parallel write" );
00238     double write_time = MPI_Wtime();
00239     // Do parallel write
00240     clock_t t = clock();
00241     std::ostringstream opts;
00242     opts << "PARALLEL=WRITE_PART";
00243     if( debug_level > 0 ) opts << ";DEBUG_IO=" << debug_level;
00244     rval = moab.write_file( output_file_name, "MOAB", opts.str().c_str() );
00245     t    = clock() - t;
00246     if( MB_SUCCESS != rval )
00247     {
00248         std::string msg;
00249         moab.get_last_error( msg );
00250         std::cerr << "File creation failed with error code: " << moab.get_error_string( rval ) << std::endl;
00251         std::cerr << "\t\"" << msg << '"' << std::endl;
00252         return (int)rval;
00253     }
00254     write_time = MPI_Wtime() - write_time;
00255 
00256     double times[3] = { gen_time, res_time, write_time };
00257     double max[3]   = { 0, 0, 0 };
00258     MPI_Reduce( times, max, 3, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
00259 
00260     // Clean up and summarize
00261     if( 0 == rank )
00262     {
00263         double sec = (double)t / CLOCKS_PER_SEC;
00264         std::cout << "Wrote " << hexes.size() * size << " hexes in " << sec << " seconds." << std::endl;
00265 
00266         if( !keep_output_file )
00267         {
00268             TPRINT( "Removing written file" );
00269             remove( output_file_name );
00270         }
00271 
00272         std::cout << "Wall time: generate: " << max[0] << ", resovle shared: " << max[1] << ", write_file: " << max[2]
00273                   << std::endl;
00274     }
00275 
00276     TPRINT( "Finalizing MPI" );
00277     return MPI_Finalize();
00278 }
00279 
00280 #define IDX( i, j, k ) ( ( num_interval + 1 ) * ( ( num_interval + 1 ) * ( k ) + ( j ) ) + ( i ) )
00281 
00282 ErrorCode generate_mesh( Interface& moab, int num_interval )
00283 {
00284     int rank, size;
00285     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00286     MPI_Comm_size( MPI_COMM_WORLD, &size );
00287 
00288     ErrorCode rval;
00289     Tag global_id = moab.globalId_tag();
00290 
00291     // Each processor will own one cube of mesh within
00292     // an 3D grid of cubes.  Calculate the dimensions of
00293     // that grid in numbers of cubes.
00294     int root = 1;
00295     while( root * root * root < size )
00296         ++root;
00297     int num_x_blocks = root;
00298     int num_y_blocks = root - 1;
00299     int num_z_blocks = root - 1;
00300     if( num_x_blocks * num_y_blocks * num_z_blocks < size ) ++num_y_blocks;
00301     if( num_x_blocks * num_y_blocks * num_z_blocks < size ) ++num_z_blocks;
00302 
00303     // calculate position of this processor in grid
00304     int my_z_block = rank / ( num_x_blocks * num_y_blocks );
00305     int rem        = rank % ( num_x_blocks * num_y_blocks );
00306     int my_y_block = rem / num_x_blocks;
00307     int my_x_block = rem % num_x_blocks;
00308 
00309     // Each processor's cube of mesh will be num_iterval^3 elements
00310     // and will be 1.0 units on a side
00311 
00312     // create vertices
00313     const int num_x_vtx = num_interval * num_x_blocks + 1;
00314     const int num_y_vtx = num_interval * num_y_blocks + 1;
00315     const int x_offset  = my_x_block * num_interval;
00316     const int y_offset  = my_y_block * num_interval;
00317     const int z_offset  = my_z_block * num_interval;
00318     double step         = 1.0 / num_interval;
00319     std::vector< EntityHandle > vertices( ( num_interval + 1 ) * ( num_interval + 1 ) * ( num_interval + 1 ) );
00320     std::vector< EntityHandle >::iterator v = vertices.begin();
00321     for( int k = 0; k <= num_interval; ++k )
00322     {
00323         for( int j = 0; j <= num_interval; ++j )
00324         {
00325             for( int i = 0; i <= num_interval; ++i )
00326             {
00327                 double coords[] = { my_x_block + i * step, my_y_block + j * step, my_z_block + k * step };
00328                 EntityHandle h;
00329                 rval = moab.create_vertex( coords, h );
00330                 if( MB_SUCCESS != rval ) return rval;
00331 
00332                 int id = 1 + x_offset + i + ( y_offset + j ) * num_x_vtx + ( z_offset + k ) * num_x_vtx * num_y_vtx;
00333                 rval   = moab.tag_set_data( global_id, &h, 1, &id );
00334                 if( MB_SUCCESS != rval ) return rval;
00335 
00336                 assert( v != vertices.end() );
00337                 *v++ = h;
00338             }
00339         }
00340     }
00341 
00342     // create hexes
00343     for( int k = 0; k < num_interval; ++k )
00344     {
00345         for( int j = 0; j < num_interval; ++j )
00346         {
00347             for( int i = 0; i < num_interval; ++i )
00348             {
00349                 assert( IDX( i + 1, j + 1, k + 1 ) < (int)vertices.size() );
00350                 const EntityHandle conn[] = { vertices[IDX( i, j, k )],
00351                                               vertices[IDX( i + 1, j, k )],
00352                                               vertices[IDX( i + 1, j + 1, k )],
00353                                               vertices[IDX( i, j + 1, k )],
00354                                               vertices[IDX( i, j, k + 1 )],
00355                                               vertices[IDX( i + 1, j, k + 1 )],
00356                                               vertices[IDX( i + 1, j + 1, k + 1 )],
00357                                               vertices[IDX( i, j + 1, k + 1 )] };
00358                 EntityHandle elem;
00359                 rval = moab.create_element( MBHEX, conn, 8, elem );
00360                 if( MB_SUCCESS != rval ) return rval;
00361             }
00362         }
00363     }
00364     /*
00365       // create interface quads
00366     for (int j = 0; j < num_interval; ++j) {
00367       for (int i = 0; i < num_interval; ++i) {
00368         EntityHandle h;
00369 
00370         const EntityHandle conn1[] = { vertices[IDX(i,  j,  0)],
00371                                          vertices[IDX(i+1,j,  0)],
00372                                          vertices[IDX(i+1,j+1,0)],
00373                                          vertices[IDX(i,  j+1,0)] };
00374         rval = moab.create_element( MBQUAD, conn1, 4, h );
00375         if (MB_SUCCESS != rval)
00376           return rval;
00377 
00378         const EntityHandle conn2[] = { vertices[IDX(i,  j,  num_interval)],
00379                                          vertices[IDX(i+1,j,  num_interval)],
00380                                          vertices[IDX(i+1,j+1,num_interval)],
00381                                          vertices[IDX(i,  j+1,num_interval)] };
00382         rval = moab.create_element( MBQUAD, conn2, 4, h );
00383         if (MB_SUCCESS != rval)
00384           return rval;
00385       }
00386     }
00387     for (int k = 0; k < num_interval; ++k) {
00388       for (int i = 0; i < num_interval; ++i) {
00389         EntityHandle h;
00390 
00391         const EntityHandle conn1[] = { vertices[IDX(i,  0,k  )],
00392                                          vertices[IDX(i+1,0,k  )],
00393                                          vertices[IDX(i+1,0,k+1)],
00394                                          vertices[IDX(i,  0,k+1)] };
00395         rval = moab.create_element( MBQUAD, conn1, 4, h );
00396         if (MB_SUCCESS != rval)
00397           return rval;
00398 
00399         const EntityHandle conn2[] = { vertices[IDX(i,  num_interval,k  )],
00400                                          vertices[IDX(i+1,num_interval,k  )],
00401                                          vertices[IDX(i+1,num_interval,k+1)],
00402                                          vertices[IDX(i,  num_interval,k+1)] };
00403         rval = moab.create_element( MBQUAD, conn2, 4, h );
00404         if (MB_SUCCESS != rval)
00405           return rval;
00406       }
00407     }
00408     for (int k = 0; k < num_interval; ++k) {
00409       for (int j = 0; j < num_interval; ++j) {
00410         EntityHandle h;
00411 
00412         const EntityHandle conn1[] = { vertices[IDX(0,j,  k  )],
00413                                          vertices[IDX(0,j+1,k  )],
00414                                          vertices[IDX(0,j+1,k+1)],
00415                                          vertices[IDX(0,j,  k+1)] };
00416         rval = moab.create_element( MBQUAD, conn1, 4, h );
00417         if (MB_SUCCESS != rval)
00418           return rval;
00419 
00420         const EntityHandle conn2[] = { vertices[IDX(num_interval,j,  k  )],
00421                                          vertices[IDX(num_interval,j+1,k  )],
00422                                          vertices[IDX(num_interval,j+1,k+1)],
00423                                          vertices[IDX(num_interval,j,  k+1)] };
00424         rval = moab.create_element( MBQUAD, conn2, 4, h );
00425         if (MB_SUCCESS != rval)
00426           return rval;
00427       }
00428     }
00429     */
00430     return MB_SUCCESS;
00431 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines