MGridGenGAMGAgglomeration.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 
27 #include "fvMesh.H"
29 #include "processorLduInterface.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
36 
38  (
41  lduMesh
42  );
43 }
44 
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 void Foam::MGridGenGAMGAgglomeration::swap
49 (
50  const lduInterfacePtrsList& interfaces,
51  const labelUList& cellValues,
52  PtrList<labelList>& nbrValues
53 ) const
54 {
55  // Initialise transfer of restrict addressing on the interface
56  forAll(interfaces, inti)
57  {
58  if (interfaces.set(inti))
59  {
60  interfaces[inti].initInternalFieldTransfer
61  (
63  cellValues
64  );
65  }
66  }
67 
68  if (Pstream::parRun())
69  {
71  }
72 
73  // Get the interface agglomeration
74  nbrValues.setSize(interfaces.size());
75  forAll(interfaces, inti)
76  {
77  if (interfaces.set(inti))
78  {
79  nbrValues.set
80  (
81  inti,
82  new labelList
83  (
84  interfaces[inti].internalFieldTransfer
85  (
87  cellValues
88  )
89  )
90  );
91  }
92  }
93 }
94 
95 
96 void Foam::MGridGenGAMGAgglomeration::getNbrAgglom
97 (
98  const lduAddressing& addr,
99  const lduInterfacePtrsList& interfaces,
100  const PtrList<labelList>& nbrGlobalAgglom,
101  labelList& cellToNbrAgglom
102 ) const
103 {
104  cellToNbrAgglom.setSize(addr.size());
105  cellToNbrAgglom = -1;
106 
107  forAll(interfaces, inti)
108  {
109  if (interfaces.set(inti))
110  {
111  if (isA<processorLduInterface>(interfaces[inti]))
112  {
113  const processorLduInterface& pldui =
114  refCast<const processorLduInterface>(interfaces[inti]);
115 
116  if (pldui.myProcNo() > pldui.neighbProcNo())
117  {
118  const labelUList& faceCells =
119  interfaces[inti].faceCells();
120  const labelList& nbrData = nbrGlobalAgglom[inti];
121 
122  forAll(faceCells, i)
123  {
124  cellToNbrAgglom[faceCells[i]] = nbrData[i];
125  }
126  }
127  }
128  }
129  }
130 }
131 
132 
133 void Foam::MGridGenGAMGAgglomeration::detectSharedFaces
134 (
135  const lduMesh& mesh,
136  const labelList& value,
137  labelHashSet& sharedFaces
138 ) const
139 {
140  const lduAddressing& addr = mesh.lduAddr();
141  const labelUList& lower = addr.lowerAddr();
142  const labelUList& upper = addr.upperAddr();
143 
144  sharedFaces.clear();
145  sharedFaces.resize(addr.lowerAddr().size()/100);
146 
147  // Detect any faces in between same value
148  forAll(lower, facei)
149  {
150  label lowerData = value[lower[facei]];
151  label upperData = value[upper[facei]];
152 
153  if (lowerData != -1 && lowerData == upperData)
154  {
155  sharedFaces.insert(facei);
156  }
157  }
158 }
159 
160 
161 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
162 
164 (
165  const lduMesh& mesh,
166  const dictionary& controlDict
167 )
168 :
170  fvMesh_(refCast<const fvMesh>(mesh))
171 {
172  // Min, max size of agglomerated cells
173  label minSize(controlDict.lookup<label>("minSize"));
174  label maxSize(controlDict.lookup<label>("maxSize"));
175 
176  // Number of iterations applied to improve agglomeration consistency across
177  // processor boundaries
178  label nProcConsistencyIter
179  (
180  controlDict.lookup<label>("nProcConsistencyIter")
181  );
182 
183  // Start geometric agglomeration from the cell volumes and areas of the mesh
184  scalarField* VPtr = const_cast<scalarField*>(&fvMesh_.cellVolumes());
185 
186  scalarField magFaceAreas(sqrt(3.0)*fvMesh_.magFaceAreas());
187  SubField<scalar> magSf(magFaceAreas, fvMesh_.nInternalFaces());
188 
189  scalarField* magSfPtr = const_cast<scalarField*>
190  (
191  &magSf.operator const scalarField&()
192  );
193 
194  // Create the boundary area cell field
195  scalarField* magSbPtr(new scalarField(fvMesh_.nCells(), 0));
196 
197  {
198  scalarField& magSb = *magSbPtr;
199 
200  const labelList& own = fvMesh_.faceOwner();
201  const vectorField& Sf = fvMesh_.faceAreas();
202 
203  forAll(Sf, facei)
204  {
205  if (!fvMesh_.isInternalFace(facei))
206  {
207  magSb[own[facei]] += mag(Sf[facei]);
208  }
209  }
210  }
211 
212  // Agglomerate until the required number of cells in the coarsest level
213  // is reached
214 
215  label nCreatedLevels = 0;
216 
217  while (nCreatedLevels < maxLevels_ - 1)
218  {
219  label nCoarseCells = -1;
220 
221  tmp<labelField> finalAgglomPtr = agglomerate
222  (
223  nCoarseCells,
224  minSize,
225  maxSize,
226  meshLevel(nCreatedLevels).lduAddr(),
227  *VPtr,
228  *magSfPtr,
229  *magSbPtr
230  );
231 
232  // Adjust weights only
233  for (int i=0; i<nProcConsistencyIter; i++)
234  {
235  const lduMesh& mesh = meshLevel(nCreatedLevels);
236  const lduAddressing& addr = mesh.lduAddr();
237  const lduInterfacePtrsList interfaces = mesh.interfaces();
238 
239  const labelField& agglom = finalAgglomPtr();
240 
241  // Global numbering
242  const globalIndex globalNumbering(nCoarseCells);
243 
244  labelField globalAgglom(addr.size());
245  forAll(agglom, celli)
246  {
247  globalAgglom[celli] = globalNumbering.toGlobal(agglom[celli]);
248  }
249 
250  // Get the interface agglomeration
251  PtrList<labelList> nbrGlobalAgglom;
252  swap(interfaces, globalAgglom, nbrGlobalAgglom);
253 
254 
255  // Get the interface agglomeration on a cell basis (-1 for all
256  // other cells)
257  labelList cellToNbrAgglom;
258  getNbrAgglom(addr, interfaces, nbrGlobalAgglom, cellToNbrAgglom);
259 
260 
261  // Mark all faces in between cells with same nbragglomeration
262  labelHashSet sharedFaces(addr.size()/100);
263  detectSharedFaces(mesh, cellToNbrAgglom, sharedFaces);
264 
265 
266  //- Note: in-place update of weights is more effective it seems?
267  // Should not be. fluke?
268  // scalarField weights(*faceWeightsPtr);
269  scalarField weights = *magSfPtr;
270  forAllConstIter(labelHashSet, sharedFaces, iter)
271  {
272  label facei= iter.key();
273  weights[facei] *= 2.0;
274  }
275 
276  // Redo the agglomeration using the new weights
277  finalAgglomPtr = agglomerate
278  (
279  nCoarseCells,
280  minSize,
281  maxSize,
282  meshLevel(nCreatedLevels).lduAddr(),
283  *VPtr,
284  weights,
285  *magSbPtr
286  );
287  }
288 
289  if (continueAgglomerating(finalAgglomPtr().size(), nCoarseCells))
290  {
291  nCells_[nCreatedLevels] = nCoarseCells;
292  restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
293  }
294  else
295  {
296  break;
297  }
298 
299  agglomerateLduAddressing(nCreatedLevels);
300 
301  // Agglomerate the cell volumes field for the next level
302  {
303  scalarField* aggVPtr
304  (
305  new scalarField(meshLevels_[nCreatedLevels].size())
306  );
307 
308  // Restrict but no parallel agglomeration (not supported)
309  restrictField(*aggVPtr, *VPtr, nCreatedLevels, false);
310 
311  if (nCreatedLevels)
312  {
313  delete VPtr;
314  }
315 
316  VPtr = aggVPtr;
317  }
318 
319  // Agglomerate the face areas field for the next level
320  {
321  scalarField* aggMagSfPtr
322  (
323  new scalarField
324  (
325  meshLevels_[nCreatedLevels].upperAddr().size(),
326  0
327  )
328  );
329 
330  restrictFaceField(*aggMagSfPtr, *magSfPtr, nCreatedLevels);
331 
332  if (nCreatedLevels)
333  {
334  delete magSfPtr;
335  }
336 
337  magSfPtr = aggMagSfPtr;
338  }
339 
340  // Agglomerate the cell boundary areas field for the next level
341  {
342  scalarField* aggMagSbPtr
343  (
344  new scalarField(meshLevels_[nCreatedLevels].size())
345  );
346 
347  // Restrict but no parallel agglomeration (not supported)
348  restrictField(*aggMagSbPtr, *magSbPtr, nCreatedLevels, false);
349 
350  delete magSbPtr;
351  magSbPtr = aggMagSbPtr;
352  }
353 
354  nCreatedLevels++;
355  }
356 
357  // Shrink the storage of the levels to those created
358  compactLevels(nCreatedLevels);
359 
360  // Delete temporary geometry storage
361  if (nCreatedLevels)
362  {
363  delete VPtr;
364  delete magSfPtr;
365  }
366  delete magSbPtr;
367 }
368 
369 
370 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
Macros for easy insertion into run-time selection tables.
Geometric agglomerated algebraic multigrid agglomeration class.
void restrictField(Field< Type > &cf, const Field< Type > &ff, const label fineLevelIndex, const bool procAgglom) const
Restrict (integrate by summation) cell field.
void agglomerateLduAddressing(const label fineLevelIndex)
Assemble coarse mesh addressing.
const label maxLevels_
Max number of levels.
void compactLevels(const label nCreatedLevels)
Shrink the number of levels to that specified.
PtrList< lduPrimitiveMesh > meshLevels_
Hierarchy of mesh addressing.
bool continueAgglomerating(const label nCells, const label nCoarseCells) const
Check the need for further agglomeration.
labelList nCells_
The number of cells in each level.
void restrictFaceField(Field< Type > &cf, const Field< Type > &ff, const label fineLevelIndex) const
Restrict (integrate by summation) face field.
PtrList< labelField > restrictAddressing_
Cell restriction addressing array.
const lduMesh & meshLevel(const label leveli) const
Return LDU mesh of given level.
Agglomerate using the MGridGen algorithm.
MGridGenGAMGAgglomeration(const lduMesh &mesh, const dictionary &controlDict)
Construct given mesh and controls.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: PtrList.H:75
Pre-declare related SubField type.
Definition: SubField.H:63
static void waitRequests(const label start=0)
Wait until all requests (from start onwards) have finished.
Definition: UPstream.C:147
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:99
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:64
label toGlobal(const label i) const
From local to global.
Definition: globalIndexI.H:82
The class contains the addressing required by the lduMatrix: upper, lower and losort.
label size() const
Return number of equations.
Abstract base class for meshes which provide LDU addressing for the construction of lduMatrix and LDU...
Definition: lduMesh.H:60
virtual lduInterfacePtrsList interfaces() const =0
Return a list of pointers for each patch.
virtual const lduAddressing & lduAddr() const =0
Return ldu addressing.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1339
label nCells() const
const scalarField & cellVolumes() const
label nInternalFaces() const
const scalarField & magFaceAreas() const
const vectorField & faceAreas() const
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
A class for managing temporary objects.
Definition: tmp.H:55
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: labelList.H:56
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
addToRunTimeSelectionTable(polyPatch, mergedCyclicPolyPatch, word)
To & refCast(From &r)
Reference type cast template function.
Definition: typeInfo.H:129
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
UPtrList< const lduInterface > lduInterfacePtrsList
List of coupled interface fields to be used in coupling.
defineTypeNameAndDebug(combustionModel, 0)
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:213
UList< label > labelUList
Definition: UList.H:65
runTime controlDict().lookup("adjustTimeStep") >> adjustTimeStep