GAMGAgglomeration.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
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 \*---------------------------------------------------------------------------*/
25 
26 #include "GAMGAgglomeration.H"
27 #include "lduMesh.H"
28 #include "lduMatrix.H"
29 #include "Time.H"
30 #include "GAMGInterface.H"
31 #include "GAMGProcAgglomeration.H"
32 #include "pairGAMGAgglomeration.H"
33 #include "IOmanip.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39  defineTypeNameAndDebug(GAMGAgglomeration, 0);
40  defineRunTimeSelectionTable(GAMGAgglomeration, lduMesh);
41  defineRunTimeSelectionTable(GAMGAgglomeration, lduMatrix);
42  defineRunTimeSelectionTable(GAMGAgglomeration, geometry);
43 }
44 
45 
46 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
47 
49 {
50  nCells_.setSize(nCreatedLevels);
51  restrictAddressing_.setSize(nCreatedLevels);
52  nFaces_.setSize(nCreatedLevels);
53  faceRestrictAddressing_.setSize(nCreatedLevels);
54  faceFlipMap_.setSize(nCreatedLevels);
55  nPatchFaces_.setSize(nCreatedLevels);
56  patchFaceRestrictAddressing_.setSize(nCreatedLevels);
57  meshLevels_.setSize(nCreatedLevels);
58 
59  // Have procCommunicator_ always, even if not procAgglomerating
60  procCommunicator_.setSize(nCreatedLevels + 1);
61  if (processorAgglomerate())
62  {
63  procAgglomMap_.setSize(nCreatedLevels);
64  agglomProcIDs_.setSize(nCreatedLevels);
65  procCellOffsets_.setSize(nCreatedLevels);
66  procFaceMap_.setSize(nCreatedLevels);
67  procBoundaryMap_.setSize(nCreatedLevels);
68  procBoundaryFaceMap_.setSize(nCreatedLevels);
69 
70  procAgglomeratorPtr_().agglomerate();
71 
72 
73  }
74 
75  // Print a bit
76  if (processorAgglomerate() && debug)
77  {
78  Info<< "GAMGAgglomeration:" << nl
79  << " local agglomerator : " << type() << nl;
80  if (processorAgglomerate())
81  {
82  Info<< " processor agglomerator : "
83  << procAgglomeratorPtr_().type() << nl
84  << nl;
85  }
86 
87  Info<< setw(36) << "nCells"
88  << setw(20) << "nFaces/nCells"
89  << setw(20) << "nInterfaces"
90  << setw(20) << "nIntFaces/nCells"
91  << setw(12) << "profile"
92  << nl
93  << setw(8) << "Level"
94  << setw(8) << "nProcs"
95  << " "
96  << setw(8) << "avg"
97  << setw(8) << "max"
98  << " "
99  << setw(8) << "avg"
100  << setw(8) << "max"
101  << " "
102  << setw(8) << "avg"
103  << setw(8) << "max"
104  << " "
105  << setw(8) << "avg"
106  << setw(8) << "max"
107  //<< " "
108  << setw(12) << "avg"
109  << nl
110  << setw(8) << "-----"
111  << setw(8) << "------"
112  << " "
113  << setw(8) << "---"
114  << setw(8) << "---"
115  << " "
116  << setw(8) << "---"
117  << setw(8) << "---"
118  << " "
119  << setw(8) << "---"
120  << setw(8) << "---"
121  << " "
122  << setw(8) << "---"
123  << setw(8) << "---"
124  //<< " "
125  << setw(12) << "---"
126  //<< " "
127  << nl;
128 
129  for (label levelI = 0; levelI <= size(); levelI++)
130  {
131  label nProcs = 0;
132  label nCells = 0;
133  scalar faceCellRatio = 0;
134  label nInterfaces = 0;
135  label nIntFaces = 0;
136  scalar ratio = 0.0;
137  scalar profile = 0.0;
138 
139  if (hasMeshLevel(levelI))
140  {
141  nProcs = 1;
142 
143  const lduMesh& fineMesh = meshLevel(levelI);
144  nCells = fineMesh.lduAddr().size();
145  faceCellRatio =
146  scalar(fineMesh.lduAddr().lowerAddr().size())/nCells;
147 
148  const lduInterfacePtrsList interfaces =
149  fineMesh.interfaces();
150  forAll(interfaces, i)
151  {
152  if (interfaces.set(i))
153  {
154  nInterfaces++;
155  nIntFaces += interfaces[i].faceCells().size();
156  }
157  }
158  ratio = scalar(nIntFaces)/nCells;
159 
160  profile = fineMesh.lduAddr().band().second();
161  }
162 
163  label totNprocs = returnReduce(nProcs, sumOp<label>());
164 
165  label maxNCells = returnReduce(nCells, maxOp<label>());
166  label totNCells = returnReduce(nCells, sumOp<label>());
167 
168  scalar maxFaceCellRatio =
169  returnReduce(faceCellRatio, maxOp<scalar>());
170  scalar totFaceCellRatio =
171  returnReduce(faceCellRatio, sumOp<scalar>());
172 
173  label maxNInt = returnReduce(nInterfaces, maxOp<label>());
174  label totNInt = returnReduce(nInterfaces, sumOp<label>());
175 
176  scalar maxRatio = returnReduce(ratio, maxOp<scalar>());
177  scalar totRatio = returnReduce(ratio, sumOp<scalar>());
178 
179  scalar totProfile = returnReduce(profile, sumOp<scalar>());
180 
181  int oldPrecision = Info().precision(4);
182 
183  Info<< setw(8) << levelI
184  << setw(8) << totNprocs
185  << " "
186  << setw(8) << totNCells/totNprocs
187  << setw(8) << maxNCells
188  << " "
189  << setw(8) << totFaceCellRatio/totNprocs
190  << setw(8) << maxFaceCellRatio
191  << " "
192  << setw(8) << scalar(totNInt)/totNprocs
193  << setw(8) << maxNInt
194  << " "
195  << setw(8) << totRatio/totNprocs
196  << setw(8) << maxRatio
197  << setw(12) << totProfile/totNprocs
198  << nl;
199 
200  Info().precision(oldPrecision);
201  }
202  Info<< endl;
203  }
204 }
205 
206 
208 (
209  const label nFineCells,
210  const label nCoarseCells
211 ) const
212 {
213  const label nTotalCoarseCells = returnReduce(nCoarseCells, sumOp<label>());
214  if (nTotalCoarseCells < Pstream::nProcs()*nCellsInCoarsestLevel_)
215  {
216  return false;
217  }
218  else
219  {
220  const label nTotalFineCells = returnReduce(nFineCells, sumOp<label>());
221  return nTotalCoarseCells < nTotalFineCells;
222  }
223 }
224 
225 
226 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
227 
229 (
230  const lduMesh& mesh,
231  const dictionary& controlDict
232 )
233 :
235 
236  maxLevels_(50),
237 
238  nCellsInCoarsestLevel_
239  (
240  controlDict.lookupOrDefault<label>("nCellsInCoarsestLevel", 10)
241  ),
242  meshInterfaces_(mesh.interfaces()),
243  procAgglomeratorPtr_
244  (
245  (
246  (UPstream::nProcs(mesh.comm()) > 1)
247  && controlDict.found("processorAgglomerator")
248  )
250  (
251  controlDict.lookup("processorAgglomerator"),
252  *this,
254  )
256  ),
257 
258  nCells_(maxLevels_),
259  restrictAddressing_(maxLevels_),
260  nFaces_(maxLevels_),
261  faceRestrictAddressing_(maxLevels_),
262  faceFlipMap_(maxLevels_),
263  nPatchFaces_(maxLevels_),
264  patchFaceRestrictAddressing_(maxLevels_),
265 
266  meshLevels_(maxLevels_)
267 {
268  procCommunicator_.setSize(maxLevels_ + 1, -1);
269  if (processorAgglomerate())
270  {
271  procAgglomMap_.setSize(maxLevels_);
272  agglomProcIDs_.setSize(maxLevels_);
273  procCellOffsets_.setSize(maxLevels_);
274  procFaceMap_.setSize(maxLevels_);
275  procBoundaryMap_.setSize(maxLevels_);
276  procBoundaryFaceMap_.setSize(maxLevels_);
277  }
278 }
279 
280 
282 (
283  const lduMesh& mesh,
284  const dictionary& controlDict
285 )
286 {
287  if
288  (
290  (
291  GAMGAgglomeration::typeName
292  )
293  )
294  {
295  const word agglomeratorType
296  (
297  controlDict.lookupOrDefault<word>("agglomerator", "faceAreaPair")
298  );
299 
300  libs.open
301  (
302  controlDict,
303  "geometricGAMGAgglomerationLibs",
304  lduMeshConstructorTablePtr_
305  );
306 
307  lduMeshConstructorTable::iterator cstrIter =
308  lduMeshConstructorTablePtr_->find(agglomeratorType);
309 
310  if (cstrIter == lduMeshConstructorTablePtr_->end())
311  {
313  << "Unknown GAMGAgglomeration type "
314  << agglomeratorType << ".\n"
315  << "Valid matrix GAMGAgglomeration types are :"
316  << lduMatrixConstructorTablePtr_->sortedToc() << endl
317  << "Valid geometric GAMGAgglomeration types are :"
318  << lduMeshConstructorTablePtr_->sortedToc()
319  << exit(FatalError);
320  }
321 
322  return store(cstrIter()(mesh, controlDict).ptr());
323  }
324  else
325  {
326  return mesh.thisDb().lookupObject<GAMGAgglomeration>
327  (
328  GAMGAgglomeration::typeName
329  );
330  }
331 }
332 
333 
335 (
336  const lduMatrix& matrix,
337  const dictionary& controlDict
338 )
339 {
340  const lduMesh& mesh = matrix.mesh();
341 
342  if
343  (
345  (
346  GAMGAgglomeration::typeName
347  )
348  )
349  {
350  const word agglomeratorType
351  (
352  controlDict.lookupOrDefault<word>("agglomerator", "faceAreaPair")
353  );
354 
355  libs.open
356  (
357  controlDict,
358  "algebraicGAMGAgglomerationLibs",
359  lduMatrixConstructorTablePtr_
360  );
361 
362  if
363  (
364  !lduMatrixConstructorTablePtr_
365  || !lduMatrixConstructorTablePtr_->found(agglomeratorType)
366  )
367  {
368  return New(mesh, controlDict);
369  }
370  else
371  {
372  lduMatrixConstructorTable::iterator cstrIter =
373  lduMatrixConstructorTablePtr_->find(agglomeratorType);
374 
375  return store(cstrIter()(matrix, controlDict).ptr());
376  }
377  }
378  else
379  {
380  return mesh.thisDb().lookupObject<GAMGAgglomeration>
381  (
382  GAMGAgglomeration::typeName
383  );
384  }
385 }
386 
387 
389 (
390  const lduMesh& mesh,
391  const scalarField& cellVolumes,
392  const vectorField& faceAreas,
393  const dictionary& controlDict
394 )
395 {
396  const word agglomeratorType
397  (
398  controlDict.lookupOrDefault<word>("agglomerator", "faceAreaPair")
399  );
400 
401  libs.open
402  (
403  controlDict,
404  "geometricGAMGAgglomerationLibs",
405  geometryConstructorTablePtr_
406  );
407 
408  geometryConstructorTable::iterator cstrIter =
409  geometryConstructorTablePtr_->find(agglomeratorType);
410 
411  if (cstrIter == geometryConstructorTablePtr_->end())
412  {
414  << "Unknown GAMGAgglomeration type "
415  << agglomeratorType << ".\n"
416  << "Valid geometric GAMGAgglomeration types are :"
417  << geometryConstructorTablePtr_->sortedToc()
418  << exit(FatalError);
419  }
420 
422  (
423  cstrIter()
424  (
425  mesh,
426  cellVolumes,
427  faceAreas,
429  )
430  );
431 }
432 
433 
434 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
435 
437 {}
438 
439 
440 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
441 
443 (
444  const label i
445 ) const
446 {
447  if (i == 0)
448  {
449  return mesh_;
450  }
451  else
452  {
453  return meshLevels_[i - 1];
454  }
455 }
456 
457 
459 {
460  if (i == 0)
461  {
462  return true;
463  }
464  else
465  {
466  return meshLevels_.set(i - 1);
467  }
468 }
469 
470 
472 (
473  const label i
474 ) const
475 {
476  if (i == 0)
477  {
478  return meshInterfaces_;
479  }
480  else
481  {
482  return meshLevels_[i - 1].rawInterfaces();
483  }
484 }
485 
486 
488 {
489  if (hasMeshLevel(i))
490  {
491  meshLevels_.set(i - 1, nullptr);
492 
493  if (i < nCells_.size())
494  {
495  nCells_[i] = -555;
496  restrictAddressing_.set(i, nullptr);
497  nFaces_[i] = -666;
498  faceRestrictAddressing_.set(i, nullptr);
499  faceFlipMap_.set(i, nullptr);
500  nPatchFaces_.set(i, nullptr);
501  patchFaceRestrictAddressing_.set(i, nullptr);
502  }
503  }
504 }
505 
506 
508 (
509  const label leveli
510 ) const
511 {
512  return procAgglomMap_[leveli];
513 }
514 
515 
517 (
518  const label leveli
519 ) const
520 {
521  return agglomProcIDs_[leveli];
522 }
523 
524 
526 {
527  return procCommunicator_[leveli] != -1;
528 }
529 
530 
532 {
533  return procCommunicator_[leveli];
534 }
535 
536 
538 (
539  const label leveli
540 ) const
541 {
542  return procCellOffsets_[leveli];
543 }
544 
545 
547 (
548  const label leveli
549 ) const
550 {
551  return procFaceMap_[leveli];
552 }
553 
554 
556 (
557  const label leveli
558 ) const
559 {
560  return procBoundaryMap_[leveli];
561 }
562 
563 
565 (
566  const label leveli
567 ) const
568 {
569  return procBoundaryFaceMap_[leveli];
570 }
571 
572 
574 (
575  labelList& newRestrict,
576  label& nNewCoarse,
577  const lduAddressing& fineAddressing,
578  const labelUList& restrict,
579  const label nCoarse
580 )
581 {
582  if (fineAddressing.size() != restrict.size())
583  {
585  << "nCells:" << fineAddressing.size()
586  << " agglom:" << restrict.size()
587  << abort(FatalError);
588  }
589 
590  // Seed (master) for every region
591  labelList master(identity(fineAddressing.size()));
592 
593  // Now loop and transport master through region
594  const labelUList& lower = fineAddressing.lowerAddr();
595  const labelUList& upper = fineAddressing.upperAddr();
596 
597  while (true)
598  {
599  label nChanged = 0;
600 
601  forAll(lower, facei)
602  {
603  label own = lower[facei];
604  label nei = upper[facei];
605 
606  if (restrict[own] == restrict[nei])
607  {
608  // coarse-mesh-internal face
609 
610  if (master[own] < master[nei])
611  {
612  master[nei] = master[own];
613  nChanged++;
614  }
615  else if (master[own] > master[nei])
616  {
617  master[own] = master[nei];
618  nChanged++;
619  }
620  }
621  }
622 
623  reduce(nChanged, sumOp<label>());
624 
625  if (nChanged == 0)
626  {
627  break;
628  }
629  }
630 
631 
632  // Count number of regions/masters per coarse cell
633  labelListList coarseToMasters(nCoarse);
634  nNewCoarse = 0;
635  forAll(restrict, celli)
636  {
637  labelList& masters = coarseToMasters[restrict[celli]];
638 
639  if (findIndex(masters, master[celli]) == -1)
640  {
641  masters.append(master[celli]);
642  nNewCoarse++;
643  }
644  }
645 
646  if (nNewCoarse > nCoarse)
647  {
648  // WarningInFunction
649  // << "Have " << nCoarse
650  // << " agglomerated cells but " << nNewCoarse
651  // << " disconnected regions" << endl;
652 
653  // Keep coarseToMasters[0] the original coarse, allocate new ones
654  // for the others
655  labelListList coarseToNewCoarse(coarseToMasters.size());
656 
657  nNewCoarse = nCoarse;
658 
659  forAll(coarseToMasters, coarseI)
660  {
661  const labelList& masters = coarseToMasters[coarseI];
662 
663  labelList& newCoarse = coarseToNewCoarse[coarseI];
664  newCoarse.setSize(masters.size());
665  newCoarse[0] = coarseI;
666  for (label i=1; i<newCoarse.size(); i++)
667  {
668  newCoarse[i] = nNewCoarse++;
669  }
670  }
671 
672  newRestrict.setSize(fineAddressing.size());
673  forAll(restrict, celli)
674  {
675  label coarseI = restrict[celli];
676 
677  label index = findIndex(coarseToMasters[coarseI], master[celli]);
678  newRestrict[celli] = coarseToNewCoarse[coarseI][index];
679  }
680 
681  return false;
682  }
683  else
684  {
685  return true;
686  }
687 }
688 
689 
690 // ************************************************************************* //
virtual lduInterfacePtrsList interfaces() const =0
Return a list of pointers for each patch.
bool hasProcMesh(const label fineLeveli) const
Check that level has combined mesh.
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:663
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
GAMGAgglomeration(const lduMesh &mesh, const dictionary &controlDict)
Construct given mesh and controls.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
runTime controlDict().lookup("adjustTimeStep") >> adjustTimeStep
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
const lduMesh & meshLevel(const label leveli) const
Return LDU mesh of given level.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const labelListList & boundaryMap(const label fineLeveli) const
Mapping from processor to procMesh boundary.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
~GAMGAgglomeration()
Destructor.
Abstract base class for meshes which provide LDU addressing for the construction of lduMatrix and LDU...
Definition: lduMesh.H:59
bool foundObject(const word &name) const
Is the named Type found?
static const GAMGAgglomeration & New(const lduMesh &mesh, const dictionary &controlDict)
Return the selected geometric agglomerator.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
virtual label comm() const =0
Return communicator used for parallel communication.
fvMesh & mesh
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type.
virtual const labelUList & lowerAddr() const =0
Return lower addressing.
static bool checkRestriction(labelList &newRestrict, label &nNewCoarse, const lduAddressing &fineAddressing, const labelUList &restrict, const label nCoarse)
Given restriction determines if coarse cells are connected.
void clearLevel(const label leveli)
dlLibraryTable libs
Table of loaded dynamic libraries.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
static autoPtr< GAMGProcAgglomeration > New(const word &type, GAMGAgglomeration &agglom, const dictionary &controlDict)
Return the selected agglomerator.
label procCommunicator(const label fineLeveli) const
Communicator for current level or -1.
Templated abstract base-class for optional mesh objects used to automate their allocation to the mesh...
Definition: MeshObject.H:87
bool continueAgglomerating(const label nCells, const label nCoarseCells) const
Check the need for further agglomeration.
A class for handling words, derived from string.
Definition: word.H:59
virtual const labelUList & upperAddr() const =0
Return upper addressing.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
const labelList & agglomProcIDs(const label fineLeveli) const
Set of processors to agglomerate. Element 0 is the.
virtual const lduAddressing & lduAddr() const =0
Return ldu addressing.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
Tuple2< label, scalar > band() const
Calculate bandwidth and profile of addressing.
Istream and Ostream manipulators taking arguments.
static const char nl
Definition: Ostream.H:260
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
bool open(const fileName &libName, const bool verbose=true)
Open the named library, optionally with warnings if problems occur.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
void setSize(const label)
Reset size of List.
Definition: List.C:281
lduMatrix is a general matrix class in which the coefficients are stored as three arrays...
Definition: lduMatrix.H:79
const labelList & cellOffsets(const label fineLeveli) const
Mapping from processor to procMesh cells.
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
const labelListListList & boundaryFaceMap(const label fineLeveli) const
Mapping from processor to procMesh boundary face.
const labelList & procAgglomMap(const label fineLeveli) const
Mapping from processor to agglomerated processor (global, all.
virtual const objectRegistry & thisDb() const
Return the object registry.
Definition: lduMesh.C:40
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
const lduInterfacePtrsList & interfaceLevel(const label leveli) const
Return LDU interface addressing of given level.
messageStream Info
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
The class contains the addressing required by the lduMatrix: upper, lower and losort.
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
Geometric agglomerated algebraic multigrid agglomeration class.
bool hasMeshLevel(const label leveli) const
Do we have mesh for given level?
const lduMesh & mesh() const
Return the LDU mesh from which the addressing is obtained.
Definition: lduMatrix.H:544
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:864
label size() const
Return number of equations.
void compactLevels(const label nCreatedLevels)
Shrink the number of levels to that specified.
const labelListList & faceMap(const label fineLeveli) const
Mapping from processor to procMesh face.