FreeFOAM The Cross-Platform CFD Toolkit
scotchDecomp.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24  From scotch forum:
25 
26  By: Francois PELLEGRINI RE: Graph mapping 'strategy' string [ reply ]
27  2008-08-22 10:09 Strategy handling in Scotch is a bit tricky. In order
28  not to be confused, you must have a clear view of how they are built.
29  Here are some rules:
30 
31  1- Strategies are made up of "methods" which are combined by means of
32  "operators".
33 
34  2- A method is of the form "m{param=value,param=value,...}", where "m"
35  is a single character (this is your first error: "f" is a method name,
36  not a parameter name).
37 
38  3- There exist different sort of strategies : bipartitioning strategies,
39  mapping strategies, ordering strategies, which cannot be mixed. For
40  instance, you cannot build a bipartitioning strategy and feed it to a
41  mapping method (this is your second error).
42 
43  To use the "mapCompute" routine, you must create a mapping strategy, not
44  a bipartitioning one, and so use stratGraphMap() and not
45  stratGraphBipart(). Your mapping strategy should however be based on the
46  "recursive bipartitioning" method ("b"). For instance, a simple (and
47  hence not very efficient) mapping strategy can be :
48 
49  "b{sep=f}"
50 
51  which computes mappings with the recursive bipartitioning method "b",
52  this latter using the Fiduccia-Mattheyses method "f" to compute its
53  separators.
54 
55  If you want an exact partition (see your previous post), try
56  "b{sep=fx}".
57 
58  However, these strategies are not the most efficient, as they do not
59  make use of the multi-level framework.
60 
61  To use the multi-level framework, try for instance:
62 
63  "b{sep=m{vert=100,low=h,asc=f}x}"
64 
65  The current default mapping strategy in Scotch can be seen by using the
66  "-vs" option of program gmap. It is, to date:
67 
68  b
69  {
70  job=t,
71  map=t,
72  poli=S,
73  sep=
74  (
75  m
76  {
77  asc=b
78  {
79  bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
80  org=f{move=80,pass=-1,bal=0.005},
81  width=3
82  },
83  low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
84  type=h,
85  vert=80,
86  rat=0.8
87  }
88  | m
89  {
90  asc=b
91  {
92  bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
93  org=f{move=80,pass=-1,bal=0.005},
94  width=3
95  },
96  low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
97  type=h,
98  vert=80,
99  rat=0.8
100  }
101  )
102  }
103 
104 
105 \*---------------------------------------------------------------------------*/
106 
107 #include "scotchDecomp.H"
109 #include <OpenFOAM/floatScalar.H>
110 #include <OpenFOAM/Time.H>
112 #include <OpenFOAM/OFstream.H>
113 
114 #include "scotch.h"
115 
116 
117 // Hack: scotch generates floating point errors so need to switch of error
118 // trapping!
119 #if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
120 # define LINUX
121 #endif
122 
123 #if defined(LINUX) && defined(__GNUC__)
124 # define LINUX_GNUC
125 #endif
126 
127 #ifdef LINUX_GNUC
128 # ifndef __USE_GNU
129 # define __USE_GNU
130 # endif
131 # include <fenv.h>
132 #endif
133 
134 
135 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
136 
137 namespace Foam
138 {
139  defineTypeNameAndDebug(scotchDecomp, 0);
140 
142  (
143  decompositionMethod,
144  scotchDecomp,
145  dictionaryMesh
146  );
147 }
148 
149 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
150 
151 void Foam::scotchDecomp::check(const int retVal, const char* str)
152 {
153  if (retVal)
154  {
155  FatalErrorIn("scotchDecomp::decompose(..)")
156  << "Call to scotch routine " << str << " failed."
157  << exit(FatalError);
158  }
159 }
160 
161 
162 // Call scotch with options from dictionary.
163 Foam::label Foam::scotchDecomp::decompose
164 (
165  const List<int>& adjncy,
166  const List<int>& xadj,
167  const scalarField& cWeights,
168 
169  List<int>& finalDecomp
170 )
171 {
172  // Dump graph
173  if (decompositionDict_.found("scotchCoeffs"))
174  {
175  const dictionary& scotchCoeffs =
176  decompositionDict_.subDict("scotchCoeffs");
177 
178  if (scotchCoeffs.found("writeGraph"))
179  {
180  Switch writeGraph(scotchCoeffs.lookup("writeGraph"));
181 
182  if (writeGraph)
183  {
184  OFstream str(mesh_.time().path() / mesh_.name() + ".grf");
185 
186  Info<< "Dumping Scotch graph file to " << str.name() << endl
187  << "Use this in combination with gpart." << endl;
188 
189  label version = 0;
190  str << version << nl;
191  // Numer of vertices
192  str << xadj.size()-1 << ' ' << adjncy.size() << nl;
193  // Numbering starts from 0
194  label baseval = 0;
195  // Has weights?
196  label hasEdgeWeights = 0;
197  label hasVertexWeights = 0;
198  label numericflag = 10*hasEdgeWeights+hasVertexWeights;
199  str << baseval << ' ' << numericflag << nl;
200  for (label cellI = 0; cellI < xadj.size()-1; cellI++)
201  {
202  label start = xadj[cellI];
203  label end = xadj[cellI+1];
204  str << end-start;
205 
206  for (label i = start; i < end; i++)
207  {
208  str << ' ' << adjncy[i];
209  }
210  str << nl;
211  }
212  }
213  }
214  }
215 
216 
217  // Strategy
218  // ~~~~~~~~
219 
220  // Default.
221  SCOTCH_Strat stradat;
222  check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
223 
224  if (decompositionDict_.found("scotchCoeffs"))
225  {
226  const dictionary& scotchCoeffs =
227  decompositionDict_.subDict("scotchCoeffs");
228 
229 
230  string strategy;
231  if (scotchCoeffs.readIfPresent("strategy", strategy))
232  {
233  if (debug)
234  {
235  Info<< "scotchDecomp : Using strategy " << strategy << endl;
236  }
237  SCOTCH_stratGraphMap(&stradat, strategy.c_str());
238  //fprintf(stdout, "S\tStrat=");
239  //SCOTCH_stratSave(&stradat, stdout);
240  //fprintf(stdout, "\n");
241  }
242  }
243 
244 
245  // Graph
246  // ~~~~~
247 
248  List<int> velotab;
249 
250 
251  // Check for externally provided cellweights and if so initialise weights
252  scalar minWeights = gMin(cWeights);
253  if (cWeights.size() > 0)
254  {
255  if (minWeights <= 0)
256  {
257  WarningIn
258  (
259  "scotchDecomp::decompose"
260  "(const pointField&, const scalarField&)"
261  ) << "Illegal minimum weight " << minWeights
262  << endl;
263  }
264 
265  if (cWeights.size() != xadj.size()-1)
266  {
268  (
269  "scotchDecomp::decompose"
270  "(const pointField&, const scalarField&)"
271  ) << "Number of cell weights " << cWeights.size()
272  << " does not equal number of cells " << xadj.size()-1
273  << exit(FatalError);
274  }
275 
276  // Convert to integers.
277  velotab.setSize(cWeights.size());
278  forAll(velotab, i)
279  {
280  velotab[i] = int(cWeights[i]/minWeights);
281  }
282  }
283 
284 
285 
286  SCOTCH_Graph grafdat;
287  check(SCOTCH_graphInit(&grafdat), "SCOTCH_graphInit");
288  check
289  (
290  SCOTCH_graphBuild
291  (
292  &grafdat,
293  0, // baseval, c-style numbering
294  xadj.size()-1, // vertnbr, nCells
295  xadj.begin(), // verttab, start index per cell into adjncy
296  &xadj[1], // vendtab, end index ,,
297  velotab.begin(), // velotab, vertex weights
298  NULL, // vlbltab
299  adjncy.size(), // edgenbr, number of arcs
300  adjncy.begin(), // edgetab
301  NULL // edlotab, edge weights
302  ),
303  "SCOTCH_graphBuild"
304  );
305  check(SCOTCH_graphCheck(&grafdat), "SCOTCH_graphCheck");
306 
307 
308  // Architecture
309  // ~~~~~~~~~~~~
310  // (fully connected network topology since using switch)
311 
312  SCOTCH_Arch archdat;
313  check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
314 
315  List<label> processorWeights;
316  if (decompositionDict_.found("scotchCoeffs"))
317  {
318  const dictionary& scotchCoeffs =
319  decompositionDict_.subDict("scotchCoeffs");
320 
321  scotchCoeffs.readIfPresent("processorWeights", processorWeights);
322  }
323  if (processorWeights.size())
324  {
325  if (debug)
326  {
327  Info<< "scotchDecomp : Using procesor weights " << processorWeights
328  << endl;
329  }
330  check
331  (
332  SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
333  "SCOTCH_archCmpltw"
334  );
335  }
336  else
337  {
338  check
339  (
340  SCOTCH_archCmplt(&archdat, nProcessors_),
341  "SCOTCH_archCmplt"
342  );
343  }
344 
345 
346  //SCOTCH_Mapping mapdat;
347  //SCOTCH_graphMapInit(&grafdat, &mapdat, &archdat, NULL);
348  //SCOTCH_graphMapCompute(&grafdat, &mapdat, &stradat); /* Perform mapping */
349  //SCOTCH_graphMapExit(&grafdat, &mapdat);
350 
351 
352  // Hack:switch off fpu error trapping
353 # ifdef LINUX_GNUC
354  int oldExcepts = fedisableexcept
355  (
356  FE_DIVBYZERO
357  | FE_INVALID
358  | FE_OVERFLOW
359  );
360 # endif
361 
362  finalDecomp.setSize(xadj.size()-1);
363  finalDecomp = 0;
364  check
365  (
366  SCOTCH_graphMap
367  (
368  &grafdat,
369  &archdat,
370  &stradat, // const SCOTCH_Strat *
371  finalDecomp.begin() // parttab
372  ),
373  "SCOTCH_graphMap"
374  );
375 
376 # ifdef LINUX_GNUC
377  feenableexcept(oldExcepts);
378 # endif
379 
380 
381 
382  //finalDecomp.setSize(xadj.size()-1);
383  //check
384  //(
385  // SCOTCH_graphPart
386  // (
387  // &grafdat,
388  // nProcessors_, // partnbr
389  // &stradat, // const SCOTCH_Strat *
390  // finalDecomp.begin() // parttab
391  // ),
392  // "SCOTCH_graphPart"
393  //);
394 
395  // Release storage for graph
396  SCOTCH_graphExit(&grafdat);
397  // Release storage for strategy
398  SCOTCH_stratExit(&stradat);
399  // Release storage for network topology
400  SCOTCH_archExit(&archdat);
401 
402  return 0;
403 }
404 
405 
406 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
407 
408 Foam::scotchDecomp::scotchDecomp
409 (
410  const dictionary& decompositionDict,
411  const polyMesh& mesh
412 )
413 :
414  decompositionMethod(decompositionDict),
415  mesh_(mesh)
416 {}
417 
418 
419 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
420 
421 Foam::labelList Foam::scotchDecomp::decompose
422 (
423  const pointField& points,
424  const scalarField& pointWeights
425 )
426 {
427  if (points.size() != mesh_.nCells())
428  {
430  (
431  "scotchDecomp::decompose(const pointField&, const scalarField&)"
432  )
433  << "Can use this decomposition method only for the whole mesh"
434  << endl
435  << "and supply one coordinate (cellCentre) for every cell." << endl
436  << "The number of coordinates " << points.size() << endl
437  << "The number of cells in the mesh " << mesh_.nCells()
438  << exit(FatalError);
439  }
440 
441  // Make Metis CSR (Compressed Storage Format) storage
442  // adjncy : contains neighbours (= edges in graph)
443  // xadj(celli) : start of information in adjncy for celli
444  List<int> adjncy;
445  List<int> xadj;
446  calcCSR(mesh_, adjncy, xadj);
447 
448  // Decompose using default weights
449  List<int> finalDecomp;
450  decompose(adjncy, xadj, pointWeights, finalDecomp);
451 
452  // Copy back to labelList
453  labelList decomp(finalDecomp.size());
454  forAll(decomp, i)
455  {
456  decomp[i] = finalDecomp[i];
457  }
458  return decomp;
459 }
460 
461 
462 Foam::labelList Foam::scotchDecomp::decompose
463 (
464  const labelList& agglom,
465  const pointField& agglomPoints,
466  const scalarField& pointWeights
467 )
468 {
469  if (agglom.size() != mesh_.nCells())
470  {
472  (
473  "parMetisDecomp::decompose(const labelList&, const pointField&)"
474  ) << "Size of cell-to-coarse map " << agglom.size()
475  << " differs from number of cells in mesh " << mesh_.nCells()
476  << exit(FatalError);
477  }
478 
479  // Make Metis CSR (Compressed Storage Format) storage
480  // adjncy : contains neighbours (= edges in graph)
481  // xadj(celli) : start of information in adjncy for celli
482  List<int> adjncy;
483  List<int> xadj;
484  {
485  // Get cellCells on coarse mesh.
486  labelListList cellCells;
487  calcCellCells
488  (
489  mesh_,
490  agglom,
491  agglomPoints.size(),
492  cellCells
493  );
494 
495  calcCSR(cellCells, adjncy, xadj);
496  }
497 
498  // Decompose using weights
499  List<int> finalDecomp;
500  decompose(adjncy, xadj, pointWeights, finalDecomp);
501 
502  // Rework back into decomposition for original mesh_
503  labelList fineDistribution(agglom.size());
504 
505  forAll(fineDistribution, i)
506  {
507  fineDistribution[i] = finalDecomp[agglom[i]];
508  }
509 
510  return fineDistribution;
511 }
512 
513 
514 Foam::labelList Foam::scotchDecomp::decompose
515 (
516  const labelListList& globalCellCells,
517  const pointField& cellCentres,
518  const scalarField& cWeights
519 )
520 {
521  if (cellCentres.size() != globalCellCells.size())
522  {
524  (
525  "scotchDecomp::decompose"
526  "(const labelListList&, const pointField&, const scalarField&)"
527  ) << "Inconsistent number of cells (" << globalCellCells.size()
528  << ") and number of cell centres (" << cellCentres.size()
529  << ")." << exit(FatalError);
530  }
531 
532 
533  // Make Metis CSR (Compressed Storage Format) storage
534  // adjncy : contains neighbours (= edges in graph)
535  // xadj(celli) : start of information in adjncy for celli
536 
537  List<int> adjncy;
538  List<int> xadj;
539  calcCSR(globalCellCells, adjncy, xadj);
540 
541  // Decompose using weights
542  List<int> finalDecomp;
543  decompose(adjncy, xadj, cWeights, finalDecomp);
544 
545  // Copy back to labelList
546  labelList decomp(finalDecomp.size());
547  forAll(decomp, i)
548  {
549  decomp[i] = finalDecomp[i];
550  }
551  return decomp;
552 }
553 
554 
556 (
557  const polyMesh& mesh,
558  List<int>& adjncy,
559  List<int>& xadj
560 )
561 {
562  // Make Metis CSR (Compressed Storage Format) storage
563  // adjncy : contains neighbours (= edges in graph)
564  // xadj(celli) : start of information in adjncy for celli
565 
566  xadj.setSize(mesh.nCells()+1);
567 
568  // Initialise the number of internal faces of the cells to twice the
569  // number of internal faces
570  label nInternalFaces = 2*mesh.nInternalFaces();
571 
572  // Check the boundary for coupled patches and add to the number of
573  // internal faces
574  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
575 
576  forAll(pbm, patchi)
577  {
578  if (isA<cyclicPolyPatch>(pbm[patchi]))
579  {
580  nInternalFaces += pbm[patchi].size();
581  }
582  }
583 
584  // Create the adjncy array the size of the total number of internal and
585  // coupled faces
586  adjncy.setSize(nInternalFaces);
587 
588  // Fill in xadj
589  // ~~~~~~~~~~~~
590  label freeAdj = 0;
591 
592  for (label cellI = 0; cellI < mesh.nCells(); cellI++)
593  {
594  xadj[cellI] = freeAdj;
595 
596  const labelList& cFaces = mesh.cells()[cellI];
597 
598  forAll(cFaces, i)
599  {
600  label faceI = cFaces[i];
601 
602  if
603  (
604  mesh.isInternalFace(faceI)
605  || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
606  )
607  {
608  freeAdj++;
609  }
610  }
611  }
612  xadj[mesh.nCells()] = freeAdj;
613 
614 
615  // Fill in adjncy
616  // ~~~~~~~~~~~~~~
617 
618  labelList nFacesPerCell(mesh.nCells(), 0);
619 
620  // Internal faces
621  for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
622  {
623  label own = mesh.faceOwner()[faceI];
624  label nei = mesh.faceNeighbour()[faceI];
625 
626  adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
627  adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
628  }
629 
630  // Coupled faces. Only cyclics done.
631  forAll(pbm, patchi)
632  {
633  if (isA<cyclicPolyPatch>(pbm[patchi]))
634  {
635  const unallocLabelList& faceCells = pbm[patchi].faceCells();
636 
637  label sizeby2 = faceCells.size()/2;
638 
639  for (label facei=0; facei<sizeby2; facei++)
640  {
641  label own = faceCells[facei];
642  label nei = faceCells[facei + sizeby2];
643 
644  adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
645  adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
646  }
647  }
648  }
649 }
650 
651 
652 // From cell-cell connections to Metis format (like CompactListList)
654 (
655  const labelListList& cellCells,
656  List<int>& adjncy,
657  List<int>& xadj
658 )
659 {
660  // Count number of internal faces
661  label nConnections = 0;
662 
663  forAll(cellCells, coarseI)
664  {
665  nConnections += cellCells[coarseI].size();
666  }
667 
668  // Create the adjncy array as twice the size of the total number of
669  // internal faces
670  adjncy.setSize(nConnections);
671 
672  xadj.setSize(cellCells.size()+1);
673 
674 
675  // Fill in xadj
676  // ~~~~~~~~~~~~
677  label freeAdj = 0;
678 
679  forAll(cellCells, coarseI)
680  {
681  xadj[coarseI] = freeAdj;
682 
683  const labelList& cCells = cellCells[coarseI];
684 
685  forAll(cCells, i)
686  {
687  adjncy[freeAdj++] = cCells[i];
688  }
689  }
690  xadj[cellCells.size()] = freeAdj;
691 }
692 
693 
694 
695 // ************************ vim: set sw=4 sts=4 et: ************************ //