GAMGSolver.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) 2011-2015 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 {
33  defineTypeNameAndDebug(GAMGSolver, 0);
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);
220  const lduInterfaceFieldPtrsList& interfaces =
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  const lduMesh& coarsestMesh =
257  matrixLevels_[coarsestLevel].mesh();
258 
259  label coarseComm = coarsestMesh.comm();
260  label oldWarn = UPstream::warnComm;
261  UPstream::warnComm = coarseComm;
262 
263  coarsestLUMatrixPtr_.set
264  (
265  new LUscalarMatrix
266  (
267  matrixLevels_[coarsestLevel],
268  interfaceLevelsBouCoeffs_[coarsestLevel],
269  interfaceLevels_[coarsestLevel]
270  )
271  );
272 
273  UPstream::warnComm = oldWarn;
274  }
275  }
276  }
277  else
278  {
280  (
281  "GAMGSolver::GAMGSolver"
282  "("
283  "const word& fieldName,"
284  "const lduMatrix& matrix,"
285  "const FieldField<Field, scalar>& interfaceBouCoeffs,"
286  "const FieldField<Field, scalar>& interfaceIntCoeffs,"
287  "const lduInterfaceFieldPtrsList& interfaces,"
288  "const dictionary& solverControls"
289  ")"
290  ) << "No coarse levels created, either matrix too small for GAMG"
291  " or nCellsInCoarsestLevel too large.\n"
292  " Either choose another solver of reduce "
293  "nCellsInCoarsestLevel."
294  << exit(FatalError);
295  }
296 }
297 
298 
299 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
300 
302 {
303  if (!cacheAgglomeration_)
304  {
305  delete &agglomeration_;
306  }
307 }
308 
309 
310 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
311 
312 void Foam::GAMGSolver::readControls()
313 {
315 
316  controlDict_.readIfPresent("cacheAgglomeration", cacheAgglomeration_);
317  controlDict_.readIfPresent("nPreSweeps", nPreSweeps_);
318  controlDict_.readIfPresent
319  (
320  "preSweepsLevelMultiplier",
321  preSweepsLevelMultiplier_
322  );
323  controlDict_.readIfPresent("maxPreSweeps", maxPreSweeps_);
324  controlDict_.readIfPresent("nPostSweeps", nPostSweeps_);
325  controlDict_.readIfPresent
326  (
327  "postSweepsLevelMultiplier",
328  postSweepsLevelMultiplier_
329  );
330  controlDict_.readIfPresent("maxPostSweeps", maxPostSweeps_);
331  controlDict_.readIfPresent("nFinestSweeps", nFinestSweeps_);
332  controlDict_.readIfPresent("interpolateCorrection", interpolateCorrection_);
333  controlDict_.readIfPresent("scaleCorrection", scaleCorrection_);
334  controlDict_.readIfPresent("directSolveCoarsest", directSolveCoarsest_);
335 
336  if (debug)
337  {
338  Pout<< "GAMGSolver settings :"
339  << " cacheAgglomeration:" << cacheAgglomeration_
340  << " nPreSweeps:" << nPreSweeps_
341  << " preSweepsLevelMultiplier:" << preSweepsLevelMultiplier_
342  << " maxPreSweeps:" << maxPreSweeps_
343  << " nPostSweeps:" << nPostSweeps_
344  << " postSweepsLevelMultiplier:" << postSweepsLevelMultiplier_
345  << " maxPostSweeps:" << maxPostSweeps_
346  << " nFinestSweeps:" << nFinestSweeps_
347  << " interpolateCorrection:" << interpolateCorrection_
348  << " scaleCorrection:" << scaleCorrection_
349  << " directSolveCoarsest:" << directSolveCoarsest_
350  << endl;
351  }
352 }
353 
354 
355 const Foam::lduMatrix& Foam::GAMGSolver::matrixLevel(const label i) const
356 {
357  if (i == 0)
358  {
359  return matrix_;
360  }
361  else
362  {
363  return matrixLevels_[i - 1];
364  }
365 }
366 
367 
368 const Foam::lduInterfaceFieldPtrsList& Foam::GAMGSolver::interfaceLevel
369 (
370  const label i
371 ) const
372 {
373  if (i == 0)
374  {
375  return interfaces_;
376  }
377  else
378  {
379  return interfaceLevels_[i - 1];
380  }
381 }
382 
383 
385 Foam::GAMGSolver::interfaceBouCoeffsLevel
386 (
387  const label i
388 ) const
389 {
390  if (i == 0)
391  {
392  return interfaceBouCoeffs_;
393  }
394  else
395  {
396  return interfaceLevelsBouCoeffs_[i - 1];
397  }
398 }
399 
400 
402 Foam::GAMGSolver::interfaceIntCoeffsLevel
403 (
404  const label i
405 ) const
406 {
407  if (i == 0)
408  {
409  return interfaceIntCoeffs_;
410  }
411  else
412  {
413  return interfaceLevelsIntCoeffs_[i - 1];
414  }
415 }
416 
417 
418 // ************************************************************************* //
Input from memory buffer stream.
Definition: IStringStream.H:49
virtual void readControls()
Read the control parameters from the controlDict_.
Generic field type.
Definition: FieldField.H:51
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:61
virtual label comm() const =0
Return communicator used for parallel communication.
lduMatrix::solver::addsymMatrixConstructorToTable< GAMGSolver > addGAMGSolverMatrixConstructorToTable_
Definition: GAMGSolver.C:36
scalarField & lower()
Definition: lduMatrix.C:165
static const GAMGAgglomeration & New(const lduMesh &mesh, const dictionary &controlDict)
Return the selected geometric agglomerator.
string str() const
Return the string.
A class for handling words, derived from string.
Definition: word.H:59
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Abstract base class for meshes which provide LDU addressing for the construction of lduMatrix and LDU...
Definition: lduMesh.H:51
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
Namespace for OpenFOAM.
runTime write()
Foam::LUscalarMatrix.
static const char nl
Definition: Ostream.H:260
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.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static label warnComm
Debugging: warn for use of any communicator differing from warnComm.
Definition: UPstream.H:264
virtual ~GAMGSolver()
Destructor.
Definition: GAMGSolver.C:301
lduMatrix is a general matrix class in which the coefficients are stored as three arrays...
Definition: lduMatrix.H:77
#define forAll(list, i)
Definition: UList.H:421
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
Output to memory buffer stream.
Definition: OStringStream.H:49
scalarField & diag()
Definition: lduMatrix.C:183
fileName::Type type(const fileName &)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:589
bool set(const label) const
Is element set.
Definition: UPtrListI.H:80
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
lduMatrix::solver::addasymMatrixConstructorToTable< GAMGSolver > addGAMGAsymSolverMatrixConstructorToTable_
Definition: GAMGSolver.C:39
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:31
error FatalError
bool symmetric() const
Definition: lduMatrix.H:600
Abstract base-class for lduMatrix solvers.
Definition: lduMatrix.H:91
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53