Sierra Toolkit Version of the Day
UseCase_Rebal_2.cpp
00001 /*------------------------------------------------------------------------*/
00002 /*                 Copyright 2010 Sandia Corporation.                     */
00003 /*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
00004 /*  license for use of this work by or on behalf of the U.S. Government.  */
00005 /*  Export of this program may require a license from the                 */
00006 /*  United States Government.                                             */
00007 /*------------------------------------------------------------------------*/
00008 
00009 #include <use_cases/UseCase_Rebal_2.hpp>
00010 
00011 #include <stk_util/parallel/Parallel.hpp>
00012 #include <stk_util/parallel/ParallelReduce.hpp>
00013 
00014 #include <stk_mesh/base/Comm.hpp>
00015 #include <stk_mesh/base/FieldData.hpp>
00016 #include <stk_mesh/base/GetEntities.hpp>
00017 #include <stk_mesh/base/Types.hpp>
00018 #include <stk_mesh/fem/CreateAdjacentEntities.hpp>
00019 
00020 #include <stk_rebalance/Rebalance.hpp>
00021 #include <stk_rebalance/Partition.hpp>
00022 #include <stk_rebalance/ZoltanPartition.hpp>
00023 
00024 #include <stk_rebalance_utils/RebalanceUtils.hpp>
00025 
00026 //----------------------------------------------------------------------
00027 
00028 using namespace stk::mesh::fixtures;
00029 
00030 typedef stk::mesh::Field<double> ScalarField ;
00031 
00032 namespace stk {
00033 namespace rebalance {
00034 namespace use_cases {
00035 
00036 namespace {
00037 
00038   void sum_element_weights_through_relations( stk::mesh::EntityVector & elements, 
00039       ScalarField & field, const std::vector<stk::mesh::EntityRank> & ranks )
00040   {
00041     for( size_t i = 0; i < ranks.size(); ++i )
00042     {
00043       for( size_t ielem = 0; ielem < elements.size(); ++ielem )
00044       {
00045         stk::mesh::Entity * elem = elements[ielem];
00046         double * elem_weight = field_data( field , *elem );
00047         const stk::mesh::PairIterRelation rel = elem->relations( ranks[i] );
00048         const unsigned num_entities = rel.size();
00049 
00050         for ( unsigned j = 0 ; j < num_entities ; ++j ) 
00051         {
00052           stk::mesh::Entity & entity = * rel[j].entity();
00053           const double * entity_weight = field_data( field , entity );
00054           elem_weight[0] += entity_weight[0];
00055         }
00056       }
00057     }
00058   }
00059 
00060 } // namespace
00061 
00062 bool test_heavy_nodes( stk::ParallelMachine pm )
00063 {
00064   const size_t nx = 3;
00065   const size_t ny = 3;
00066   const size_t nz = 3;
00067 
00068   const unsigned p_size = stk::parallel_machine_size(pm);
00069   const unsigned p_rank = stk::parallel_machine_rank(pm);
00070 
00071   stk::mesh::fixtures::HexFixture fixture(pm, nx, ny, nz);
00072 
00073   stk::mesh::fem::FEMMetaData & fem_meta  = fixture.m_fem_meta;
00074   stk::mesh::BulkData & bulk  = fixture.m_bulk_data;
00075 
00076   // Put weights field on all entities
00077   const stk::mesh::EntityRank element_rank = fem_meta.element_rank();
00078   const stk::mesh::EntityRank face_rank = fem_meta.face_rank();
00079   const stk::mesh::EntityRank edge_rank = fem_meta.edge_rank();
00080   const stk::mesh::EntityRank node_rank = fem_meta.node_rank();
00081   ScalarField & weight_field( fem_meta.declare_field< ScalarField >( "entity_weights" ) );
00082   stk::mesh::put_field(weight_field , element_rank , fem_meta.universal_part() );
00083   stk::mesh::put_field(weight_field , face_rank , fem_meta.universal_part() );
00084   stk::mesh::put_field(weight_field , edge_rank , fem_meta.universal_part() );
00085   stk::mesh::put_field(weight_field , node_rank , fem_meta.universal_part() );
00086 
00087   fem_meta.commit();
00088 
00089   // Configure our mesh on proc 0
00090   std::vector<stk::mesh::EntityId> my_element_ids;
00091   if( 0 == p_rank )
00092   {
00093     for ( unsigned i = 0 ; i < nx*ny*nz; ++i )
00094       my_element_ids.push_back(i+1);
00095   }
00096 
00097   fixture.generate_mesh(my_element_ids);
00098 
00099   // create faces and edges for the mesh
00100   stk::mesh::PartVector add_parts;
00101   stk::mesh::create_adjacent_entities(bulk, add_parts);
00102 
00103   bulk.modification_begin();
00104 
00105   // Assign entity weights
00106   if( 0 == p_rank ) 
00107   {
00108     // Get the faces on the x=0 plane and give them a characteristic weight
00109     stk::mesh::EntityVector selected_nodes;
00110     stk::mesh::EntityVector selected_faces;
00111     stk::mesh::EntityVector one_face;
00112     for ( unsigned j = 0 ; j < ny; ++j ) 
00113       for ( unsigned k = 0 ; k < nz; ++k ) 
00114       {
00115         selected_nodes.clear();
00116         selected_nodes.push_back( fixture.node(0, j,   k  ) );
00117         selected_nodes.push_back( fixture.node(0, j+1, k  ) );
00118         selected_nodes.push_back( fixture.node(0, j,   k+1) );
00119         selected_nodes.push_back( fixture.node(0, j+1, k+1) );
00120         stk::mesh::get_entities_through_relations(selected_nodes, face_rank, one_face);
00121         selected_faces.push_back(one_face[0]);
00122       }
00123 
00124     for( size_t iface = 0; iface < selected_faces.size(); ++iface )
00125     {
00126       stk::mesh::Entity * face = selected_faces[iface];
00127       double * const weight = stk::mesh::field_data( weight_field, *face );
00128       weight[0] = 10.0;
00129     }
00130 
00131     // Get the edges on the boundary of the x=0 plane and give them a characteristic weight
00132     stk::mesh::EntityVector selected_edges;
00133     stk::mesh::EntityVector one_edge;
00134     for ( unsigned j = 0 ; j < ny; ++j ) 
00135     {
00136       selected_nodes.clear();
00137       selected_nodes.push_back( fixture.node(0, j,   0) );
00138       selected_nodes.push_back( fixture.node(0, j+1, 0) );
00139       stk::mesh::get_entities_through_relations(selected_nodes, edge_rank, one_edge);
00140       selected_edges.push_back(one_edge[0]);
00141       selected_nodes.clear();
00142       selected_nodes.push_back( fixture.node(0, j,   nz) );
00143       selected_nodes.push_back( fixture.node(0, j+1, nz) );
00144       stk::mesh::get_entities_through_relations(selected_nodes, edge_rank, one_edge);
00145       selected_edges.push_back(one_edge[0]);
00146     }
00147     for ( unsigned k = 0 ; k < nz; ++k ) 
00148     {
00149       selected_nodes.clear();
00150       selected_nodes.push_back( fixture.node(0, 0, k) );
00151       selected_nodes.push_back( fixture.node(0, 0, k+1) );
00152       stk::mesh::get_entities_through_relations(selected_nodes, edge_rank, one_edge);
00153       selected_edges.push_back(one_edge[0]);
00154       selected_nodes.clear();
00155       selected_nodes.push_back( fixture.node(0, ny, k) );
00156       selected_nodes.push_back( fixture.node(0, ny, k+1) );
00157       stk::mesh::get_entities_through_relations(selected_nodes, edge_rank, one_edge);
00158       selected_edges.push_back(one_edge[0]);
00159     }
00160     for( size_t iedge = 0; iedge < selected_edges.size(); ++iedge )
00161     {
00162       stk::mesh::Entity * edge = selected_edges[iedge];
00163       double * const weight = stk::mesh::field_data( weight_field, *edge );
00164       weight[0] = 100.0;
00165     }
00166 
00167     // Finally, give the corner nodes of the x=0 plane a characteristic weight
00168     selected_nodes.clear();
00169     double * weight = stk::mesh::field_data( weight_field, *fixture.node(0, 0, 0) );
00170     weight[0] = 1000.0;
00171     weight = stk::mesh::field_data( weight_field, *fixture.node(0, ny, 0) );
00172     weight[0] = 1000.0;
00173     weight = stk::mesh::field_data( weight_field, *fixture.node(0, 0, nz) );
00174     weight[0] = 1000.0;
00175     weight = stk::mesh::field_data( weight_field, *fixture.node(0, ny, nz) );
00176     weight[0] = 1000.0;
00177 
00178     // Assign element weights
00179     for( size_t i = 0; i < my_element_ids.size(); ++i )
00180     {
00181       stk::mesh::Entity * elem = bulk.get_entity(element_rank, my_element_ids[i]);
00182       double * const e_weight = stk::mesh::field_data( weight_field , *elem );
00183       *e_weight = 1.0;
00184     }
00185     //
00186     // Get the elements on the x=0 plane and sum in weights from relations
00187     selected_nodes.clear();
00188     for ( unsigned j = 0 ; j < ny+1; ++j )
00189       for ( unsigned k = 0 ; k < nz+1; ++k ) 
00190         selected_nodes.push_back( fixture.node(0, j, k) );
00191 
00192     std::vector<stk::mesh::EntityRank> ranks;
00193     ranks.push_back(face_rank);
00194     ranks.push_back(edge_rank);
00195     ranks.push_back(node_rank);
00196     stk::mesh::EntityVector selected_elems;
00197     for ( unsigned j = 0 ; j < ny; ++j ) 
00198       for ( unsigned k = 0 ; k < nz; ++k ) 
00199       {
00200         selected_nodes.clear();
00201         selected_nodes.push_back( fixture.node(0, j,   k  ) );
00202         selected_nodes.push_back( fixture.node(0, j+1, k  ) );
00203         selected_nodes.push_back( fixture.node(0, j,   k+1) );
00204         selected_nodes.push_back( fixture.node(0, j+1, k+1) );
00205         stk::mesh::get_entities_through_relations(selected_nodes, element_rank, one_face);
00206         selected_elems.push_back(one_face[0]);
00207       }
00208     sum_element_weights_through_relations(selected_elems, weight_field, ranks);
00209   }
00210 
00211   bulk.modification_end();
00212 
00213   // Use Zoltan to determine new partition
00214   Teuchos::ParameterList emptyList;
00215   stk::rebalance::Zoltan zoltan_partition(pm, fixture.m_spatial_dimension, emptyList);
00216 
00217   stk::rebalance::rebalance(bulk, fem_meta.universal_part(), &fixture.m_coord_field, &weight_field, zoltan_partition);
00218 
00219   const double imbalance_threshold = ( 3 == p_size )? 1.45 // Zoltan does not do so well for 3 procs
00220                                                     : 1.1; // ... but does pretty well for 2 and 4 procs
00221   const bool do_rebal = imbalance_threshold < stk::rebalance::check_balance(bulk, &weight_field, element_rank);
00222 
00223   if( 0 == p_rank )
00224     std::cerr << std::endl 
00225       << "imbalance_threshold after rebalance = " << imbalance_threshold << ", " << do_rebal << std::endl;
00226 
00227 
00228   // Check that we satisfy our threshhold
00229   bool result = !do_rebal;
00230 
00231   // And verify that all dependent entities are on the same proc as their parent element
00232   {
00233     stk::mesh::EntityVector entities;
00234     stk::mesh::Selector selector = fem_meta.locally_owned_part();
00235 
00236     get_selected_entities(selector, bulk.buckets(node_rank), entities);
00237     result &= verify_dependent_ownership(element_rank, entities);
00238 
00239     get_selected_entities(selector, bulk.buckets(edge_rank), entities);
00240     result &= verify_dependent_ownership(element_rank, entities);
00241 
00242     get_selected_entities(selector, bulk.buckets(face_rank), entities);
00243     result &= verify_dependent_ownership(element_rank, entities);
00244   }
00245 
00246   return result;
00247 }
00248 
00249 } //namespace use_cases
00250 } //namespace rebalance
00251 } //namespace stk
00252 
00253 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends