Sierra Toolkit  Version of the Day
UseCase_Rebal_4.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_4.hpp>
10 
11 #include <stk_util/parallel/Parallel.hpp>
12 #include <stk_util/parallel/ParallelReduce.hpp>
13 
14 #include <stk_mesh/base/FieldData.hpp>
15 #include <stk_mesh/base/GetEntities.hpp>
16 #include <stk_mesh/base/MetaData.hpp>
17 #include <stk_mesh/base/BulkData.hpp>
18 
19 #include <stk_mesh/fem/CoordinateSystems.hpp>
20 #include <stk_mesh/fem/FEMMetaData.hpp>
21 #include <stk_mesh/fem/FEMHelpers.hpp>
22 #include <stk_mesh/fem/CreateAdjacentEntities.hpp>
23 
24 #include <stk_mesh/fixtures/HexFixture.hpp>
25 
29 
30 #include <stk_rebalance_utils/RebalanceUtils.hpp>
31 
32 //----------------------------------------------------------------------
33 
34 using namespace stk_classic::mesh::fixtures;
35 
37 
38 namespace stk_classic {
39 namespace rebalance {
40 namespace use_cases {
41 
42 class GreedySideset : public Partition {
43  public :
44  struct MeshInfo {
45  std::vector<mesh::Entity *> mesh_entities;
46  const VectorField * nodal_coord_ref ;
47  const ScalarField * elem_weight_ref;
48  std::vector<unsigned> dest_proc_ids ;
49 
51  MeshInfo():
52  nodal_coord_ref(NULL),
53  elem_weight_ref(NULL) {}
54 
56  ~MeshInfo() {}
57  };
58  explicit GreedySideset(ParallelMachine pm,
59  const stk_classic::mesh::PartVector & surfaces,
60  mesh::BulkData & bulk_data);
61  virtual ~GreedySideset();
62  virtual void reset_dest_proc_data();
63  virtual void set_mesh_info ( const std::vector<mesh::Entity *> &mesh_entities,
64  const VectorField * nodal_coord_ref,
65  const ScalarField * elem_weight_ref=NULL);
66  virtual void determine_new_partition(bool &RebalancingNeeded);
67  virtual unsigned num_elems() const;
68  virtual int get_new_partition(stk_classic::mesh::EntityProcVec &new_partition);
69  virtual bool partition_dependents_needed()const;
70  bool find_mesh_entity(const mesh::Entity * entity, unsigned & moid) const;
71  unsigned destination_proc(const unsigned moid) const;
72  void set_destination_proc(const unsigned moid, const unsigned proc );
73  MeshInfo mesh_information_;
74  unsigned total_number_entities_;
75  const stk_classic::mesh::PartVector & surfaces_;
76  mesh::BulkData & bulk_data_;
77 };
78 
79 GreedySideset::GreedySideset(ParallelMachine pm,
80  const stk_classic::mesh::PartVector & surfaces,
81  mesh::BulkData & bulk_data) :
82  stk_classic::rebalance::Partition(pm),
83  mesh_information_(),
84  surfaces_(surfaces),
85  bulk_data_(bulk_data) {}
86 GreedySideset::~GreedySideset() {}
87 void GreedySideset::reset_dest_proc_data() {
88  const int proc = parallel_machine_rank(comm_);
89  const unsigned size = mesh_information_.mesh_entities.size();
90  mesh_information_.dest_proc_ids.assign(size, proc);
91 }
92 void GreedySideset::set_mesh_info ( const std::vector<mesh::Entity *> &mesh_entities,
93  const VectorField * nodal_coord_ref,
94  const ScalarField * elem_weight_ref){
95  MeshInfo mesh_info;
96 
97  /* Keep track of the total number of elements. */
98  total_number_entities_ = mesh_entities.size();
99 
100  mesh_info.mesh_entities = mesh_entities;
101  mesh_info.nodal_coord_ref = nodal_coord_ref;
102  mesh_info.elem_weight_ref = elem_weight_ref;
103 
109  mesh_info.dest_proc_ids.assign(mesh_entities.size(), stk_classic::parallel_machine_rank(comm_));
110 
111  mesh_information_ = mesh_info;
112 }
113 
114 unsigned GreedySideset::num_elems() const {return total_number_entities_ ;}
115 int GreedySideset::get_new_partition(stk_classic::mesh::EntityProcVec &new_partition){
116 std::vector<mesh::Entity*>::iterator i=mesh_information_.mesh_entities.begin();
117 std::vector<unsigned> ::iterator j=mesh_information_.dest_proc_ids.begin();
118  for (;i != mesh_information_.mesh_entities.end(),
119  j != mesh_information_.dest_proc_ids.end();
120  ++i,++j) {
121  mesh::Entity * mesh_entity = *i;
122  unsigned proc = *j;
123  mesh::EntityProc et(mesh_entity, proc);
124  new_partition.push_back(et);
125  }
126  return 0;
127 }
128 bool GreedySideset::partition_dependents_needed()const{return true;}
129 
130 bool GreedySideset::find_mesh_entity(const mesh::Entity * entity, unsigned & moid) const
131 {
132  unsigned len = mesh_information_.mesh_entities.size();
133  for(moid = 0; moid < len; ++moid)
134  {
135  if(mesh_information_.mesh_entities[moid] == entity) return true;
136  }
137  return false;
138 }
139 unsigned GreedySideset::destination_proc(const unsigned moid) const
140 {
141  return mesh_information_.dest_proc_ids[ moid ];
142 }
143 void GreedySideset::set_destination_proc(const unsigned moid,
144  const unsigned proc )
145 {
146  mesh_information_.dest_proc_ids[ moid ] = proc;
147 }
148 
149 
150 void GreedySideset::determine_new_partition(bool &RebalancingNeeded) {
151 
152  reset_dest_proc_data();
153 
155  const stk_classic::mesh::EntityRank side_rank = fem_meta.side_rank();
156  const stk_classic::mesh::EntityRank elem_rank = fem_meta.element_rank();
157 
158  // Select active ghosted side faces.
159  stk_classic::mesh::Selector selector(!fem_meta.locally_owned_part() &
160  stk_classic::mesh::selectIntersection(surfaces_));
161 
162  mesh::EntityVector sides;
163  mesh::get_selected_entities(selector, bulk_data_.buckets(side_rank), sides);
164 
165  const unsigned p_rank = bulk_data_.parallel_rank();
166  size_t local_changes = 0;
167  const unsigned nSide = sides.size();
168  for(unsigned iSide = 0; iSide < nSide; ++iSide)
169  {
170  const mesh::Entity & side = *sides[iSide];
171  const unsigned sideProc = side.owner_rank();
172  ThrowRequireMsg(sideProc!=p_rank,
173  "When iterating Non-locally owned sides, found a locally owned side.");
174 
175  stk_classic::mesh::PairIterRelation iElem = side.relations(elem_rank);
176  for ( ; iElem.first != iElem.second; ++iElem.first ) {
177  const mesh::Entity & elem = *iElem.first->entity();
178  unsigned moid;
179  const bool mesh_entity_found = find_mesh_entity(&elem, moid);
180  if (mesh_entity_found) {
181  const unsigned elemProc = elem.owner_rank();
182  ThrowRequireMsg(elemProc==p_rank,
183  "When iterating locally owned elements, found a non-locally owned element.");
184  const unsigned destProc = destination_proc(moid);
185  ThrowRequireMsg(destProc==p_rank || destProc==sideProc,
186  " Sanity check failed: "
187  "It's possible that an element is connected to "
188  "two sides that are owned by different procs. We don't "
189  "yet handle that situation here but we can, at least, "
190  "detect it. ");
191  if(elemProc != sideProc)
192  {
193  ++local_changes;
194  set_destination_proc(moid, sideProc);
195  }
196  }
197  }
198  }
199  size_t global_changes = 0;
200  stk_classic::all_reduce_sum (comm_, &local_changes, &global_changes, 1);
201  RebalancingNeeded = global_changes > 0;
202 }
203 
204 enum { nx = 3, ny = 3 };
205 
206 bool test_greedy_sideset ( stk_classic::ParallelMachine comm )
207 {
208  unsigned spatial_dimension = 2;
209  std::vector<std::string> rank_names = stk_classic::mesh::fem::entity_rank_names(spatial_dimension);
211  fem_meta.FEM_initialize(spatial_dimension, rank_names);
213  stk_classic::mesh::BulkData bulk_data( meta_data , comm , 100 );
214  const stk_classic::mesh::EntityRank element_rank = fem_meta.element_rank();
215  const stk_classic::mesh::EntityRank node_rank = fem_meta.node_rank();
216 
217  stk_classic::mesh::fem::CellTopology quad_top(shards::getCellTopologyData<shards::Quadrilateral<4> >());
218  stk_classic::mesh::fem::CellTopology line_top(shards::getCellTopologyData<shards::Line<2> >());
219  stk_classic::mesh::Part & quad_part( fem_meta.declare_part("quad", quad_top ) );
220  stk_classic::mesh::Part & side_part( fem_meta.declare_part("line", line_top ) );
221  VectorField & coord_field( fem_meta.declare_field< VectorField >( "coordinates" ) );
222  ScalarField & weight_field( fem_meta.declare_field< ScalarField >( "element_weights" ) );
223 
224  stk_classic::mesh::put_field( coord_field , node_rank , fem_meta.universal_part() );
225  stk_classic::mesh::put_field(weight_field , element_rank , fem_meta.universal_part() );
226 
227  fem_meta.commit();
228  const unsigned p_rank = bulk_data.parallel_rank();
229  bulk_data.modification_begin();
230 
231  if ( !p_rank ) {
232 
233  std::vector<std::vector<stk_classic::mesh::Entity*> > quads(nx);
234  for ( unsigned ix = 0 ; ix < nx ; ++ix ) quads[ix].resize(ny);
235 
236  const unsigned nnx = nx + 1 ;
237  for ( unsigned iy = 0 ; iy < ny ; ++iy ) {
238  for ( unsigned ix = 0 ; ix < nx ; ++ix ) {
239  stk_classic::mesh::EntityId elem = 1 + ix + iy * nx ;
240  stk_classic::mesh::EntityId nodes[4] ;
241  nodes[0] = 1 + ix + iy * nnx ;
242  nodes[1] = 2 + ix + iy * nnx ;
243  nodes[2] = 2 + ix + ( iy + 1 ) * nnx ;
244  nodes[3] = 1 + ix + ( iy + 1 ) * nnx ;
245 
246  stk_classic::mesh::Entity &q = stk_classic::mesh::fem::declare_element( bulk_data , quad_part , elem , nodes );
247  quads[ix][iy] = &q;
248  }
249  }
250 
251  for ( unsigned iy = 0 ; iy < ny ; ++iy ) {
252  for ( unsigned ix = 0 ; ix < nx ; ++ix ) {
253  stk_classic::mesh::EntityId elem = 1 + ix + iy * nx ;
254  stk_classic::mesh::Entity * e = bulk_data.get_entity( element_rank, elem );
255  double * const e_weight = stk_classic::mesh::field_data( weight_field , *e );
256  *e_weight = 1.0;
257  }
258  }
259  for ( unsigned iy = 0 ; iy <= ny ; ++iy ) {
260  for ( unsigned ix = 0 ; ix <= nx ; ++ix ) {
261  stk_classic::mesh::EntityId nid = 1 + ix + iy * nnx ;
262  stk_classic::mesh::Entity * n = bulk_data.get_entity( node_rank, nid );
263  double * const coord = stk_classic::mesh::field_data( coord_field , *n );
264  coord[0] = .1*ix;
265  coord[1] = .1*iy;
266  coord[2] = 0;
267  }
268  }
269  }
270 
271  bulk_data.modification_end();
272 
273  // create some sides and faces to rebalance.
275  stk_classic::mesh::create_adjacent_entities(bulk_data, add_parts);
276 
277  bulk_data.modification_begin();
278 
279  const stk_classic::mesh::PartVector surfaces(1, &side_part);
280  {
281  const stk_classic::mesh::PartVector empty_remove_parts;
283  const stk_classic::mesh::EntityRank side_rank = fmeta.side_rank();
285  mesh::EntityVector sides;
286  mesh::get_selected_entities(selector2, bulk_data.buckets(side_rank), sides);
287 
288  const unsigned nSide = sides.size();
289  for(unsigned iSide = 0; iSide < nSide; ++iSide)
290  {
291  mesh::Entity & side = *sides[iSide];
292  if (side.identifier()==7) {
293  bulk_data.change_entity_parts(side, surfaces, empty_remove_parts);
294  }
295  }
296  }
297  bulk_data.modification_end();
298 
299  // Zoltan partition is specialized form a virtual base class, stk_classic::rebalance::Partition.
300  // Other specializations are possible.
301  Teuchos::ParameterList emptyList;
302  stk_classic::rebalance::Zoltan zoltan_partition(comm, spatial_dimension, emptyList);
303  stk_classic::mesh::Selector selector3(fem_meta.locally_owned_part());
304  stk_classic::rebalance::rebalance(bulk_data, selector3, &coord_field, NULL, zoltan_partition);
305  {
306  const int print_stats = 1;
307  int nentity = 0;
308  double entity_wgt = 0;
309  int ncuts = 0;
310  double cut_wgt = 0;
311  int nboundary = 0;
312  int nadj = 0;
313  const int ierr = zoltan_partition.evaluate (print_stats, &nentity, &entity_wgt, &ncuts, &cut_wgt, &nboundary, &nadj);
314  std::cout <<" Information returned from the Zoltan evaluate function:"<<std::endl;
315  std::cout <<" Error Code: :"<<ierr <<std::endl;
316  std::cout <<" Number of entities: :"<<nentity <<std::endl;
317  std::cout <<" Number of cuts: :"<<ncuts <<std::endl;
318  std::cout <<" Cut Weight: :"<<cut_wgt <<std::endl;
319  std::cout <<" Number on Boundary: :"<<nboundary <<std::endl;
320  std::cout <<" Number Adjancent: :"<<nadj <<std::endl;
321  {
323  const stk_classic::mesh::EntityRank side_rank = fmeta.side_rank();
324  const stk_classic::mesh::EntityRank elem_rank = fmeta.element_rank();
325  const mesh::Entity *s = bulk_data.get_entity(side_rank,7);
326  if (s) {
327  const mesh::Entity & side = *s;
328  if (p_rank == side.owner_rank()) {
329  stk_classic::mesh::PairIterRelation iElem = side.relations(elem_rank);
330  for ( ; iElem.first != iElem.second; ++iElem.first ) {
331  const mesh::Entity & elem = *iElem.first->entity();
332  const unsigned elemProc = elem.owner_rank();
333  if (elemProc!=p_rank) {
334  std::cout <<p_rank<<" Good: Found element of of side 7 not owned."
335  <<" Element "<<elemProc
336  <<" on processor "<<p_rank
337  <<" owned by "<<elemProc<<std::endl;
338  }
339  }
340  }
341  }
342  }
343  }
344 
345  const double imbalance_threshold = 1.5;
346  bool do_rebal = imbalance_threshold < stk_classic::rebalance::check_balance(bulk_data, NULL, element_rank);
347 
348  if( !p_rank )
349  std::cout << std::endl
350  << "Use Case 4: imbalance_threshold after rebalance 1 = " << imbalance_threshold <<", "<<do_rebal << std::endl;
351 
352  {
353  stk_classic::rebalance::use_cases::GreedySideset greedy_sideset(comm, surfaces, bulk_data);
354  stk_classic::mesh::Selector selector4(fem_meta.locally_owned_part());
355  stk_classic::rebalance::rebalance(bulk_data, selector4, &coord_field, NULL, greedy_sideset);
356  }
357 
358  do_rebal = imbalance_threshold < stk_classic::rebalance::check_balance(bulk_data, NULL, element_rank);
359 
360  if( !p_rank )
361  std::cout << std::endl
362  << "Use Case 4: imbalance_threshold after rebalance 2 = " << imbalance_threshold <<", "<<do_rebal << std::endl;
363  {
365  const stk_classic::mesh::EntityRank side_rank = fmeta.side_rank();
366  const stk_classic::mesh::EntityRank elem_rank = fmeta.element_rank();
367  mesh::Entity *s = bulk_data.get_entity(side_rank,7);
368  if (s) {
369  mesh::Entity & side = *s;
370  if (p_rank == side.owner_rank()) {
371  stk_classic::mesh::PairIterRelation iElem = side.relations(elem_rank);
372  for ( ; iElem.first != iElem.second; ++iElem.first ) {
373  const mesh::Entity & elem = *iElem.first->entity();
374  const unsigned elemProc = elem.owner_rank();
375  if (elemProc!=p_rank) {
376  std::cerr <<p_rank<<" Error: Found element of of side 7 not owned:"<<elemProc<<std::endl;
377  }
378  ThrowRequireMsg(elemProc==p_rank,
379  "Use case 4 error check failed. Found element of of side 7 not owned.");
380  }
381  }
382  }
383  }
384  // Check that we satisfy our threshhold
385  bool result = !do_rebal ;
386 
387  // And verify that all dependent entities are on the same proc as their parent element
388  {
389  stk_classic::mesh::EntityVector entities;
390  stk_classic::mesh::Selector selector5 = fem_meta.locally_owned_part();
391 
392  get_selected_entities(selector5, bulk_data.buckets(node_rank), entities);
393  result &= verify_dependent_ownership(element_rank, entities);
394  get_selected_entities(selector5, bulk_data.buckets(fem_meta.side_rank()), entities);
395  result &= verify_dependent_ownership(element_rank, entities);
396  }
397 
398 
399 
400  return result;
401 }
402 
403 } //namespace use_cases
404 } //namespace rebalance
405 } //namespace stk_classic
406 
407 
Part & locally_owned_part() const
Subset for the problem domain that is owned by the local process. Ghost entities are not members of t...
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
The manager of an integrated collection of parts and fields.
Definition: MetaData.hpp:56
EntityRank side_rank() const
Returns the side rank which changes depending on spatial dimension.
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
Entity & declare_element(BulkData &mesh, Part &part, const EntityId elem_id, const EntityId node_id[])
Declare an element member of a Part with a CellTopology and nodes conformal to that topology...
Definition: FEMHelpers.cpp:72
void all_reduce_sum(ParallelMachine comm, const double *local, double *global, unsigned count)
Parallel summation to all processors.
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
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.
std::pair< Entity *, unsigned > EntityProc
Pairing of an entity with a processor rank.
Definition: Types.hpp:111
An application-defined subset of a problem domain.
Definition: Part.hpp:49
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
Part & declare_part(const std::string &name, fem::CellTopology cell_topology)
Declare a part with a given cell topology.
Manager for an integrated collection of entities, entity relations, and buckets of field data...
Definition: BulkData.hpp:49
static MetaData & get_meta_data(FEMMetaData &fem_meta)
Getter for MetaData off of a FEMMetaData object.
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.
MPI_Comm ParallelMachine
Definition: Parallel.hpp:32
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
static FEMMetaData & get(const MetaData &meta)
Getter for FEMMetaData off of a MetaData object.
void FEM_initialize(size_t spatial_dimension, const std::vector< std::string > &in_entity_rank_names=std::vector< std::string >())
Initialize the spatial dimension and an optional list of entity rank names associated with each rank...
Definition: FEMMetaData.cpp:61