Sierra Toolkit  Version of the Day
UseCase_Rebal_2.cpp
1 /*------------------------------------------------------------------------*/
2 /* Copyright 2010 Sandia Corporation. */
3 /* Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive */
4 /* license for use of this work by or on behalf of the U.S. Government. */
5 /* Export of this program may require a license from the */
6 /* United States Government. */
7 /*------------------------------------------------------------------------*/
8 
9 #include <use_cases/UseCase_Rebal_2.hpp>
10 
11 #include <stk_util/parallel/Parallel.hpp>
12 #include <stk_util/parallel/ParallelReduce.hpp>
13 
14 #include <stk_mesh/base/Comm.hpp>
15 #include <stk_mesh/base/FieldData.hpp>
16 #include <stk_mesh/base/GetEntities.hpp>
17 #include <stk_mesh/base/Types.hpp>
18 #include <stk_mesh/fem/CreateAdjacentEntities.hpp>
19 
23 
24 #include <stk_rebalance_utils/RebalanceUtils.hpp>
25 
26 //----------------------------------------------------------------------
27 
28 using namespace stk_classic::mesh::fixtures;
29 
31 
32 namespace stk_classic {
33 namespace rebalance {
34 namespace use_cases {
35 
36 namespace {
37 
38  void sum_element_weights_through_relations( stk_classic::mesh::EntityVector & elements,
39  ScalarField & field, const std::vector<stk_classic::mesh::EntityRank> & ranks )
40  {
41  for( size_t i = 0; i < ranks.size(); ++i )
42  {
43  for( size_t ielem = 0; ielem < elements.size(); ++ielem )
44  {
45  stk_classic::mesh::Entity * elem = elements[ielem];
46  double * elem_weight = field_data( field , *elem );
47  const stk_classic::mesh::PairIterRelation rel = elem->relations( ranks[i] );
48  const unsigned num_entities = rel.size();
49 
50  for ( unsigned j = 0 ; j < num_entities ; ++j )
51  {
52  stk_classic::mesh::Entity & entity = * rel[j].entity();
53  const double * entity_weight = field_data( field , entity );
54  elem_weight[0] += entity_weight[0];
55  }
56  }
57  }
58  }
59 
60 } // namespace
61 
62 bool test_heavy_nodes( stk_classic::ParallelMachine pm )
63 {
64  const size_t nx = 3;
65  const size_t ny = 3;
66  const size_t nz = 3;
67 
68  const unsigned p_size = stk_classic::parallel_machine_size(pm);
69  const unsigned p_rank = stk_classic::parallel_machine_rank(pm);
70 
71  stk_classic::mesh::fixtures::HexFixture fixture(pm, nx, ny, nz);
72 
73  stk_classic::mesh::fem::FEMMetaData & fem_meta = fixture.m_fem_meta;
74  stk_classic::mesh::BulkData & bulk = fixture.m_bulk_data;
75 
76  // Put weights field on all entities
77  const stk_classic::mesh::EntityRank element_rank = fem_meta.element_rank();
78  const stk_classic::mesh::EntityRank face_rank = fem_meta.face_rank();
79  const stk_classic::mesh::EntityRank edge_rank = fem_meta.edge_rank();
80  const stk_classic::mesh::EntityRank node_rank = fem_meta.node_rank();
81  ScalarField & weight_field( fem_meta.declare_field< ScalarField >( "entity_weights" ) );
82  stk_classic::mesh::put_field(weight_field , element_rank , fem_meta.universal_part() );
83  stk_classic::mesh::put_field(weight_field , face_rank , fem_meta.universal_part() );
84  stk_classic::mesh::put_field(weight_field , edge_rank , fem_meta.universal_part() );
85  stk_classic::mesh::put_field(weight_field , node_rank , fem_meta.universal_part() );
86 
87  fem_meta.commit();
88 
89  // Configure our mesh on proc 0
90  std::vector<stk_classic::mesh::EntityId> my_element_ids;
91  if( 0 == p_rank )
92  {
93  for ( unsigned i = 0 ; i < nx*ny*nz; ++i )
94  my_element_ids.push_back(i+1);
95  }
96 
97  fixture.generate_mesh(my_element_ids);
98 
99  // create faces and edges for the mesh
101  stk_classic::mesh::create_adjacent_entities(bulk, add_parts);
102 
103  bulk.modification_begin();
104 
105  // Assign entity weights
106  if( 0 == p_rank )
107  {
108  // Get the faces on the x=0 plane and give them a characteristic weight
109  stk_classic::mesh::EntityVector selected_nodes;
110  stk_classic::mesh::EntityVector selected_faces;
111  stk_classic::mesh::EntityVector one_face;
112  for ( unsigned j = 0 ; j < ny; ++j )
113  for ( unsigned k = 0 ; k < nz; ++k )
114  {
115  selected_nodes.clear();
116  selected_nodes.push_back( fixture.node(0, j, k ) );
117  selected_nodes.push_back( fixture.node(0, j+1, k ) );
118  selected_nodes.push_back( fixture.node(0, j, k+1) );
119  selected_nodes.push_back( fixture.node(0, j+1, k+1) );
120  stk_classic::mesh::get_entities_through_relations(selected_nodes, face_rank, one_face);
121  selected_faces.push_back(one_face[0]);
122  }
123 
124  for( size_t iface = 0; iface < selected_faces.size(); ++iface )
125  {
126  stk_classic::mesh::Entity * face = selected_faces[iface];
127  double * const weight = stk_classic::mesh::field_data( weight_field, *face );
128  weight[0] = 10.0;
129  }
130 
131  // Get the edges on the boundary of the x=0 plane and give them a characteristic weight
132  stk_classic::mesh::EntityVector selected_edges;
133  stk_classic::mesh::EntityVector one_edge;
134  for ( unsigned j = 0 ; j < ny; ++j )
135  {
136  selected_nodes.clear();
137  selected_nodes.push_back( fixture.node(0, j, 0) );
138  selected_nodes.push_back( fixture.node(0, j+1, 0) );
139  stk_classic::mesh::get_entities_through_relations(selected_nodes, edge_rank, one_edge);
140  selected_edges.push_back(one_edge[0]);
141  selected_nodes.clear();
142  selected_nodes.push_back( fixture.node(0, j, nz) );
143  selected_nodes.push_back( fixture.node(0, j+1, nz) );
144  stk_classic::mesh::get_entities_through_relations(selected_nodes, edge_rank, one_edge);
145  selected_edges.push_back(one_edge[0]);
146  }
147  for ( unsigned k = 0 ; k < nz; ++k )
148  {
149  selected_nodes.clear();
150  selected_nodes.push_back( fixture.node(0, 0, k) );
151  selected_nodes.push_back( fixture.node(0, 0, k+1) );
152  stk_classic::mesh::get_entities_through_relations(selected_nodes, edge_rank, one_edge);
153  selected_edges.push_back(one_edge[0]);
154  selected_nodes.clear();
155  selected_nodes.push_back( fixture.node(0, ny, k) );
156  selected_nodes.push_back( fixture.node(0, ny, k+1) );
157  stk_classic::mesh::get_entities_through_relations(selected_nodes, edge_rank, one_edge);
158  selected_edges.push_back(one_edge[0]);
159  }
160  for( size_t iedge = 0; iedge < selected_edges.size(); ++iedge )
161  {
162  stk_classic::mesh::Entity * edge = selected_edges[iedge];
163  double * const weight = stk_classic::mesh::field_data( weight_field, *edge );
164  weight[0] = 100.0;
165  }
166 
167  // Finally, give the corner nodes of the x=0 plane a characteristic weight
168  selected_nodes.clear();
169  double * weight = stk_classic::mesh::field_data( weight_field, *fixture.node(0, 0, 0) );
170  weight[0] = 1000.0;
171  weight = stk_classic::mesh::field_data( weight_field, *fixture.node(0, ny, 0) );
172  weight[0] = 1000.0;
173  weight = stk_classic::mesh::field_data( weight_field, *fixture.node(0, 0, nz) );
174  weight[0] = 1000.0;
175  weight = stk_classic::mesh::field_data( weight_field, *fixture.node(0, ny, nz) );
176  weight[0] = 1000.0;
177 
178  // Assign element weights
179  for( size_t i = 0; i < my_element_ids.size(); ++i )
180  {
181  stk_classic::mesh::Entity * elem = bulk.get_entity(element_rank, my_element_ids[i]);
182  double * const e_weight = stk_classic::mesh::field_data( weight_field , *elem );
183  *e_weight = 1.0;
184  }
185  //
186  // Get the elements on the x=0 plane and sum in weights from relations
187  selected_nodes.clear();
188  for ( unsigned j = 0 ; j < ny+1; ++j )
189  for ( unsigned k = 0 ; k < nz+1; ++k )
190  selected_nodes.push_back( fixture.node(0, j, k) );
191 
192  std::vector<stk_classic::mesh::EntityRank> ranks;
193  ranks.push_back(face_rank);
194  ranks.push_back(edge_rank);
195  ranks.push_back(node_rank);
196  stk_classic::mesh::EntityVector selected_elems;
197  for ( unsigned j = 0 ; j < ny; ++j )
198  for ( unsigned k = 0 ; k < nz; ++k )
199  {
200  selected_nodes.clear();
201  selected_nodes.push_back( fixture.node(0, j, k ) );
202  selected_nodes.push_back( fixture.node(0, j+1, k ) );
203  selected_nodes.push_back( fixture.node(0, j, k+1) );
204  selected_nodes.push_back( fixture.node(0, j+1, k+1) );
205  stk_classic::mesh::get_entities_through_relations(selected_nodes, element_rank, one_face);
206  selected_elems.push_back(one_face[0]);
207  }
208  sum_element_weights_through_relations(selected_elems, weight_field, ranks);
209  }
210 
211  bulk.modification_end();
212 
213  // Use Zoltan to determine new partition
214  Teuchos::ParameterList emptyList;
215  stk_classic::rebalance::Zoltan zoltan_partition(pm, fixture.m_spatial_dimension, emptyList);
216 
217  stk_classic::rebalance::rebalance(bulk, fem_meta.universal_part(), &fixture.m_coord_field, &weight_field, zoltan_partition);
218 
219  const double imbalance_threshold = ( 3 == p_size )? 1.45 // Zoltan does not do so well for 3 procs
220  : 1.1; // ... but does pretty well for 2 and 4 procs
221  const bool do_rebal = imbalance_threshold < stk_classic::rebalance::check_balance(bulk, &weight_field, element_rank);
222 
223  if( 0 == p_rank )
224  std::cerr << std::endl
225  << "imbalance_threshold after rebalance = " << imbalance_threshold << ", " << do_rebal << std::endl;
226 
227 
228  // Check that we satisfy our threshhold
229  bool result = !do_rebal;
230 
231  // And verify that all dependent entities are on the same proc as their parent element
232  {
233  stk_classic::mesh::EntityVector entities;
234  stk_classic::mesh::Selector selector = fem_meta.locally_owned_part();
235 
236  get_selected_entities(selector, bulk.buckets(node_rank), entities);
237  result &= verify_dependent_ownership(element_rank, entities);
238 
239  get_selected_entities(selector, bulk.buckets(edge_rank), entities);
240  result &= verify_dependent_ownership(element_rank, entities);
241 
242  get_selected_entities(selector, bulk.buckets(face_rank), entities);
243  result &= verify_dependent_ownership(element_rank, entities);
244  }
245 
246  return result;
247 }
248 
249 } //namespace use_cases
250 } //namespace rebalance
251 } //namespace stk_classic
252 
253 
Part & locally_owned_part() const
Subset for the problem domain that is owned by the local process. Ghost entities are not members of t...
EntityRank face_rank() const
Returns the face rank which changes depending on spatial dimension.
FEMMetaData is a class that implements a Finite Element Method skin on top of the Sierra Tool Kit Met...
Definition: FEMMetaData.hpp:54
bool rebalance(mesh::BulkData &bulk_data, const mesh::Selector &selector, const VectorField *coord_ref, const ScalarField *elem_weight_ref, Partition &partition, const stk_classic::mesh::EntityRank rank=stk_classic::mesh::InvalidEntityRank)
Rebalance with a Partition object.
Definition: Rebalance.cpp:164
FieldTraits< field_type >::data_type * field_data(const field_type &f, const Bucket::iterator i)
Pointer to the field data array.
Definition: FieldData.hpp:116
EntityRank element_rank() const
Returns the element rank which is always equal to spatial dimension.
Part & universal_part() const
Universal subset for the problem domain. All other parts are a subset of the universal part...
This is a class for selecting buckets based on a set of meshparts and set logic.
Definition: Selector.hpp:112
const std::vector< Bucket * > & buckets(EntityRank rank) const
Query all buckets of a given entity rank.
Definition: BulkData.hpp:195
field_type & put_field(field_type &field, EntityRank entity_rank, const Part &part, const void *init_value=NULL)
Declare a field to exist for a given entity type and Part.
Entity * get_entity(EntityRank entity_rank, EntityId entity_id) const
Get entity with a given key.
Definition: BulkData.hpp:211
void get_selected_entities(const Selector &selector, const std::vector< Bucket * > &input_buckets, std::vector< Entity * > &entities)
Get entities in selected buckets (selected by the given selector instance), and sorted by ID...
Definition: GetEntities.cpp:77
unsigned parallel_machine_rank(ParallelMachine parallel_machine)
Member function parallel_machine_rank ...
Definition: Parallel.cpp:29
bool modification_end()
Parallel synchronization of modifications and transition to the guaranteed parallel consistent state...
bool modification_begin()
Begin a modification phase during which the mesh bulk data could become parallel inconsistent. This is a parallel synchronous call. The first time this method is called the mesh meta data is verified to be committed and parallel consistent. An exception is thrown if this verification fails.
Definition: BulkData.cpp:172
PairIterRelation relations() const
All Entity relations for which this entity is a member. The relations are ordered from lowest entity-...
Definition: Entity.hpp:161
unsigned parallel_machine_size(ParallelMachine parallel_machine)
Member function parallel_machine_size ...
Definition: Parallel.cpp:18
Manager for an integrated collection of entities, entity relations, and buckets of field data...
Definition: BulkData.hpp:49
void commit()
Commit the part and field declarations so that the meta data manager can be used to create mesh bulk ...
field_type & declare_field(const std::string &name, unsigned number_of_states=1)
Declare a field of the given field_type, test name, and number of states.
A fundamental unit within the discretization of a problem domain, including but not limited to nodes...
Definition: Entity.hpp:120
Sierra Toolkit.
EntityRank edge_rank() const
Returns the edge rank which changes depending on spatial dimension.
MPI_Comm ParallelMachine
Definition: Parallel.hpp:32
void get_entities_through_relations(const std::vector< Entity *> &entities, std::vector< Entity *> &entities_related)
Query which mesh entities have a relation to all of the input mesh entities.
Definition: Relation.cpp:156
For partitioning of mesh entities over a processing grid.
Static functions for dynamic load balancing.
EntityRank node_rank() const
Returns the node rank, which is always zero.
std::vector< Part *> PartVector
Collections of parts are frequently maintained as a vector of Part pointers.
Definition: Types.hpp:31