Flow123d  master-f44eb46
partitioning.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file partitioning.cc
15  * @brief
16  */
17 
18 #include "system/index_types.hh"
19 #include "mesh/partitioning.hh"
20 #include "la/sparse_graph.hh"
21 #include "la/distribution.hh"
22 #include "mesh/mesh.h"
23 #include "mesh/accessors.hh"
24 #include "mesh/range_wrapper.hh"
25 #include "mesh/neighbours.h"
26 
27 #include "petscao.h"
28 
29 
30 namespace IT = Input::Type;
32  return IT::Selection("GraphType",
33  "Different algorithms to make the sparse graph with weighted edges\n"
34  "from the multidimensional mesh. Main difference is dealing with \n"
35  "neighboring of elements of different dimension.")
36  .add_value(any_neighboring, "any_neighboring", "Add an edge for any pair of neighboring elements.")
37  .add_value(any_weight_lower_dim_cuts, "any_weight_lower_dim_cuts", "Same as before and assign higher weight to cuts of lower dimension in order to make them stick to one face.")
38  .add_value(same_dimension_neighboring, "same_dimension_neighboring", "Add an edge for any pair of neighboring elements of the same dimension (bad for matrix multiply).")
39  .close();
40 }
41 
43  return IT::Selection("PartTool", "Select the partitioning tool to use.")
44  .add_value(PETSc, "PETSc", "Use PETSc interface to various partitioning tools.")
45  .add_value(METIS, "METIS", "Use direct interface to Metis.")
46  .close();
47 }
48 
50  static IT::Record input_type = IT::Record("Partition","Setting for various types of mesh partitioning." )
51  .declare_key("tool", Partitioning::get_tool_sel(), IT::Default("\"METIS\""), "Software package used for partitioning. See corresponding selection.")
52  .declare_key("graph_type", Partitioning::get_graph_type_sel(), IT::Default("\"any_neighboring\""), "Algorithm for generating graph and its weights from a multidimensional mesh.")
53  .allow_auto_conversion("graph_type") // mainly in order to allow Default value for the whole record Partition
54  .close();
55  input_type.finish();
56 
57  return input_type;
58 }
59 
60 
62 : mesh_(mesh), in_(in), graph_(NULL), loc_part_(NULL), init_el_ds_(NULL)
63 {
65 }
66 
67 
68 
70  if (loc_part_) delete [] loc_part_;
71  loc_part_ = NULL;
72  if (init_el_ds_) delete init_el_ds_;
73  init_el_ds_ = NULL;
74  if (graph_) delete graph_;
75  graph_ = NULL;
76 }
77 
79  ASSERT_PTR(init_el_ds_).error("NULL initial distribution.");
80  return init_el_ds_;
81 }
82 
83 
84 
86  ASSERT_PTR(loc_part_).error("NULL local partitioning.");
87  return loc_part_;
88 }
89 
90 
91 
93 
94  Distribution edistr = graph_->get_distr();
95 
96  unsigned int e_idx;
97  unsigned int i_neigh;
98  int i_s, n_s;
99 
100  // Add nigbouring edges only for "any_*" graph types
101  bool neigh_on = ( in_.val<PartitionGraphType>("graph_type") != same_dimension_neighboring );
102 
103  for (auto ele : mesh_->elements_range()) {
104  // skip non-local elements
105  if ( !edistr.is_local( ele.idx() ) )
106  continue;
107 
108  // for all connected elements
109  for (unsigned int si=0; si<ele->n_sides(); si++) {
110  Edge edg = ele.side(si)->edge();
111 
112  for (unsigned int li=0; li<edg.n_sides(); li++) {
113  ASSERT(edg.side(li)->is_valid()).error("NULL side of edge.");
114  e_idx = edg.side(li)->element().idx();
115 
116  // for elements of connected elements, excluding element itself
117  if ( e_idx != ele.idx() ) {
118  graph_->set_edge( ele.idx(), e_idx );
119  }
120  }
121  }
122 
123  // include connections from lower dim. edge
124  // to the higher dimension
125  if (neigh_on) {
126  for (i_neigh = 0; i_neigh < ele->n_neighs_vb(); i_neigh++) {
127  n_s = ele->neigh_vb[i_neigh]->edge().n_sides();
128  for (i_s = 0; i_s < n_s; i_s++) {
129  e_idx = ele->neigh_vb[i_neigh]->edge().side(i_s)->element().idx();
130  graph_->set_edge( ele.idx(), e_idx );
131  graph_->set_edge( e_idx, ele.idx() );
132  }
133  }
134  }
135  }
136  graph_->finalize();
137 }
138 
139 
140 
142 
143 
144  // prepare dual graph
145  switch ( in_.val<PartitionTool>("tool")) {
146  case PETSc:
147  init_el_ds_ = new Distribution(DistributionBlock(), mesh_->n_elements(), mesh_->get_comm() ); // initial distr.
149 
150  break;
151  case METIS:
152  init_el_ds_ = new Distribution(DistributionLocalized(), mesh_->n_elements(), mesh_->get_comm() ); // initial distr.
154  break;
155  }
156  int mesh_size = mesh_->n_elements();
157  int num_of_procs = init_el_ds_->np();
158  if (mesh_size < num_of_procs) { // check if decomposing is possible
159  THROW( ExcDecomposeMesh() << EI_NElems( mesh_size ) << EI_NProcs( num_of_procs ) );
160  }
161 
163 
164  // compute partitioning
167  delete graph_; graph_ = NULL;
168 }
169 
170 
171 
172 
173 /**
174  * Old UGLY, PETSC dependent method for getting new numbering after partitioning.
175  *
176  * n_ids - given maximal ID used in id_4_old
177  * id_4_old - given array of size init_el_ds_.size() - assign ID to an old index
178  *
179  * new_ds - new distribution of elements according to current distributed partitioning loc_part_
180  * id_4_loc - IDs for local elements in new distribution, has size new_ds->lsize()
181  * new_4_id - for given ID, the new index, -1 for unknown IDs
182  *
183  */
184 void Partitioning::id_maps(int n_ids, LongIdx *id_4_old,
185  const Distribution &old_ds, LongIdx *loc_part,
186  Distribution * &new_ds, LongIdx * &id_4_loc, LongIdx * &new_4_id) {
187 
188  IS part, new_numbering;
189  unsigned int size = old_ds.size(); // whole size of distr. array
190  int new_counts[old_ds.np()];
191  AO new_old_ao;
192  int *old_4_new;
193  int i_loc;
194 
195  // make distribution and numbering
196  ISCreateGeneral(PETSC_COMM_WORLD, old_ds.lsize(), loc_part, PETSC_COPY_VALUES, &part); // global IS part.
197  ISPartitioningCount(part, old_ds.np(), new_counts); // new size of each proc
198 
199  new_ds = new Distribution((unsigned int *) new_counts, PETSC_COMM_WORLD); // new distribution
200  ISPartitioningToNumbering(part, &new_numbering); // new numbering
201  ISDestroy(&part);
202 
203  old_4_new = new int [size];
204  id_4_loc = new LongIdx [ new_ds->lsize() ];
205  new_4_id = new LongIdx [ n_ids + 1 ];
206 
207  // create whole new->old mapping on each proc
208  AOCreateBasicIS(new_numbering, PETSC_NULL, &new_old_ao); // app ordering= new; petsc ordering = old
209  ISDestroy(&new_numbering);
210  for (unsigned int i = 0; i < size; i++)
211  old_4_new[i] = i;
212  AOApplicationToPetsc(new_old_ao, size, old_4_new);
213  AODestroy(&(new_old_ao));
214 
215  // compute id_4_loc
216  i_loc = 0;
217 
218  for (unsigned int i_new = new_ds->begin(); i_new < new_ds->end(); i_new++) {
219  id_4_loc[i_loc++] = id_4_old[old_4_new[i_new]];
220  }
221  // compute row_4_id
222  for (i_loc = 0; i_loc <= n_ids; i_loc++)
223  new_4_id[i_loc] = -1; // ensure that all ids are initialized
224  for (unsigned int i_new = 0; i_new < size; i_new++)
225  new_4_id[id_4_old[old_4_new[i_new]]] = i_new;
226 
227 
228  delete [] old_4_new;
229 }
230 
231 
232 void Partitioning::id_maps(int n_ids, LongIdx *id_4_old, Distribution * &new_ds, LongIdx * &id_4_loc, LongIdx * &new_4_id) {
233  Partitioning::id_maps(n_ids, id_4_old, *init_el_ds_, loc_part_, new_ds, id_4_loc, new_4_id);
234 }
235 
236 
237 
238 shared_ptr< vector<int> > Partitioning::subdomain_id_field_data() {
239  ASSERT_PERMANENT(loc_part_).error("Partition is not yet computed.\n");
240  if (!seq_part_) {
241  unsigned int seq_size=mesh_->get_el_ds()->lsize();
242  seq_part_ = make_shared< vector<int> >(seq_size, mesh_->get_el_ds()->myp());
243  }
244 
245  return seq_part_;
246 
247 }
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
bool is_local(unsigned int idx) const
identify local index
unsigned int myp() const
get my processor
unsigned int begin(int proc) const
get starting local index
unsigned int lsize(int proc) const
get local size
unsigned int size() const
get global size
unsigned int np() const
get num of processors
unsigned int n_sides() const
Returns number of sides aligned with the edge.
Definition: accessors.hh:329
SideIter side(const unsigned int i) const
Gets side iterator of the i -th side.
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Definition: accessors.hh:215
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
const Ret val(const string &key) const
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
Record type proxy class.
Definition: type_record.hh:182
FinishStatus finish(FinishStatus finish_type=FinishStatus::regular_) override
Finish declaration of the Record type.
Definition: type_record.cc:243
virtual Record & allow_auto_conversion(const string &from_key)
Allows shorter input of the Record providing only value of the from_key given as the parameter.
Definition: type_record.cc:133
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
Template for classes storing finite set of named values.
Selection & add_value(const int value, const std::string &key, const std::string &description="", TypeBase::attribute_map attributes=TypeBase::attribute_map())
Adds one new value with name given by key to the Selection.
const Selection & close() const
Close the Selection, no more values can be added.
Distribution * get_el_ds() const
Definition: mesh.h:119
unsigned int n_elements() const
Definition: mesh.h:111
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
Definition: mesh.cc:1174
Definition: mesh.h:362
MPI_Comm get_comm() const
Definition: mesh.h:453
static const Input::Type::Selection & get_tool_sel()
Definition: partitioning.cc:42
const LongIdx * get_loc_part() const
Definition: partitioning.cc:85
~Partitioning()
Destructor.
Definition: partitioning.cc:69
void make_partition()
shared_ptr< vector< int > > subdomain_id_field_data()
Mesh * mesh_
The input mesh.
static const Input::Type::Record & get_input_type()
Definition: partitioning.cc:49
LongIdx * loc_part_
Partition numbers for local elements in original distribution of elements given be init_el_ds_.
static const Input::Type::Selection & get_graph_type_sel()
Input specification objects.
Definition: partitioning.cc:31
void make_element_connection_graph()
Definition: partitioning.cc:92
@ any_weight_lower_dim_cuts
Same as before and assign higher weight to cuts of lower dimension in order to make them stick to one...
@ same_dimension_neighboring
Add edge for any pair of neighboring elements of same dimension (bad for matrix multiply)
@ any_neighboring
Add edge for any pair of neighboring elements.
shared_ptr< vector< int > > seq_part_
Sequential partitioning for output.
const Distribution * get_init_distr() const
Definition: partitioning.cc:78
@ METIS
Use direct interface to Metis.
@ PETSc
Use PETSc interface to various partitioing tools.
Input::Record in_
Input Record accessor.
SparseGraph * graph_
Graph used to partitioning the mesh.
Partitioning(Mesh *mesh, Input::Record in)
Definition: partitioning.cc:61
Distribution * init_el_ds_
Original distribution of elements. Depends on type of partitioner.
void id_maps(int n_ids, LongIdx *id_4_old, Distribution *&new_ds, LongIdx *&id_4_loc, LongIdx *&new_4_id)
bool is_valid() const
Returns true if the side has assigned element.
Definition: accessors.hh:452
ElementAccessor< 3 > element() const
Returns iterator to the element of the side.
Edge edge() const
Returns pointer to the edge connected to the side.
Distribution get_distr()
void set_edge(const int a, const int b, int weight=1)
Definition: sparse_graph.cc:66
void finalize()
Make sparse graph structures: rows, adj.
Definition: sparse_graph.cc:98
virtual void partition(int *loc_part)=0
Support classes for parallel programing.
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: index_types.hh:24
Implementation of range helper class.
Distributed sparse graphs, partitioning.