Panzer  Version of the Day
Panzer_LocalPartitioningUtilities.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Panzer: A partial differential equation assembly
5 // engine for strongly coupled complex multiphysics systems
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39 // Eric C. Cyr (eccyr@sandia.gov)
40 // ***********************************************************************
41 // @HEADER
42 
44 
45 #include "Teuchos_RCP.hpp"
46 #include "Teuchos_Comm.hpp"
47 #include "Teuchos_Assert.hpp"
48 
49 #include "Panzer_Workset_Builder.hpp"
51 
52 #include "Phalanx_KokkosDeviceTypes.hpp"
53 
54 #include <unordered_set>
55 #include <unordered_map>
56 
57 namespace panzer
58 {
59 
60 namespace partitioning_utilities
61 {
62 
63 template<typename LO, typename GO>
64 void
66  const std::vector<LO> & owned_parent_cells,
68 {
69 
70  // The goal of this function is to fill a LocalMeshInfoBase (sub_info) with
71  // a subset of cells from a given parent LocalMeshInfoBase (parent_info)
72 
73  // Note: owned_parent_cells are the owned cells for sub_info in the parent_info's indexing scheme
74  // We need to generate sub_info's ghosts and figure out the virtual cells
75 
76  // Note: We will only handle a single ghost layer
77 
78  // Note: We assume owned_parent_cells are owned cells of the parent
79  // i.e. owned_parent_indexes cannot refer to ghost or virtual cells in parent_info
80 
81  // Note: This function works with inter-face connectivity. NOT node connectivity.
82 
83  const int num_owned_cells = owned_parent_cells.size();
84  TEUCHOS_TEST_FOR_EXCEPT_MSG(num_owned_cells == 0, "panzer::partitioning_utilities::setupSubLocalMeshInfo : Input parent subcells must exist (owned_parent_cells)");
85 
86  const int num_parent_owned_cells = parent_info.num_owned_cells;
87  TEUCHOS_TEST_FOR_EXCEPT_MSG(num_parent_owned_cells <= 0, "panzer::partitioning_utilities::setupSubLocalMeshInfo : Input parent info must contain owned cells");
88 
89  const int num_parent_ghstd_cells = parent_info.num_ghstd_cells;
90 
91  const int num_parent_total_cells = parent_info.num_owned_cells + parent_info.num_ghstd_cells + parent_info.num_virtual_cells;
92 
93  // Just as a precaution, make sure the parent_info is setup properly
94  TEUCHOS_ASSERT(static_cast<int>(parent_info.cell_to_faces.extent(0)) == num_parent_total_cells);
95  const int num_faces_per_cell = parent_info.cell_to_faces.extent(1);
96 
97  // The first thing to do is construct a vector containing the parent cell indexes of all
98  // owned, ghstd, and virtual cells
99  std::vector<LO> ghstd_parent_cells;
100  std::vector<LO> virtual_parent_cells;
101  {
102 
103  // We grab all of the owned cells and put their global indexes into sub_info
104  // We also put all of the owned cell indexes in the parent's indexing scheme into a set to use for lookups
105  std::unordered_set<LO> owned_parent_cells_set(owned_parent_cells.begin(), owned_parent_cells.end());
106 
107  // We need to create a list of ghstd and virtual cells
108  // We do this by running through sub_cell_indexes
109  // and looking at the neighbors to find neighbors that are not owned
110 
111  // Virtual cells are defined as cells with indexes outside of the range of owned_cells and ghstd_cells
112  const int virtual_parent_cell_offset = num_parent_owned_cells + num_parent_ghstd_cells;
113 
114  std::unordered_set<LO> ghstd_parent_cells_set;
115  std::unordered_set<LO> virtual_parent_cells_set;
116  for(int i=0;i<num_owned_cells;++i){
117  const LO parent_cell_index = owned_parent_cells[i];
118  for(int local_face_index=0;local_face_index<num_faces_per_cell;++local_face_index){
119  const LO parent_face = parent_info.cell_to_faces(parent_cell_index, local_face_index);
120 
121  // Sidesets can have owned cells that border the edge of the domain (i.e. parent_face == -1)
122  // If we are at the edge of the domain, we can ignore this face.
123  if(parent_face < 0)
124  continue;
125 
126  // Find the side index for neighbor cell with respect to the face
127  const LO neighbor_parent_side = (parent_info.face_to_cells(parent_face,0) == parent_cell_index) ? 1 : 0;
128 
129  // Get the neighbor cell index in the parent's indexing scheme
130  const LO neighbor_parent_cell = parent_info.face_to_cells(parent_face, neighbor_parent_side);
131 
132  // If the face exists, then the neighbor should exist
133  TEUCHOS_ASSERT(neighbor_parent_cell >= 0);
134 
135  // We can easily check if this is a virtual cell
136  if(neighbor_parent_cell >= virtual_parent_cell_offset){
137  virtual_parent_cells_set.insert(neighbor_parent_cell);
138  } else if(neighbor_parent_cell >= num_parent_owned_cells){
139  // This is a quick check for a ghost cell
140  // This branch only exists to cut down on the number of times the next branch (much slower) is called
141  ghstd_parent_cells_set.insert(neighbor_parent_cell);
142  } else {
143  // There is still potential for this to be a ghost cell with respect to 'our' cells
144  // The only way to check this is with a super slow lookup call
145  if(owned_parent_cells_set.find(neighbor_parent_cell) == owned_parent_cells_set.end()){
146  // The neighbor cell is not owned by 'us', therefore it is a ghost
147  ghstd_parent_cells_set.insert(neighbor_parent_cell);
148  }
149  }
150  }
151  }
152 
153  // We now have a list of the owned, ghstd, and virtual cells in the parent's indexing scheme.
154  // We will take the 'unordered_set's ordering for the the sub-indexing scheme
155 
156  ghstd_parent_cells.assign(ghstd_parent_cells_set.begin(), ghstd_parent_cells_set.end());
157  virtual_parent_cells.assign(virtual_parent_cells_set.begin(), virtual_parent_cells_set.end());
158 
159  }
160 
161  const int num_ghstd_cells = ghstd_parent_cells.size();
162  const int num_virtual_cells = virtual_parent_cells.size();
163  const int num_real_cells = num_owned_cells + num_ghstd_cells;
164  const int num_total_cells = num_real_cells + num_virtual_cells;
165 
166  std::vector<LO> all_parent_cells;
167  all_parent_cells.insert(all_parent_cells.end(), owned_parent_cells.begin(), owned_parent_cells.end());
168  all_parent_cells.insert(all_parent_cells.end(), ghstd_parent_cells.begin(), ghstd_parent_cells.end());
169  all_parent_cells.insert(all_parent_cells.end(), virtual_parent_cells.begin(), virtual_parent_cells.end());
170 
171  sub_info.num_owned_cells = owned_parent_cells.size();
172  sub_info.num_ghstd_cells = ghstd_parent_cells.size();
173  sub_info.num_virtual_cells = virtual_parent_cells.size();
174 
175  // We now have the indexing order for our sub_info
176 
177  // Just as a precaution, make sure the parent_info is setup properly
178  TEUCHOS_ASSERT(static_cast<int>(parent_info.cell_vertices.extent(0)) == num_parent_total_cells);
179  TEUCHOS_ASSERT(static_cast<int>(parent_info.local_cells.extent(0)) == num_parent_total_cells);
180  TEUCHOS_ASSERT(static_cast<int>(parent_info.global_cells.extent(0)) == num_parent_total_cells);
181 
182  const int num_vertices_per_cell = parent_info.cell_vertices.extent(1);
183  const int num_dims = parent_info.cell_vertices.extent(2);
184 
185  // Fill owned, ghstd, and virtual cells: global indexes, local indexes and vertices
186  sub_info.global_cells = Kokkos::View<GO*>("global_cells", num_total_cells);
187  sub_info.local_cells = Kokkos::View<LO*>("local_cells", num_total_cells);
188  sub_info.cell_vertices = Kokkos::View<double***, PHX::Device>("cell_vertices", num_total_cells, num_vertices_per_cell, num_dims);
189  for(int cell=0;cell<num_total_cells;++cell){
190  const LO parent_cell = all_parent_cells[cell];
191  sub_info.global_cells(cell) = parent_info.global_cells(parent_cell);
192  sub_info.local_cells(cell) = parent_info.local_cells(parent_cell);
193  for(int vertex=0;vertex<num_vertices_per_cell;++vertex){
194  for(int dim=0;dim<num_dims;++dim){
195  sub_info.cell_vertices(cell,vertex,dim) = parent_info.cell_vertices(parent_cell,vertex,dim);
196  }
197  }
198  }
199 
200  // Now for the difficult part
201 
202  // We need to create a new face indexing scheme from the old face indexing scheme
203 
204  // Create an auxiliary list with all cells - note this preserves indexing
205 
206  struct face_t{
207  face_t(LO c0, LO c1, LO sc0, LO sc1)
208  {
209  cell_0=c0;
210  cell_1=c1;
211  subcell_index_0=sc0;
212  subcell_index_1=sc1;
213  }
214  LO cell_0;
215  LO cell_1;
216  LO subcell_index_0;
217  LO subcell_index_1;
218  };
219 
220 
221  // First create the faces
222  std::vector<face_t> faces;
223  {
224 
225  // faces_set: cell_0, subcell_index_0, cell_1, subcell_index_1
226  std::unordered_map<LO,std::unordered_map<LO, std::pair<LO,LO> > > faces_set;
227  for(int owned_cell=0;owned_cell<num_owned_cells;++owned_cell){
228  const LO owned_parent_cell = owned_parent_cells[owned_cell];
229  for(int local_face=0;local_face<num_faces_per_cell;++local_face){
230  const LO parent_face = parent_info.cell_to_faces(owned_parent_cell,local_face);
231 
232  // Skip faces at the edge of the domain
233  if(parent_face<0)
234  continue;
235 
236  // Get the cell on the other side of the face
237  const LO neighbor_side = (parent_info.face_to_cells(parent_face,0) == owned_parent_cell) ? 1 : 0;
238 
239  const LO neighbor_parent_cell = parent_info.face_to_cells(parent_face, neighbor_side);
240  const LO neighbor_subcell_index = parent_info.face_to_lidx(parent_face, neighbor_side);
241 
242  // Convert parent cell index into sub cell index
243  auto itr = std::find(all_parent_cells.begin(), all_parent_cells.end(), neighbor_parent_cell);
244  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr == all_parent_cells.end(), "panzer_stk::setupSubLocalMeshInfo : Neighbor cell was not found in owned, ghosted, or virtual cells");
245 
246  const LO neighbor_cell = std::distance(all_parent_cells.begin(), itr);
247 
248  LO cell_0, cell_1, subcell_index_0, subcell_index_1;
249  if(owned_cell < neighbor_cell){
250  cell_0 = owned_cell;
251  subcell_index_0 = local_face;
252  cell_1 = neighbor_cell;
253  subcell_index_1 = neighbor_subcell_index;
254  } else {
255  cell_1 = owned_cell;
256  subcell_index_1 = local_face;
257  cell_0 = neighbor_cell;
258  subcell_index_0 = neighbor_subcell_index;
259  }
260 
261  // Add this interface to the set of faces - smaller cell index is 'left' (or '0') side of face
262  faces_set[cell_0][subcell_index_0] = std::pair<LO,LO>(cell_1, subcell_index_1);
263  }
264  }
265 
266  for(const auto & cell_pair : faces_set){
267  const LO cell_0 = cell_pair.first;
268  for(const auto & subcell_pair : cell_pair.second){
269  const LO subcell_index_0 = subcell_pair.first;
270  const LO cell_1 = subcell_pair.second.first;
271  const LO subcell_index_1 = subcell_pair.second.second;
272  faces.push_back(face_t(cell_0,cell_1,subcell_index_0,subcell_index_1));
273  }
274  }
275  }
276 
277  const int num_faces = faces.size();
278 
279  sub_info.face_to_cells = Kokkos::View<LO*[2]>("face_to_cells", num_faces);
280  sub_info.face_to_lidx = Kokkos::View<LO*[2]>("face_to_lidx", num_faces);
281  sub_info.cell_to_faces = Kokkos::View<LO**>("cell_to_faces", num_total_cells, num_faces_per_cell);
282 
283  // Default the system with invalid cell index
284  Kokkos::deep_copy(sub_info.cell_to_faces, -1);
285 
286  for(int face_index=0;face_index<num_faces;++face_index){
287  const face_t & face = faces[face_index];
288 
289  sub_info.face_to_cells(face_index,0) = face.cell_0;
290  sub_info.face_to_cells(face_index,1) = face.cell_1;
291 
292  sub_info.cell_to_faces(face.cell_0,face.subcell_index_0) = face_index;
293  sub_info.cell_to_faces(face.cell_1,face.subcell_index_1) = face_index;
294 
295  sub_info.face_to_lidx(face_index,0) = face.subcell_index_0;
296  sub_info.face_to_lidx(face_index,1) = face.subcell_index_1;
297 
298  }
299 
300  // Complete.
301 
302 }
303 
304 
305 
306 
307 
308 template<typename LO, typename GO>
309 void
311  const int splitting_size,
312  std::vector<panzer::LocalMeshPartition<LO,GO> > & partitions)
313 {
314 
315  // Make sure the splitting size makes sense
316  TEUCHOS_ASSERT(splitting_size != 0);
317 
318  // This is not a partitioning scheme.
319  // This just breaks the mesh_info into equally sized chunks and ignores connectivity
320  // This means that the cells in the partition probably won't be nearby each other - leads to an excess of ghost cells
321 
322  const LO num_owned_cells = mesh_info.num_owned_cells;
323 
324  if(splitting_size < 0){
325 
326  // Just one chunk
327  std::vector<LO> partition_cells;
328  partition_cells.reserve(mesh_info.num_owned_cells);
329  for(LO i=0;i<mesh_info.num_owned_cells;++i){
330  partition_cells.push_back(i);
331  }
332 
333  // It seems this can happen
334  if(partition_cells.size() > 0){
335  partitions.push_back(panzer::LocalMeshPartition<LO,GO>());
336  setupSubLocalMeshInfo(mesh_info,partition_cells,partitions.back());
337  }
338 
339  } else {
340 
341  std::vector<LO> partition_cells;
342  partition_cells.reserve(splitting_size);
343 
344  // There should be at least one partition
345  const LO num_partitions = mesh_info.num_owned_cells / splitting_size + ((mesh_info.num_owned_cells % splitting_size == 0) ? 0 : 1);
346 
347  LO partition_start_index = 0;
348  for(LO partition=0;partition<num_partitions;++partition){
349 
350  // Make sure end of partition maxes out at num_owned_cells
351  const LO partition_end_index = std::min(partition_start_index + splitting_size, num_owned_cells);
352 
353  partition_cells.resize(partition_end_index - partition_start_index,-1);
354  for(int i=partition_start_index;i<partition_end_index;++i){
355  partition_cells[i] = i;
356  }
357 
358  // It seems this can happen
359  if(partition_cells.size() > 0){
360  partitions.push_back(panzer::LocalMeshPartition<LO,GO>());
361  setupSubLocalMeshInfo(mesh_info,partition_cells,partitions.back());
362  }
363 
364  partition_start_index = partition_end_index;
365  }
366 
367  }
368 
369 }
370 
371 template<typename LO, typename GO>
372 void
374  const size_t /* requested_partition_size */,
375  std::vector<panzer::LocalMeshPartition<LO,GO>>& /* partitions */)
376 {
377  // Not yet sure how to do this
378  TEUCHOS_ASSERT(false);
379 }
380 
381 }
382 
383 template <typename LO, typename GO>
384 void
386  const panzer::WorksetDescriptor & description,
387  std::vector<panzer::LocalMeshPartition<LO,GO> > & partitions)
388 {
389 
390  // We have to make sure that the partitioning is possible
391  TEUCHOS_ASSERT(description.getWorksetSize() != panzer::WorksetSizeType::CLASSIC_MODE);
392  TEUCHOS_ASSERT(description.getWorksetSize() != 0);
393 
394  // This could just return, but it would be difficult to debug why no partitions were returned
395  TEUCHOS_ASSERT(description.requiresPartitioning());
396 
397  const std::string & element_block_name = description.getElementBlock();
398 
399  // We have two processes for in case this is a sideset or element block
400  if(description.useSideset()){
401 
402  // If the element block doesn't exist, there are no partitions to create
403  if(mesh_info.sidesets.find(element_block_name) == mesh_info.sidesets.end())
404  return;
405  const auto & sideset_map = mesh_info.sidesets.at(element_block_name);
406 
407  const std::string & sideset_name = description.getSideset();
408 
409  // If the sideset doesn't exist, there are no partitions to create
410  if(sideset_map.find(sideset_name) == sideset_map.end())
411  return;
412 
413  const panzer::LocalMeshSidesetInfo<LO,GO> & sideset_info = sideset_map.at(sideset_name);
414 
415  // Partitioning is not important for sidesets
416  panzer::partitioning_utilities::splitMeshInfo<LO,GO>(sideset_info, description.getWorksetSize(), partitions);
417 
418  for(auto & partition : partitions){
419  partition.sideset_name = sideset_name;
420  partition.element_block_name = element_block_name;
421  partition.cell_topology = sideset_info.cell_topology;
422  }
423 
424  } else {
425 
426  // If the element block doesn't exist, there are no partitions to create
427  if(mesh_info.element_blocks.find(element_block_name) == mesh_info.element_blocks.end())
428  return;
429 
430  // Grab the element block we're interested in
431  const panzer::LocalMeshBlockInfo<LO,GO> & block_info = mesh_info.element_blocks.at(element_block_name);
432 
434  // We only have one partition describing the entire local mesh
435  panzer::partitioning_utilities::splitMeshInfo(block_info, -1, partitions);
436  } else {
437  // We need to partition local mesh
438 
439  // TODO: This needs to be used, but first needs to be written
440  //panzer::partitioning_utilities::partitionMeshInfo<LO,GO>(block_info, description.getWorksetSize(), partitions);
441 
442  // FIXME: Until the above function is fixed, we will use this hack - this will lead to horrible partitions
443  panzer::partitioning_utilities::splitMeshInfo(block_info, description.getWorksetSize(), partitions);
444 
445  }
446 
447  for(auto & partition : partitions){
448  partition.element_block_name = element_block_name;
449  partition.cell_topology = block_info.cell_topology;
450  }
451  }
452 
453 }
454 
455 }
456 
457 
458 template
459 void
460 panzer::partitioning_utilities::setupSubLocalMeshInfo<int,int>(const panzer::LocalMeshInfoBase<int,int> & parent_info,
461  const std::vector<int> & owned_parent_cells,
463 
464 #ifndef PANZER_ORDINAL64_IS_INT
465 template
466 void
467 panzer::partitioning_utilities::setupSubLocalMeshInfo<int,panzer::Ordinal64>(const panzer::LocalMeshInfoBase<int,panzer::Ordinal64> & parent_info,
468  const std::vector<int> & owned_parent_cells,
470 #endif
471 
472 template
473 void
474 panzer::generateLocalMeshPartitions<int,int>(const panzer::LocalMeshInfo<int,int> & mesh_info,
475  const panzer::WorksetDescriptor & description,
476  std::vector<panzer::LocalMeshPartition<int,int> > & partitions);
477 
478 #ifndef PANZER_ORDINAL64_IS_INT
479 template
480 void
481 panzer::generateLocalMeshPartitions<int,panzer::Ordinal64>(const panzer::LocalMeshInfo<int,panzer::Ordinal64> & mesh_info,
482  const panzer::WorksetDescriptor & description,
483  std::vector<panzer::LocalMeshPartition<int,panzer::Ordinal64> > & partitions);
484 #endif
Backwards compatibility mode that ignores the worksetSize in the WorksetDescriptor.
std::map< std::string, LocalMeshBlockInfo< LO, GO > > element_blocks
Kokkos::View< GO * > global_cells
Teuchos::RCP< const shards::CellTopology > cell_topology
const std::string & getSideset(const int block=0) const
Get sideset name.
void setupSubLocalMeshInfo(const panzer::LocalMeshInfoBase< LO, GO > &parent_info, const std::vector< LO > &owned_parent_cells, panzer::LocalMeshInfoBase< LO, GO > &sub_info)
void splitMeshInfo(const panzer::LocalMeshInfoBase< LO, GO > &mesh_info, const int splitting_size, std::vector< panzer::LocalMeshPartition< LO, GO > > &partitions)
std::map< std::string, std::map< std::string, LocalMeshSidesetInfo< LO, GO > > > sidesets
Kokkos::View< LO * > local_cells
bool requiresPartitioning() const
Do we need to partition the local mesh prior to generating worksets.
int getWorksetSize() const
Get the requested workset size (default -2 (workset size is set elsewhere), -1 (largest possible work...
bool useSideset() const
This descriptor is for a side set.
Kokkos::View< LO *[2]> face_to_cells
Kokkos::View< double ***, PHX::Device > cell_vertices
Kokkos::View< LO ** > cell_to_faces
Kokkos::View< LO *[2]> face_to_lidx
void partitionMeshInfo(const panzer::LocalMeshInfoBase< LO, GO > &, const size_t, std::vector< panzer::LocalMeshPartition< LO, GO >> &)
void generateLocalMeshPartitions(const panzer::LocalMeshInfo< LO, GO > &mesh_info, const panzer::WorksetDescriptor &description, std::vector< panzer::LocalMeshPartition< LO, GO > > &partitions)
const std::string & getElementBlock(const int block=0) const
Get element block name.
Workset size is set to the total number of local elements in the MPI process.