perftool.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation. Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */
 
// MOAB performance tests building mapped mesh with nodes and
// hexes created one at a time. This also creates the node to hex adjacencies.
 
#include <cmath>
#include <iostream>
#include <iomanip>
#include <ctime>
#include <cassert>
#include <list>
#include "moab/Core.hpp"
#include "moab/Skinner.hpp"
#include "moab/ReadUtilIface.hpp"
 
using namespace moab;
 
double LENGTH = 1.0;
const int DEFAULT_INTERVALS = 50;
 
void create_regular_mesh( int interval, int dimension );
void skin_common( int interval, int dim, int blocks, bool use_adj );
void skin( int intervals, int dim, int num )
{
 std::cout << "Skinning w/out adjacencies:" << std::endl;
 skin_common( intervals, dim, num, false );
}
void skin_adj( int intervals, int dim, int num )
{
 std::cout << "Skinning with adjacencies:" << std::endl;
 skin_common( intervals, dim, num, true );
}
 
void tag_time( TagType storage, bool direct, int intervals, int dim, int blocks );
 
void dense_tag( int intervals, int dim, int blocks )
{
 std::cout << "Dense Tag Time:";
 tag_time( MB_TAG_DENSE, false, intervals, dim, blocks );
}
void sparse_tag( int intervals, int dim, int blocks )
{
 std::cout << "Sparse Tag Time:";
 tag_time( MB_TAG_SPARSE, false, intervals, dim, blocks );
}
void direct_tag( int intervals, int dim, int blocks )
{
 std::cout << "Direct Tag Time:";
 tag_time( MB_TAG_DENSE, true, intervals, dim, blocks );
}
 
typedef void ( *test_func_t )( int, int, int );
const struct
{
 std::string testName;
 test_func_t testFunc;
 std::string testDesc;
} TestList[] = {
 { "skin", &skin, "Test time to get skin mesh w/out adjacencies" },
 { "skin_adj", &skin_adj, "Test time to get skin mesh with adjacencies" },
 { "sparse", &sparse_tag, "Sparse tag data manipulation" },
 { "dense", &dense_tag, "Dense tag data manipulation" },
 { "direct", &direct_tag, "Dense tag data manipulation using direct data access" },
};
const int TestListSize = sizeof( TestList ) / sizeof( TestList[0] );
 
void usage( const char* argv0, bool error = true )
{
 std::ostream& str = error ? std::cerr : std::cout;
 str << "Usage: " << argv0
 << " [-i <ints_per_side>] [-d <dimension>] [-n <test_specifc_int>] <test_name> "
 "[<test_name2> ...]"
 << std::endl;
 str << " " << argv0 << " [-h|-l]" << std::endl;
 if( error ) return;
 str << " -i : specify interverals per side (num hex = ints^3, default: " << DEFAULT_INTERVALS << std::endl;
 str << " -d : specify element dimension, default: 3" << std::endl;
 str << " -n : specify an integer value that for which the meaning is test-specific" << std::endl;
 str << " -h : print this help text." << std::endl;
 str << " -l : list available tests" << std::endl;
}
 
void list_tests()
{
 unsigned max_test_name = 0, max_test_desc = 0;
 for( int i = 0; i < TestListSize; ++i )
 {
 if( TestList[i].testName.size() > max_test_name ) max_test_name = TestList[i].testName.size();
 if( TestList[i].testDesc.size() > max_test_desc ) max_test_desc = TestList[i].testDesc.size();
 }
 std::cout << std::setw( max_test_name ) << "NAME"
 << " " << std::setw( max_test_desc ) << std::left << "DESCRIPTION" << std::endl
 << std::setfill( '-' ) << std::setw( max_test_name ) << ""
 << " " << std::setfill( '-' ) << std::setw( max_test_desc ) << "" << std::setfill( ' ' ) << std::endl;
 for( int i = 0; i < TestListSize; ++i )
 std::cout << std::setw( max_test_name ) << TestList[i].testName << " : " << std::setw( max_test_desc )
 << std::left << TestList[i].testDesc << std::endl;
}
 
int main( int argc, char* argv[] )
{
 int intervals = DEFAULT_INTERVALS;
 int dimension = 3;
 int number = 0;
 std::vector< test_func_t > test_list;
 std::list< int* > expected_list;
 bool did_help = false;
 for( int i = 1; i < argc; ++i )
 {
 if( !expected_list.empty() )
 {
 int* ptr = expected_list.front();
 expected_list.pop_front();
 char* endptr;
 *ptr = strtol( argv[i], &endptr, 0 );
 if( !endptr )
 {
 usage( argv[0] );
 std::cerr << "Expected integer value, got \"" << argv[i] << '"' << std::endl;
 return 1;
 }
 }
 else if( *argv[i] == '-' )
 { // flag
 for( int j = 1; argv[i][j]; ++j )
 {
 switch( argv[i][j] )
 {
 case 'i':
 expected_list.push_back( &intervals );
 break;
 case 'd':
 expected_list.push_back( &dimension );
 break;
 case 'n':
 expected_list.push_back( &number );
 break;
 case 'h':
 did_help = true;
 usage( argv[0], false );
 break;
 case 'l':
 did_help = true;
 list_tests();
 break;
 default:
 usage( argv[0] );
 std::cerr << "Invalid option: -" << argv[i][j] << std::endl;
 return 1;
 }
 }
 }
 else
 {
 int j = -1;
 do
 {
 ++j;
 if( j >= TestListSize )
 {
 usage( argv[0] );
 std::cerr << "Invalid test name: " << argv[i] << std::endl;
 return 1;
 }
 } while( TestList[j].testName != argv[i] );
 test_list.push_back( TestList[j].testFunc );
 }
 }
 
 if( !expected_list.empty() )
 {
 usage( argv[0] );
 std::cerr << "Missing final argument" << std::endl;
 return 1;
 }
 
 // error if no input
 if( test_list.empty() && !did_help )
 {
 usage( argv[0] );
 std::cerr << "No tests specified" << std::endl;
 return 1;
 }
 
 if( intervals < 1 )
 {
 std::cerr << "Invalid interval count: " << intervals << std::endl;
 return 1;
 }
 
 if( dimension < 1 || dimension > 3 )
 {
 std::cerr << "Invalid dimension: " << dimension << std::endl;
 return 1;
 }
 
 // now run the tests
 for( std::vector< test_func_t >::iterator i = test_list.begin(); i != test_list.end(); ++i )
 {
 test_func_t fptr = *i;
 fptr( intervals, dimension, number );
 }
 
 return 0;
}
 
void create_regular_mesh( Interface* gMB, int interval, int dim )
{
 if( dim < 1 || dim > 3 || interval < 1 )
 {
 std::cerr << "Invalid arguments" << std::endl;
 exit( 1 );
 }
 
 const int nvi = interval + 1;
 const int dims[3] = { nvi, dim > 1 ? nvi : 1, dim > 2 ? nvi : 1 };
 int num_vert = dims[0] * dims[1] * dims[2];
 
 ReadUtilIface* readMeshIface;
 gMB->query_interface( readMeshIface );
 
 EntityHandle vstart;
 std::vector< double* > arrays;
 ErrorCode rval = readMeshIface->get_node_coords( 3, num_vert, 1, vstart, arrays );
 if( MB_SUCCESS != rval || arrays.size() < 3 )
 {
 std::cerr << "Vertex creation failed" << std::endl;
 exit( 2 );
 }
 double *x = arrays[0], *y = arrays[1], *z = arrays[2];
 
 // Calculate vertex coordinates
 for( int k = 0; k < dims[2]; ++k )
 for( int j = 0; j < dims[1]; ++j )
 for( int i = 0; i < dims[0]; ++i )
 {
 *x = i;
 ++x;
 *y = j;
 ++y;
 *z = k;
 ++z;
 }
 
 const long vert_per_elem = 1 << dim; // 2^dim
 const long intervals[3] = { interval, dim > 1 ? interval : 1, dim > 2 ? interval : 1 };
 const long num_elem = intervals[0] * intervals[1] * intervals[2];
 const EntityType type = ( dim == 1 ) ? MBEDGE : ( dim == 2 ) ? MBQUAD : MBHEX;
 
 EntityHandle estart, *conn = 0;
 rval = readMeshIface->get_element_connect( num_elem, vert_per_elem, type, 0, estart, conn );
 if( MB_SUCCESS != rval || !conn )
 {
 std::cerr << "Element creation failed" << std::endl;
 exit( 2 );
 }
 
 // Offsets of element vertices in grid relative to corner closest to origin
 long c = dims[0] * dims[1];
 const long corners[8] = { 0, 1, 1 + dims[0], dims[0], c, c + 1, c + 1 + dims[0], c + dims[0] };
 
 // Populate element list
 EntityHandle* iter = conn;
 for( long z1 = 0; z1 < intervals[2]; ++z1 )
 for( long y1 = 0; y1 < intervals[1]; ++y1 )
 for( long x1 = 0; x1 < intervals[0]; ++x1 )
 {
 const long index = x1 + y1 * dims[0] + z1 * ( dims[0] * dims[1] );
 for( long j = 0; j < vert_per_elem; ++j, ++iter )
 *iter = index + corners[j] + vstart;
 }
 
 // notify MOAB of the new elements
 rval = readMeshIface->update_adjacencies( estart, num_elem, vert_per_elem, conn );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Element update failed" << std::endl;
 exit( 2 );
 }
}
 
void skin_common( int interval, int dim, int num, bool use_adj )
{
 Core moab;
 Interface* gMB = &moab;
 ErrorCode rval;
 double d;
 clock_t t, tt;
 
 create_regular_mesh( gMB, interval, dim );
 
 Range skin, verts, elems;
 rval = gMB->get_entities_by_dimension( 0, dim, elems );
 assert( MB_SUCCESS == rval );
 assert( !elems.empty() );
 
 Skinner tool( gMB );
 
 t = clock();
 rval = tool.find_skin( 0, elems, true, verts, 0, use_adj, false );
 t = clock() - t;
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Search for skin vertices failed" << std::endl;
 exit( 2 );
 }
 d = ( (double)t ) / CLOCKS_PER_SEC;
 std::cout << "Got " << verts.size() << " skin vertices in " << d << " seconds." << std::endl;
 
 t = 0;
 if( num < 1 ) num = 1000;
 long blocksize = elems.size() / num;
 if( !blocksize ) blocksize = 1;
 long numblocks = elems.size() / blocksize;
 Range::iterator it = elems.begin();
 for( long i = 0; i < numblocks; ++i )
 {
 verts.clear();
 Range::iterator end = it + blocksize;
 Range blockelems;
 blockelems.merge( it, end );
 it = end;
 tt = clock();
 rval = tool.find_skin( 0, blockelems, true, verts, 0, use_adj, false );
 t += clock() - tt;
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Search for skin vertices failed" << std::endl;
 exit( 2 );
 }
 }
 d = ( (double)t ) / CLOCKS_PER_SEC;
 std::cout << "Got skin vertices for " << numblocks << " blocks of " << blocksize << " elements in " << d
 << " seconds." << std::endl;
 
 for( int e = 0; e < 2; ++e )
 { // do this twice
 if( e == 1 )
 {
 // create all interior faces
 skin.clear();
 t = clock();
 gMB->get_adjacencies( elems, dim - 1, true, skin, Interface::UNION );
 t = clock() - t;
 d = ( (double)t ) / CLOCKS_PER_SEC;
 std::cout << "Created " << skin.size() << " entities of dimension-1 in " << d << " seconds" << std::endl;
 }
 
 skin.clear();
 t = clock();
 rval = tool.find_skin( 0, elems, false, skin, 0, use_adj, true );
 t = clock() - t;
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Search for skin vertices failed" << std::endl;
 exit( 2 );
 }
 d = ( (double)t ) / CLOCKS_PER_SEC;
 std::cout << "Got " << skin.size() << " skin elements in " << d << " seconds." << std::endl;
 
 t = 0;
 it = elems.begin();
 for( long i = 0; i < numblocks; ++i )
 {
 skin.clear();
 Range::iterator end = it + blocksize;
 Range blockelems;
 blockelems.merge( it, end );
 it = end;
 tt = clock();
 rval = tool.find_skin( 0, blockelems, false, skin, 0, use_adj, true );
 t += clock() - tt;
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Search for skin elements failed" << std::endl;
 exit( 2 );
 }
 }
 d = ( (double)t ) / CLOCKS_PER_SEC;
 std::cout << "Got skin elements for " << numblocks << " blocks of " << blocksize << " elements in " << d
 << " seconds." << std::endl;
 }
}
 
void tag_time( TagType storage, bool direct, int intervals, int dim, int blocks )
{
 Core moab;
 Interface& mb = moab;
 create_regular_mesh( &mb, intervals, dim );
 
 // Create tag in which to store data
 Tag tag;
 mb.tag_get_handle( "data", 1, MB_TYPE_DOUBLE, tag, storage | MB_TAG_CREAT );
 
 // Make up some arbitrary iterative calculation for timing purposes:
 // set each value v_n = (V + v_n)/2 until all values are within
 // epsilon of V.
 std::vector< double > data;
 Range verts;
 mb.get_entities_by_type( 0, MBVERTEX, verts );
 
 clock_t t = clock();
 
 // initialize
 if( direct )
 {
 Range::iterator i, j = verts.begin();
 void* ptr;
 int count;
 while( j != verts.end() )
 {
 mb.tag_iterate( tag, i, verts.end(), count, ptr );
 double* arr = reinterpret_cast< double* >( ptr );
 for( j = i + count; i != j; ++i, ++arr )
 *arr = ( 11.0 * *i + 7.0 ) / ( *i );
 }
 }
 else
 {
 data.resize( verts.size() );
 double* arr = &data[0];
 for( Range::iterator i = verts.begin(); i != verts.end(); ++i, ++arr )
 *arr = ( 11.0 * *i + 7.0 ) / ( *i );
 mb.tag_set_data( tag, verts, &data[0] );
 }
 
 // iterate
 const double v0 = acos( -1.0 ); // pi
 size_t iter_count = 0;
 double max_diff;
 const size_t num_verts = verts.size();
 do
 {
 if( direct )
 {
 max_diff = 0.0;
 Range::iterator i, j = verts.begin();
 void* ptr;
 while( j != verts.end() )
 {
 int count;
 mb.tag_iterate( tag, i, verts.end(), count, ptr );
 double* arr = reinterpret_cast< double* >( ptr );
 
 for( j = i + count; i != j; ++i, ++arr )
 {
 *arr = 0.5 * ( *arr + v0 );
 double diff = fabs( *arr - v0 );
 if( diff > max_diff ) max_diff = diff;
 }
 }
 }
 else if( blocks < 1 )
 {
 max_diff = 0.0;
 mb.tag_get_data( tag, verts, &data[0] );
 for( size_t i = 0; i < data.size(); ++i )
 {
 data[i] = 0.5 * ( v0 + data[i] );
 double diff = fabs( data[i] - v0 );
 if( diff > max_diff ) max_diff = diff;
 }
 mb.tag_set_data( tag, verts, &data[0] );
 }
 else
 {
 max_diff = 0.0;
 Range r;
 Range::iterator it = verts.begin();
 size_t step = num_verts / blocks;
 for( int j = 0; j < blocks - 1; ++j )
 {
 Range::iterator nx = it;
 nx += step;
 r.clear();
 r.merge( it, nx );
 mb.tag_get_data( tag, r, &data[0] );
 it = nx;
 
 for( size_t i = 0; i < step; ++i )
 {
 data[i] = 0.5 * ( v0 + data[i] );
 double diff = fabs( data[i] - v0 );
 if( diff > max_diff ) max_diff = diff;
 }
 mb.tag_set_data( tag, r, &data[0] );
 }
 
 r.clear();
 r.merge( it, verts.end() );
 mb.tag_get_data( tag, r, &data[0] );
 for( size_t i = 0; i < ( num_verts - ( blocks - 1 ) * step ); ++i )
 {
 data[i] = 0.5 * ( v0 + data[i] );
 double diff = fabs( data[i] - v0 );
 if( diff > max_diff ) max_diff = diff;
 }
 mb.tag_set_data( tag, r, &data[0] );
 }
 ++iter_count;
 // std::cout << iter_count << " " << max_diff << std::endl;
 } while( max_diff > 1e-6 );
 
 double secs = ( clock() - t ) / (double)CLOCKS_PER_SEC;
 std::cout << " " << iter_count << " iterations in " << secs << " seconds" << std::endl;
}