GAMGSolver.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-2018 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 "GAMGSolver.H"
27 #include "GAMGInterface.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
34 
35  lduMatrix::solver::addsymMatrixConstructorToTable<GAMGSolver>
37 
38  lduMatrix::solver::addasymMatrixConstructorToTable<GAMGSolver>
40 }
41 
42 
43 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
44 
46 (
47  const word& fieldName,
48  const lduMatrix& matrix,
49  const FieldField<Field, scalar>& interfaceBouCoeffs,
50  const FieldField<Field, scalar>& interfaceIntCoeffs,
51  const lduInterfaceFieldPtrsList& interfaces,
52  const dictionary& solverControls
53 )
54 :
56  (
57  fieldName,
58  matrix,
59  interfaceBouCoeffs,
60  interfaceIntCoeffs,
61  interfaces,
62  solverControls
63  ),
64 
65  // Default values for all controls
66  // which may be overridden by those in controlDict
67  cacheAgglomeration_(true),
68  nPreSweeps_(0),
69  preSweepsLevelMultiplier_(1),
70  maxPreSweeps_(4),
71  nPostSweeps_(2),
72  postSweepsLevelMultiplier_(1),
73  maxPostSweeps_(4),
74  nFinestSweeps_(2),
75  interpolateCorrection_(false),
76  scaleCorrection_(matrix.symmetric()),
77  directSolveCoarsest_(false),
78  agglomeration_(GAMGAgglomeration::New(matrix_, controlDict_)),
79 
80  matrixLevels_(agglomeration_.size()),
81  primitiveInterfaceLevels_(agglomeration_.size()),
82  interfaceLevels_(agglomeration_.size()),
83  interfaceLevelsBouCoeffs_(agglomeration_.size()),
84  interfaceLevelsIntCoeffs_(agglomeration_.size())
85 {
86  readControls();
87 
88  if (agglomeration_.processorAgglomerate())
89  {
90  forAll(agglomeration_, fineLevelIndex)
91  {
92  if (agglomeration_.hasMeshLevel(fineLevelIndex))
93  {
94  if
95  (
96  (fineLevelIndex+1) < agglomeration_.size()
97  && agglomeration_.hasProcMesh(fineLevelIndex+1)
98  )
99  {
100  // Construct matrix without referencing the coarse mesh so
101  // construct a dummy mesh instead. This will get overwritten
102  // by the call to procAgglomerateMatrix so is only to get
103  // it through agglomerateMatrix
104 
105 
106  const lduInterfacePtrsList& fineMeshInterfaces =
107  agglomeration_.interfaceLevel(fineLevelIndex);
108 
109  PtrList<GAMGInterface> dummyPrimMeshInterfaces
110  (
111  fineMeshInterfaces.size()
112  );
113  lduInterfacePtrsList dummyMeshInterfaces
114  (
115  dummyPrimMeshInterfaces.size()
116  );
117  forAll(fineMeshInterfaces, intI)
118  {
119  if (fineMeshInterfaces.set(intI))
120  {
121  OStringStream os;
122  refCast<const GAMGInterface>
123  (
124  fineMeshInterfaces[intI]
125  ).write(os);
126  IStringStream is(os.str());
127 
128  dummyPrimMeshInterfaces.set
129  (
130  intI,
132  (
133  fineMeshInterfaces[intI].type(),
134  intI,
135  dummyMeshInterfaces,
136  is
137  )
138  );
139  }
140  }
141 
142  forAll(dummyPrimMeshInterfaces, intI)
143  {
144  if (dummyPrimMeshInterfaces.set(intI))
145  {
146  dummyMeshInterfaces.set
147  (
148  intI,
149  &dummyPrimMeshInterfaces[intI]
150  );
151  }
152  }
153 
154  // So:
155  // - pass in incorrect mesh (= fine mesh instead of coarse)
156  // - pass in dummy interfaces
157  agglomerateMatrix
158  (
159  fineLevelIndex,
160  agglomeration_.meshLevel(fineLevelIndex),
161  dummyMeshInterfaces
162  );
163 
164 
165  const labelList& procAgglomMap =
166  agglomeration_.procAgglomMap(fineLevelIndex+1);
167  const List<label>& procIDs =
168  agglomeration_.agglomProcIDs(fineLevelIndex+1);
169 
170  procAgglomerateMatrix
171  (
172  procAgglomMap,
173  procIDs,
174  fineLevelIndex
175  );
176  }
177  else
178  {
179  agglomerateMatrix
180  (
181  fineLevelIndex,
182  agglomeration_.meshLevel(fineLevelIndex + 1),
183  agglomeration_.interfaceLevel(fineLevelIndex + 1)
184  );
185  }
186  }
187  else
188  {
189  // No mesh. Not involved in calculation anymore
190  }
191  }
192  }
193  else
194  {
195  forAll(agglomeration_, fineLevelIndex)
196  {
197  // Agglomerate on to coarse level mesh
198  agglomerateMatrix
199  (
200  fineLevelIndex,
201  agglomeration_.meshLevel(fineLevelIndex + 1),
202  agglomeration_.interfaceLevel(fineLevelIndex + 1)
203  );
204  }
205  }
206 
207 
208  if (debug)
209  {
210  for
211  (
212  label fineLevelIndex = 0;
213  fineLevelIndex <= matrixLevels_.size();
214  fineLevelIndex++
215  )
216  {
217  if (fineLevelIndex == 0 || matrixLevels_.set(fineLevelIndex-1))
218  {
219  const lduMatrix& matrix = matrixLevel(fineLevelIndex);
221  interfaceLevel(fineLevelIndex);
222 
223  Pout<< "level:" << fineLevelIndex << nl
224  << " nCells:" << matrix.diag().size() << nl
225  << " nFaces:" << matrix.lower().size() << nl
226  << " nInterfaces:" << interfaces.size()
227  << endl;
228 
229  forAll(interfaces, i)
230  {
231  if (interfaces.set(i))
232  {
233  Pout<< " " << i
234  << "\ttype:" << interfaces[i].type()
235  << endl;
236  }
237  }
238  }
239  else
240  {
241  Pout<< "level:" << fineLevelIndex << " : no matrix" << endl;
242  }
243  }
244  Pout<< endl;
245  }
246 
247 
248  if (matrixLevels_.size())
249  {
250  if (directSolveCoarsest_)
251  {
252  const label coarsestLevel = matrixLevels_.size() - 1;
253 
254  if (matrixLevels_.set(coarsestLevel))
255  {
256  coarsestLUMatrixPtr_.set
257  (
258  new LUscalarMatrix
259  (
260  matrixLevels_[coarsestLevel],
261  interfaceLevelsBouCoeffs_[coarsestLevel],
262  interfaceLevels_[coarsestLevel]
263  )
264  );
265  }
266  }
267  }
268  else
269  {
271  << "No coarse levels created, either matrix too small for GAMG"
272  " or nCellsInCoarsestLevel too large.\n"
273  " Either choose another solver of reduce "
274  "nCellsInCoarsestLevel."
275  << exit(FatalError);
276  }
277 }
278 
279 
280 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
281 
283 {
284  if (!cacheAgglomeration_)
285  {
286  delete &agglomeration_;
287  }
288 }
289 
290 
291 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
292 
293 void Foam::GAMGSolver::readControls()
294 {
296 
297  controlDict_.readIfPresent("cacheAgglomeration", cacheAgglomeration_);
298  controlDict_.readIfPresent("nPreSweeps", nPreSweeps_);
299  controlDict_.readIfPresent
300  (
301  "preSweepsLevelMultiplier",
302  preSweepsLevelMultiplier_
303  );
304  controlDict_.readIfPresent("maxPreSweeps", maxPreSweeps_);
305  controlDict_.readIfPresent("nPostSweeps", nPostSweeps_);
306  controlDict_.readIfPresent
307  (
308  "postSweepsLevelMultiplier",
309  postSweepsLevelMultiplier_
310  );
311  controlDict_.readIfPresent("maxPostSweeps", maxPostSweeps_);
312  controlDict_.readIfPresent("nFinestSweeps", nFinestSweeps_);
313  controlDict_.readIfPresent("interpolateCorrection", interpolateCorrection_);
314  controlDict_.readIfPresent("scaleCorrection", scaleCorrection_);
315  controlDict_.readIfPresent("directSolveCoarsest", directSolveCoarsest_);
316 
317  if (debug)
318  {
319  Pout<< "GAMGSolver settings :"
320  << " cacheAgglomeration:" << cacheAgglomeration_
321  << " nPreSweeps:" << nPreSweeps_
322  << " preSweepsLevelMultiplier:" << preSweepsLevelMultiplier_
323  << " maxPreSweeps:" << maxPreSweeps_
324  << " nPostSweeps:" << nPostSweeps_
325  << " postSweepsLevelMultiplier:" << postSweepsLevelMultiplier_
326  << " maxPostSweeps:" << maxPostSweeps_
327  << " nFinestSweeps:" << nFinestSweeps_
328  << " interpolateCorrection:" << interpolateCorrection_
329  << " scaleCorrection:" << scaleCorrection_
330  << " directSolveCoarsest:" << directSolveCoarsest_
331  << endl;
332  }
333 }
334 
335 
336 const Foam::lduMatrix& Foam::GAMGSolver::matrixLevel(const label i) const
337 {
338  if (i == 0)
339  {
340  return matrix_;
341  }
342  else
343  {
344  return matrixLevels_[i - 1];
345  }
346 }
347 
348 
349 const Foam::lduInterfaceFieldPtrsList& Foam::GAMGSolver::interfaceLevel
350 (
351  const label i
352 ) const
353 {
354  if (i == 0)
355  {
356  return interfaces_;
357  }
358  else
359  {
360  return interfaceLevels_[i - 1];
361  }
362 }
363 
364 
366 Foam::GAMGSolver::interfaceBouCoeffsLevel
367 (
368  const label i
369 ) const
370 {
371  if (i == 0)
372  {
373  return interfaceBouCoeffs_;
374  }
375  else
376  {
377  return interfaceLevelsBouCoeffs_[i - 1];
378  }
379 }
380 
381 
383 Foam::GAMGSolver::interfaceIntCoeffsLevel
384 (
385  const label i
386 ) const
387 {
388  if (i == 0)
389  {
390  return interfaceIntCoeffs_;
391  }
392  else
393  {
394  return interfaceLevelsIntCoeffs_[i - 1];
395  }
396 }
397 
398 
399 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Generic field type.
Definition: FieldField.H:77
Geometric agglomerated algebraic multigrid agglomeration class.
bool processorAgglomerate() const
Whether to agglomerate across processors.
const labelList & procAgglomMap(const label fineLeveli) const
Mapping from processor to agglomerated processor (global, all.
bool hasProcMesh(const label fineLeveli) const
Check that level has combined mesh.
const labelList & agglomProcIDs(const label fineLeveli) const
Set of processors to agglomerate. Element 0 is the.
const lduInterfacePtrsList & interfaceLevel(const label leveli) const
Return LDU interface addressing of given level.
bool hasMeshLevel(const label leveli) const
Do we have mesh for given level?
const lduMesh & meshLevel(const label leveli) const
Return LDU mesh of given level.
static autoPtr< GAMGInterface > New(const label index, const lduInterfacePtrsList &coarseInterfaces, const lduInterface &fineInterface, const labelField &localRestrictAddressing, const labelField &neighbourRestrictAddressing, const label fineLevelIndex, const label coarseComm)
Return a pointer to a new interface created on freestore given.
Geometric agglomerated algebraic multigrid solver.
Definition: GAMGSolver.H:73
GAMGSolver(const word &fieldName, const lduMatrix &matrix, const FieldField< Field, scalar > &interfaceBouCoeffs, const FieldField< Field, scalar > &interfaceIntCoeffs, const lduInterfaceFieldPtrsList &interfaces, const dictionary &solverControls)
Construct from lduMatrix and solver controls.
Definition: GAMGSolver.C:46
virtual ~GAMGSolver()
Destructor.
Definition: GAMGSolver.C:282
Input from memory buffer stream.
Definition: IStringStream.H:52
Class to perform the LU decomposition on a symmetric matrix.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
virtual Ostream & write(const char)
Write character.
Definition: OSstream.C:32
Output to memory buffer stream.
Definition: OStringStream.H:52
string str() const
Return the string.
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
bool set(const label) const
Is element set.
Definition: PtrListI.H:62
bool set(const label) const
Is element set.
Definition: UPtrListI.H:87
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
virtual const word & type() const =0
Runtime type information.
const lduInterfaceFieldPtrsList & interfaces() const
Definition: lduMatrix.H:241
virtual void readControls()
Read the control parameters from the controlDict_.
const lduMatrix & matrix() const
Definition: lduMatrix.H:226
lduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: lduMatrix.H:80
scalarField & lower()
Definition: lduMatrix.C:168
scalarField & diag()
Definition: lduMatrix.C:186
Abstract base class for run-time selectable region solvers.
Definition: solver.H:56
A class for handling words, derived from string.
Definition: word.H:62
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
defineTypeNameAndDebug(combustionModel, 0)
lduMatrix::solver::addasymMatrixConstructorToTable< GAMGSolver > addGAMGAsymSolverMatrixConstructorToTable_
Definition: GAMGSolver.C:39
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
error FatalError
static const char nl
Definition: Ostream.H:266
lduMatrix::solver::addsymMatrixConstructorToTable< GAMGSolver > addGAMGSolverMatrixConstructorToTable_
Definition: GAMGSolver.C:36