parallel_write_test.cpp

Go to the documentation of this file.

#include "moab/Core.hpp"
#include "moab/ParallelComm.hpp"
#include "MBTagConventions.hpp"
#include "moab_mpi.h"
#include <cstdlib>
#include <iostream>
#include <ctime>
#include <cstring>
#include <cmath>
#include <cassert>
#include <cstdio>
#include <sstream>
 
using namespace moab;
 
#define TPRINT( A ) tprint( ( A ) )
static void tprint( const char* A )
{
 int rank;
 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
 char buffer[128];
 sprintf( buffer, "%02d: %6.2f: %s\n", rank, (double)clock() / CLOCKS_PER_SEC, A );
 fputs( buffer, stderr );
}
 
const int DEFAULT_INTERVALS = 2;
const char* DEFAULT_FILE_NAME = "parallel_write_test.h5m";
 
// Create mesh for each processor that is a cube of hexes
// with the specified interval count along each edge. Cubes
// of mesh will be positioned and vertex IDs assigned such that
// there are shared entities. If the cubic root of the number
// of processors is a whole number, then each processors mesh
// will be a cube in a grid with that number of processor blocks
// along each edge. Otherwise processor blocks will be arranged
// within the subset of the grid that is the ceiling of the cubic
// root of the comm size such that there are no disjoint regions.
ErrorCode generate_mesh( Interface& moab, int intervals );
 
const char args[] = "[-i <intervals>] [-o <filename>] [-L <filename>] [-g <n>]";
void help()
{
 std::cout << "parallel_write_test " << args << std::endl
 << " -i <N> Each processor owns an NxNxN cube of hex elements (default: " << DEFAULT_INTERVALS << ")"
 << std::endl
 << " -o <name> Retain output file and name it as specified." << std::endl
 << " -L <name> Write local mesh to file name prefixed with MPI rank" << std::endl
 << " -g <n> Specify writer debug output level" << std::endl
 << " -R Skip resolve of shared entities (interface ents will be duplicated in file)" << std::endl
 << std::endl
 << "This program creates a (non-strict) subset of a regular hex mesh "
 "such that the mesh is already partitioned, and then attempts to "
 "write that mesh using MOAB's parallel HDF5 writer. The mesh size "
 "will scale with the number of processors and the number of elements "
 "per processor (the latter is a function of the value specified "
 "with the '-i' flag.)"
 << std::endl
 << std::endl
 << "Let N = ceil(cbrt(P)), where P is the number of processes. "
 "The mesh will be some subset of a cube with one corner at the "
 "origin and the other at (N,N,N). Each processor will own a "
 "non-overlapping 1x1x1 unit block of mesh within that cube. "
 "If P is a power of 3, then the entire NxNxN cube will be "
 "filled with hex elements. Otherwise, some connected subset "
 "of the cube will be meshed. Each processor is assigned a "
 "sub-block of the cube by rank where the blocks are enumerated "
 "sequentally with x increasing most rapidly and z least rapidly."
 << std::endl
 << std::endl
 << "The size of the mesh owned by each processor is controlled by "
 "the number of intervals along each edge of its block of mesh. "
 "If each block has N intervals, than each processor will have "
 "N^3 hex elements."
 << std::endl
 << std::endl;
}
 
int main( int argc, char* argv[] )
{
 int ierr = MPI_Init( &argc, &argv );
 if( ierr )
 {
 std::cerr << "MPI_Init failed with error code: " << ierr << std::endl;
 return ierr;
 }
 int rank;
 ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
 if( ierr )
 {
 std::cerr << "MPI_Comm_rank failed with error code: " << ierr << std::endl;
 return ierr;
 }
 int size;
 ierr = MPI_Comm_size( MPI_COMM_WORLD, &size );
 if( ierr )
 {
 std::cerr << "MPI_Comm_size failed with error code: " << ierr << std::endl;
 return ierr;
 }
 
 // settings controlled by CL flags
 const char* output_file_name = 0;
 const char* indiv_file_name = 0;
 int intervals = 0, debug_level = 0;
 // state for CL flag processing
 bool expect_intervals = false;
 bool expect_file_name = false;
 bool expect_indiv_file = false;
 bool skip_resolve_shared = false;
 bool expect_debug_level = false;
 // process CL args
 for( int i = 1; i < argc; ++i )
 {
 if( expect_intervals )
 {
 char* endptr = 0;
 intervals = (int)strtol( argv[i], &endptr, 0 );
 if( *endptr || intervals < 1 )
 {
 std::cerr << "Invalid block interval value: " << argv[i] << std::endl;
 return 1;
 }
 expect_intervals = false;
 }
 else if( expect_indiv_file )
 {
 indiv_file_name = argv[i];
 expect_indiv_file = false;
 }
 else if( expect_file_name )
 {
 output_file_name = argv[i];
 expect_file_name = false;
 }
 else if( expect_debug_level )
 {
 debug_level = atoi( argv[i] );
 if( debug_level < 1 )
 {
 std::cerr << "Invalid argument following -g flag: \"" << argv[i] << '"' << std::endl;
 return 1;
 }
 expect_debug_level = false;
 }
 else if( !strcmp( "-i", argv[i] ) )
 expect_intervals = true;
 else if( !strcmp( "-o", argv[i] ) )
 expect_file_name = true;
 else if( !strcmp( "-L", argv[i] ) )
 expect_indiv_file = true;
 else if( !strcmp( "-R", argv[i] ) )
 skip_resolve_shared = true;
 else if( !strcmp( "-g", argv[i] ) )
 expect_debug_level = true;
 else if( !strcmp( "-h", argv[i] ) )
 {
 help();
 return 0;
 }
 else
 {
 std::cerr << "Unexpected argument: " << argv[i] << std::endl
 << "Usage: " << argv[0] << " " << args << std::endl
 << " " << argv[0] << " -h" << std::endl
 << " Try '-h' for help." << std::endl;
 return 1;
 }
 }
 // Check for missing argument after last CL flag
 if( expect_file_name || expect_intervals || expect_indiv_file )
 {
 std::cerr << "Missing argument for '" << argv[argc - 1] << "'" << std::endl;
 return 1;
 }
 // If intervals weren't specified, use default
 if( intervals == 0 )
 {
 std::cout << "Using default interval count: " << DEFAULT_INTERVALS << std::endl;
 intervals = DEFAULT_INTERVALS;
 }
 // If no output file was specified, use default one and note that
 // we need to delete it when the test completes.
 bool keep_output_file = true;
 if( !output_file_name )
 {
 output_file_name = DEFAULT_FILE_NAME;
 keep_output_file = false;
 }
 
 // Create mesh
 TPRINT( "Generating mesh" );
 double gen_time = MPI_Wtime();
 Core mb;
 Interface& moab = mb;
 ErrorCode rval = generate_mesh( moab, intervals );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Mesh creation failed with error code: " << rval << std::endl;
 return (int)rval;
 }
 gen_time = MPI_Wtime() - gen_time;
 
 // Write out local mesh on each processor if requested.
 if( indiv_file_name )
 {
 TPRINT( "Writing individual file" );
 char buffer[64];
 int width = (int)ceil( log10( size ) );
 sprintf( buffer, "%0*d-", width, rank );
 std::string name( buffer );
 name += indiv_file_name;
 rval = moab.write_file( name.c_str() );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Failed to write file: " << name << std::endl;
 return (int)rval;
 }
 }
 
 double res_time = MPI_Wtime();
 Range hexes;
 moab.get_entities_by_type( 0, MBHEX, hexes );
 if( !skip_resolve_shared )
 {
 TPRINT( "Resolving shared entities" );
 // Negotiate shared entities using vertex global IDs
 ParallelComm* pcomm = new ParallelComm( &moab, MPI_COMM_WORLD );
 rval = pcomm->resolve_shared_ents( 0, hexes, 3, 0 );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "ParallelComm::resolve_shared_ents failed" << std::endl;
 return rval;
 }
 }
 res_time = MPI_Wtime() - res_time;
 
 TPRINT( "Beginning parallel write" );
 double write_time = MPI_Wtime();
 // Do parallel write
 clock_t t = clock();
 std::ostringstream opts;
 opts << "PARALLEL=WRITE_PART";
 if( debug_level > 0 ) opts << ";DEBUG_IO=" << debug_level;
 rval = moab.write_file( output_file_name, "MOAB", opts.str().c_str() );
 t = clock() - t;
 if( MB_SUCCESS != rval )
 {
 std::string msg;
 moab.get_last_error( msg );
 std::cerr << "File creation failed with error code: " << moab.get_error_string( rval ) << std::endl;
 std::cerr << "\t\"" << msg << '"' << std::endl;
 return (int)rval;
 }
 write_time = MPI_Wtime() - write_time;
 
 double times[3] = { gen_time, res_time, write_time };
 double max[3] = { 0, 0, 0 };
 MPI_Reduce( times, max, 3, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
 
 // Clean up and summarize
 if( 0 == rank )
 {
 double sec = (double)t / CLOCKS_PER_SEC;
 std::cout << "Wrote " << hexes.size() * size << " hexes in " << sec << " seconds." << std::endl;
 
 if( !keep_output_file )
 {
 TPRINT( "Removing written file" );
 remove( output_file_name );
 }
 
 std::cout << "Wall time: generate: " << max[0] << ", resovle shared: " << max[1] << ", write_file: " << max[2]
 << std::endl;
 }
 
 TPRINT( "Finalizing MPI" );
 return MPI_Finalize();
}
 
#define IDX( i, j, k ) ( ( num_interval + 1 ) * ( ( num_interval + 1 ) * ( k ) + ( j ) ) + ( i ) )
 
ErrorCode generate_mesh( Interface& moab, int num_interval )
{
 int rank, size;
 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
 MPI_Comm_size( MPI_COMM_WORLD, &size );
 
 ErrorCode rval;
 Tag global_id = moab.globalId_tag();
 
 // Each processor will own one cube of mesh within
 // an 3D grid of cubes. Calculate the dimensions of
 // that grid in numbers of cubes.
 int root = 1;
 while( root * root * root < size )
 ++root;
 int num_x_blocks = root;
 int num_y_blocks = root - 1;
 int num_z_blocks = root - 1;
 if( num_x_blocks * num_y_blocks * num_z_blocks < size ) ++num_y_blocks;
 if( num_x_blocks * num_y_blocks * num_z_blocks < size ) ++num_z_blocks;
 
 // calculate position of this processor in grid
 int my_z_block = rank / ( num_x_blocks * num_y_blocks );
 int rem = rank % ( num_x_blocks * num_y_blocks );
 int my_y_block = rem / num_x_blocks;
 int my_x_block = rem % num_x_blocks;
 
 // Each processor's cube of mesh will be num_iterval^3 elements
 // and will be 1.0 units on a side
 
 // create vertices
 const int num_x_vtx = num_interval * num_x_blocks + 1;
 const int num_y_vtx = num_interval * num_y_blocks + 1;
 const int x_offset = my_x_block * num_interval;
 const int y_offset = my_y_block * num_interval;
 const int z_offset = my_z_block * num_interval;
 double step = 1.0 / num_interval;
 std::vector< EntityHandle > vertices( ( num_interval + 1 ) * ( num_interval + 1 ) * ( num_interval + 1 ) );
 std::vector< EntityHandle >::iterator v = vertices.begin();
 for( int k = 0; k <= num_interval; ++k )
 {
 for( int j = 0; j <= num_interval; ++j )
 {
 for( int i = 0; i <= num_interval; ++i )
 {
 double coords[] = { my_x_block + i * step, my_y_block + j * step, my_z_block + k * step };
 EntityHandle h;
 rval = moab.create_vertex( coords, h );
 if( MB_SUCCESS != rval ) return rval;
 
 int id = 1 + x_offset + i + ( y_offset + j ) * num_x_vtx + ( z_offset + k ) * num_x_vtx * num_y_vtx;
 rval = moab.tag_set_data( global_id, &h, 1, &id );
 if( MB_SUCCESS != rval ) return rval;
 
 assert( v != vertices.end() );
 *v++ = h;
 }
 }
 }
 
 // create hexes
 for( int k = 0; k < num_interval; ++k )
 {
 for( int j = 0; j < num_interval; ++j )
 {
 for( int i = 0; i < num_interval; ++i )
 {
 assert( IDX( i + 1, j + 1, k + 1 ) < (int)vertices.size() );
 const EntityHandle conn[] = { vertices[IDX( i, j, k )],
 vertices[IDX( i + 1, j, k )],
 vertices[IDX( i + 1, j + 1, k )],
 vertices[IDX( i, j + 1, k )],
 vertices[IDX( i, j, k + 1 )],
 vertices[IDX( i + 1, j, k + 1 )],
 vertices[IDX( i + 1, j + 1, k + 1 )],
 vertices[IDX( i, j + 1, k + 1 )] };
 EntityHandle elem;
 rval = moab.create_element( MBHEX, conn, 8, elem );
 if( MB_SUCCESS != rval ) return rval;
 }
 }
 }
 /*
 // create interface quads
 for (int j = 0; j < num_interval; ++j) {
 for (int i = 0; i < num_interval; ++i) {
 EntityHandle h;
 
 const EntityHandle conn1[] = { vertices[IDX(i, j, 0)],
 vertices[IDX(i+1,j, 0)],
 vertices[IDX(i+1,j+1,0)],
 vertices[IDX(i, j+1,0)] };
 rval = moab.create_element( MBQUAD, conn1, 4, h );
 if (MB_SUCCESS != rval)
 return rval;
 
 const EntityHandle conn2[] = { vertices[IDX(i, j, num_interval)],
 vertices[IDX(i+1,j, num_interval)],
 vertices[IDX(i+1,j+1,num_interval)],
 vertices[IDX(i, j+1,num_interval)] };
 rval = moab.create_element( MBQUAD, conn2, 4, h );
 if (MB_SUCCESS != rval)
 return rval;
 }
 }
 for (int k = 0; k < num_interval; ++k) {
 for (int i = 0; i < num_interval; ++i) {
 EntityHandle h;
 
 const EntityHandle conn1[] = { vertices[IDX(i, 0,k )],
 vertices[IDX(i+1,0,k )],
 vertices[IDX(i+1,0,k+1)],
 vertices[IDX(i, 0,k+1)] };
 rval = moab.create_element( MBQUAD, conn1, 4, h );
 if (MB_SUCCESS != rval)
 return rval;
 
 const EntityHandle conn2[] = { vertices[IDX(i, num_interval,k )],
 vertices[IDX(i+1,num_interval,k )],
 vertices[IDX(i+1,num_interval,k+1)],
 vertices[IDX(i, num_interval,k+1)] };
 rval = moab.create_element( MBQUAD, conn2, 4, h );
 if (MB_SUCCESS != rval)
 return rval;
 }
 }
 for (int k = 0; k < num_interval; ++k) {
 for (int j = 0; j < num_interval; ++j) {
 EntityHandle h;
 
 const EntityHandle conn1[] = { vertices[IDX(0,j, k )],
 vertices[IDX(0,j+1,k )],
 vertices[IDX(0,j+1,k+1)],
 vertices[IDX(0,j, k+1)] };
 rval = moab.create_element( MBQUAD, conn1, 4, h );
 if (MB_SUCCESS != rval)
 return rval;
 
 const EntityHandle conn2[] = { vertices[IDX(num_interval,j, k )],
 vertices[IDX(num_interval,j+1,k )],
 vertices[IDX(num_interval,j+1,k+1)],
 vertices[IDX(num_interval,j, k+1)] };
 rval = moab.create_element( MBQUAD, conn2, 4, h );
 if (MB_SUCCESS != rval)
 return rval;
 }
 }
 */
 return MB_SUCCESS;
}