GAMGSolverSolve.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 "PCG.H"
28 #include "PBiCGStab.H"
29 #include "SubField.H"
30 
31 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
32 
34 (
35  scalarField& psi,
36  const scalarField& source,
37  const direction cmpt
38 ) const
39 {
40  // Setup class containing solver performance data
41  solverPerformance solverPerf(typeName, fieldName_);
42 
43  // Calculate A.psi used to calculate the initial residual
44  scalarField Apsi(psi.size());
45  matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
46 
47  // Create the storage for the finestCorrection which may be used as a
48  // temporary in normFactor
49  scalarField finestCorrection(psi.size());
50 
51  // Calculate normalisation factor
52  scalar normFactor = this->normFactor(psi, source, Apsi, finestCorrection);
53 
54  if (debug >= 2)
55  {
56  Pout<< " Normalisation factor = " << normFactor << endl;
57  }
58 
59  // Calculate initial finest-grid residual field
60  scalarField finestResidual(source - Apsi);
61 
62  // Calculate normalised residual for convergence test
63  solverPerf.initialResidual() = gSumMag
64  (
65  finestResidual,
66  matrix().mesh().comm()
67  )/normFactor;
68  solverPerf.finalResidual() = solverPerf.initialResidual();
69 
70 
71  // Check convergence, solve if not converged
72  if
73  (
74  minIter_ > 0
75  || !solverPerf.checkConvergence(tolerance_, relTol_)
76  )
77  {
78  // Create coarse grid correction fields
79  PtrList<scalarField> coarseCorrFields;
80 
81  // Create coarse grid sources
82  PtrList<scalarField> coarseSources;
83 
84  // Create the smoothers for all levels
86 
87  // Scratch fields if processor-agglomerated coarse level meshes
88  // are bigger than original. Usually not needed
89  scalarField scratch1;
90  scalarField scratch2;
91 
92  // Initialise the above data structures
93  initVcycle
94  (
95  coarseCorrFields,
96  coarseSources,
97  smoothers,
98  scratch1,
99  scratch2
100  );
101 
102  do
103  {
104  Vcycle
105  (
106  smoothers,
107  psi,
108  source,
109  Apsi,
110  finestCorrection,
111  finestResidual,
112 
113  (scratch1.size() ? scratch1 : Apsi),
114  (scratch2.size() ? scratch2 : finestCorrection),
115 
116  coarseCorrFields,
117  coarseSources,
118  cmpt
119  );
120 
121  // Calculate finest level residual field
122  matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
123  finestResidual = source;
124  finestResidual -= Apsi;
125 
126  solverPerf.finalResidual() = gSumMag
127  (
128  finestResidual,
129  matrix().mesh().comm()
130  )/normFactor;
131 
132  if (debug >= 2)
133  {
134  solverPerf.print(Info.masterStream(matrix().mesh().comm()));
135  }
136  } while
137  (
138  (
139  ++solverPerf.nIterations() < maxIter_
140  && !solverPerf.checkConvergence(tolerance_, relTol_)
141  )
142  || solverPerf.nIterations() < minIter_
143  );
144  }
145 
146  return solverPerf;
147 }
148 
149 
150 void Foam::GAMGSolver::Vcycle
151 (
152  const PtrList<lduMatrix::smoother>& smoothers,
153  scalarField& psi,
154  const scalarField& source,
155  scalarField& Apsi,
156  scalarField& finestCorrection,
157  scalarField& finestResidual,
158 
159  scalarField& scratch1,
160  scalarField& scratch2,
161 
162  PtrList<scalarField>& coarseCorrFields,
163  PtrList<scalarField>& coarseSources,
164  const direction cmpt
165 ) const
166 {
167  // debug = 2;
168 
169  const label coarsestLevel = matrixLevels_.size() - 1;
170 
171  // Restrict finest grid residual for the next level up.
172  agglomeration_.restrictField(coarseSources[0], finestResidual, 0, true);
173 
174  if (debug >= 2 && nPreSweeps_)
175  {
176  Pout<< "Pre-smoothing scaling factors: ";
177  }
178 
179 
180  // Residual restriction (going to coarser levels)
181  for (label leveli = 0; leveli < coarsestLevel; leveli++)
182  {
183  if (coarseSources.set(leveli + 1))
184  {
185  // If the optional pre-smoothing sweeps are selected
186  // smooth the coarse-grid field for the restricted source
187  if (nPreSweeps_)
188  {
189  coarseCorrFields[leveli] = 0.0;
190 
191  smoothers[leveli + 1].smooth
192  (
193  coarseCorrFields[leveli],
194  coarseSources[leveli],
195  cmpt,
196  min
197  (
198  nPreSweeps_ + preSweepsLevelMultiplier_*leveli,
199  maxPreSweeps_
200  )
201  );
202 
204  (
205  scratch1,
206  coarseCorrFields[leveli].size()
207  );
208 
209  // Scale coarse-grid correction field
210  // but not on the coarsest level because it evaluates to 1
211  if (scaleCorrection_ && leveli < coarsestLevel - 1)
212  {
213  scale
214  (
215  coarseCorrFields[leveli],
216  const_cast<scalarField&>
217  (
218  ACf.operator const scalarField&()
219  ),
220  matrixLevels_[leveli],
221  interfaceLevelsBouCoeffs_[leveli],
222  interfaceLevels_[leveli],
223  coarseSources[leveli],
224  cmpt
225  );
226  }
227 
228  // Correct the residual with the new solution
229  matrixLevels_[leveli].Amul
230  (
231  const_cast<scalarField&>
232  (
233  ACf.operator const scalarField&()
234  ),
235  coarseCorrFields[leveli],
236  interfaceLevelsBouCoeffs_[leveli],
237  interfaceLevels_[leveli],
238  cmpt
239  );
240 
241  coarseSources[leveli] -= ACf;
242  }
243 
244  // Residual is equal to source
245  agglomeration_.restrictField
246  (
247  coarseSources[leveli + 1],
248  coarseSources[leveli],
249  leveli + 1,
250  true
251  );
252  }
253  }
254 
255  if (debug >= 2 && nPreSweeps_)
256  {
257  Pout<< endl;
258  }
259 
260 
261  // Solve Coarsest level with either an iterative or direct solver
262  if (coarseCorrFields.set(coarsestLevel))
263  {
264  solveCoarsestLevel
265  (
266  coarseCorrFields[coarsestLevel],
267  coarseSources[coarsestLevel]
268  );
269  }
270 
271  if (debug >= 2)
272  {
273  Pout<< "Post-smoothing scaling factors: ";
274  }
275 
276  // Smoothing and prolongation of the coarse correction fields
277  // (going to finer levels)
278 
279  scalarField dummyField(0);
280 
281  for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
282  {
283  if (coarseCorrFields.set(leveli))
284  {
285  // Create a field for the pre-smoothed correction field
286  // as a sub-field of the finestCorrection which is not
287  // currently being used
288  scalarField::subField preSmoothedCoarseCorrField
289  (
290  scratch2,
291  coarseCorrFields[leveli].size()
292  );
293 
294  // Only store the preSmoothedCoarseCorrField if pre-smoothing is
295  // used
296  if (nPreSweeps_)
297  {
298  preSmoothedCoarseCorrField = coarseCorrFields[leveli];
299  }
300 
301  agglomeration_.prolongField
302  (
303  coarseCorrFields[leveli],
304  (
305  coarseCorrFields.set(leveli + 1)
306  ? coarseCorrFields[leveli + 1]
307  : dummyField // dummy value
308  ),
309  leveli + 1,
310  true
311  );
312 
313 
314  // Create A.psi for this coarse level as a sub-field of Apsi
316  (
317  scratch1,
318  coarseCorrFields[leveli].size()
319  );
320  scalarField& ACfRef =
321  const_cast<scalarField&>(ACf.operator const scalarField&());
322 
323  if (interpolateCorrection_) //&& leveli < coarsestLevel - 2)
324  {
325  if (coarseCorrFields.set(leveli+1))
326  {
328  (
329  coarseCorrFields[leveli],
330  ACfRef,
331  matrixLevels_[leveli],
332  interfaceLevelsBouCoeffs_[leveli],
333  interfaceLevels_[leveli],
334  agglomeration_.restrictAddressing(leveli + 1),
335  coarseCorrFields[leveli + 1],
336  cmpt
337  );
338  }
339  else
340  {
342  (
343  coarseCorrFields[leveli],
344  ACfRef,
345  matrixLevels_[leveli],
346  interfaceLevelsBouCoeffs_[leveli],
347  interfaceLevels_[leveli],
348  cmpt
349  );
350  }
351  }
352 
353  // Scale coarse-grid correction field
354  // but not on the coarsest level because it evaluates to 1
355  if
356  (
357  scaleCorrection_
358  && (interpolateCorrection_ || leveli < coarsestLevel - 1)
359  )
360  {
361  scale
362  (
363  coarseCorrFields[leveli],
364  ACfRef,
365  matrixLevels_[leveli],
366  interfaceLevelsBouCoeffs_[leveli],
367  interfaceLevels_[leveli],
368  coarseSources[leveli],
369  cmpt
370  );
371  }
372 
373  // Only add the preSmoothedCoarseCorrField if pre-smoothing is
374  // used
375  if (nPreSweeps_)
376  {
377  coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
378  }
379 
380  smoothers[leveli + 1].smooth
381  (
382  coarseCorrFields[leveli],
383  coarseSources[leveli],
384  cmpt,
385  min
386  (
387  nPostSweeps_ + postSweepsLevelMultiplier_*leveli,
388  maxPostSweeps_
389  )
390  );
391  }
392  }
393 
394  // Prolong the finest level correction
395  agglomeration_.prolongField
396  (
397  finestCorrection,
398  coarseCorrFields[0],
399  0,
400  true
401  );
402 
403  if (interpolateCorrection_)
404  {
406  (
407  finestCorrection,
408  Apsi,
409  matrix_,
411  interfaces_,
412  agglomeration_.restrictAddressing(0),
413  coarseCorrFields[0],
414  cmpt
415  );
416  }
417 
418  if (scaleCorrection_)
419  {
420  // Scale the finest level correction
421  scale
422  (
423  finestCorrection,
424  Apsi,
425  matrix_,
427  interfaces_,
428  finestResidual,
429  cmpt
430  );
431  }
432 
433  forAll(psi, i)
434  {
435  psi[i] += finestCorrection[i];
436  }
437 
438  smoothers[0].smooth
439  (
440  psi,
441  source,
442  cmpt,
443  nFinestSweeps_
444  );
445 }
446 
447 
448 void Foam::GAMGSolver::initVcycle
449 (
450  PtrList<scalarField>& coarseCorrFields,
451  PtrList<scalarField>& coarseSources,
452  PtrList<lduMatrix::smoother>& smoothers,
453  scalarField& scratch1,
454  scalarField& scratch2
455 ) const
456 {
457  label maxSize = matrix_.diag().size();
458 
459  coarseCorrFields.setSize(matrixLevels_.size());
460  coarseSources.setSize(matrixLevels_.size());
461  smoothers.setSize(matrixLevels_.size() + 1);
462 
463  // Create the smoother for the finest level
464  smoothers.set
465  (
466  0,
468  (
469  fieldName_,
470  matrix_,
473  interfaces_,
475  )
476  );
477 
478  forAll(matrixLevels_, leveli)
479  {
480  if (agglomeration_.nCells(leveli) >= 0)
481  {
482  label nCoarseCells = agglomeration_.nCells(leveli);
483 
484  coarseSources.set(leveli, new scalarField(nCoarseCells));
485  }
486 
487  if (matrixLevels_.set(leveli))
488  {
489  const lduMatrix& mat = matrixLevels_[leveli];
490 
491  label nCoarseCells = mat.diag().size();
492 
493  maxSize = max(maxSize, nCoarseCells);
494 
495  coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
496 
497  smoothers.set
498  (
499  leveli + 1,
501  (
502  fieldName_,
503  matrixLevels_[leveli],
504  interfaceLevelsBouCoeffs_[leveli],
505  interfaceLevelsIntCoeffs_[leveli],
506  interfaceLevels_[leveli],
508  )
509  );
510  }
511  }
512 
513  if (maxSize > matrix_.diag().size())
514  {
515  // Allocate some scratch storage
516  scratch1.setSize(maxSize);
517  scratch2.setSize(maxSize);
518  }
519 }
520 
521 
522 Foam::dictionary Foam::GAMGSolver::PCGsolverDict
523 (
524  const scalar tol,
525  const scalar relTol
526 ) const
527 {
528  dictionary dict(IStringStream("solver PCG; preconditioner DIC;")());
529  dict.add("tolerance", tol);
530  dict.add("relTol", relTol);
531 
532  return dict;
533 }
534 
535 
536 Foam::dictionary Foam::GAMGSolver::PBiCGStabSolverDict
537 (
538  const scalar tol,
539  const scalar relTol
540 ) const
541 {
542  dictionary dict(IStringStream("solver PBiCGStab; preconditioner DILU;")());
543  dict.add("tolerance", tol);
544  dict.add("relTol", relTol);
545 
546  return dict;
547 }
548 
549 
550 void Foam::GAMGSolver::solveCoarsestLevel
551 (
552  scalarField& coarsestCorrField,
553  const scalarField& coarsestSource
554 ) const
555 {
556  const label coarsestLevel = matrixLevels_.size() - 1;
557 
558  label coarseComm = matrixLevels_[coarsestLevel].mesh().comm();
559 
560  if (directSolveCoarsest_)
561  {
562  coarsestLUMatrixPtr_->solve(coarsestCorrField, coarsestSource);
563  }
564  // else if
565  //(
566  // agglomeration_.processorAgglomerate()
567  // && procMatrixLevels_.set(coarsestLevel)
568  //)
569  //{
570  // // const labelList& agglomProcIDs = agglomeration_.agglomProcIDs
571  // //(
572  // // coarsestLevel
573  // //);
574  // //
575  // // scalarField allSource;
576  // //
577  // // globalIndex cellOffsets;
578  // // if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
579  // //{
580  // // cellOffsets.offsets() =
581  // // agglomeration_.cellOffsets(coarsestLevel);
582  // //}
583  // //
584  // // cellOffsets.gather
585  // //(
586  // // coarseComm,
587  // // agglomProcIDs,
588  // // coarsestSource,
589  // // allSource
590  // //);
591  // //
592  // // scalarField allCorrField;
593  // // solverPerformance coarseSolverPerf;
594  //
595  // label solveComm = agglomeration_.procCommunicator(coarsestLevel);
596  //
597  // coarsestCorrField = 0;
598  // solverPerformance coarseSolverPerf;
599  //
600  // if (Pstream::myProcNo(solveComm) != -1)
601  // {
602  // const lduMatrix& allMatrix = procMatrixLevels_[coarsestLevel];
603  //
604  // {
605  // Pout<< "** Master:Solving on comm:" << solveComm
606  // << " with procs:" << UPstream::procID(solveComm) << endl;
607  //
608  // if (allMatrix.asymmetric())
609  // {
610  // coarseSolverPerf = PBiCGStab
611  // (
612  // "coarsestLevelCorr",
613  // allMatrix,
614  // procInterfaceLevelsBouCoeffs_[coarsestLevel],
615  // procInterfaceLevelsIntCoeffs_[coarsestLevel],
616  // procInterfaceLevels_[coarsestLevel],
617  // PBiCGStabSolverDict(tolerance_, relTol_)
618  // ).solve
619  // (
620  // coarsestCorrField,
621  // coarsestSource
622  // );
623  // }
624  // else
625  // {
626  // coarseSolverPerf = PCG
627  // (
628  // "coarsestLevelCorr",
629  // allMatrix,
630  // procInterfaceLevelsBouCoeffs_[coarsestLevel],
631  // procInterfaceLevelsIntCoeffs_[coarsestLevel],
632  // procInterfaceLevels_[coarsestLevel],
633  // PCGsolverDict(tolerance_, relTol_)
634  // ).solve
635  // (
636  // coarsestCorrField,
637  // coarsestSource
638  // );
639  // }
640  // }
641  // }
642  //
643  // Pout<< "done master solve." << endl;
644  //
645  // //// Scatter to all processors
646  // // coarsestCorrField.setSize(coarsestSource.size());
647  // // cellOffsets.scatter
648  // //(
649  // // coarseComm,
650  // // agglomProcIDs,
651  // // allCorrField,
652  // // coarsestCorrField
653  // //);
654  //
655  // if (debug >= 2)
656  // {
657  // coarseSolverPerf.print(Info.masterStream(coarseComm));
658  // }
659  //
660  // Pout<< "procAgglom: coarsestSource :" << coarsestSource << endl;
661  // Pout<< "procAgglom: coarsestCorrField:" << coarsestCorrField << endl;
662  //}
663  else
664  {
665  coarsestCorrField = 0;
666  solverPerformance coarseSolverPerf;
667 
668  if (matrixLevels_[coarsestLevel].asymmetric())
669  {
670  coarseSolverPerf = PBiCGStab
671  (
672  "coarsestLevelCorr",
673  matrixLevels_[coarsestLevel],
674  interfaceLevelsBouCoeffs_[coarsestLevel],
675  interfaceLevelsIntCoeffs_[coarsestLevel],
676  interfaceLevels_[coarsestLevel],
677  PBiCGStabSolverDict(tolerance_, relTol_)
678  ).solve
679  (
680  coarsestCorrField,
681  coarsestSource
682  );
683  }
684  else
685  {
686  coarseSolverPerf = PCG
687  (
688  "coarsestLevelCorr",
689  matrixLevels_[coarsestLevel],
690  interfaceLevelsBouCoeffs_[coarsestLevel],
691  interfaceLevelsIntCoeffs_[coarsestLevel],
692  interfaceLevels_[coarsestLevel],
693  PCGsolverDict(tolerance_, relTol_)
694  ).solve
695  (
696  coarsestCorrField,
697  coarsestSource
698  );
699  }
700 
701  if (debug >= 2)
702  {
703  coarseSolverPerf.print(Info.masterStream(coarseComm));
704  }
705  }
706 }
707 
708 
709 // ************************************************************************* //
scalar gSumMag(const FieldField< Field, Type > &f)
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
void restrictField(Field< Type > &cf, const Field< Type > &ff, const label fineLevelIndex, const bool procAgglom) const
Restrict (integrate by summation) cell field.
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
const FieldField< Field, scalar > & interfaceIntCoeffs_
Definition: lduMatrix.H:102
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
const labelType & nIterations() const
Return number of iterations.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const FieldField< Field, scalar > & interfaceBouCoeffs_
Definition: lduMatrix.H:101
uint8_t direction
Definition: direction.H:45
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve.
const Type & finalResidual() const
Return final residual.
Pre-declare related SubField type.
Definition: Field.H:60
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PCG.C:66
bool add(entry *, bool mergeEntry=false)
Add a new entry.
Definition: dictionary.C:1133
dynamicFvMesh & mesh
scalar tolerance_
Final convergence tolerance.
Definition: lduMatrix.H:118
bool checkConvergence(const Type &tolerance, const Type &relTolerance)
Check, store and return convergence.
const Type & initialResidual() const
Return initial residual.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics...
void prolongField(Field< Type > &ff, const Field< Type > &cf, const label coarseLevelIndex, const bool procAgglom) const
Prolong (interpolate by injection) cell field.
const lduMatrix & matrix_
Definition: lduMatrix.H:100
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:131
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PBiCGStab.C:69
lduInterfaceFieldPtrsList interfaces_
Definition: lduMatrix.H:103
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
dictionary controlDict_
Dictionary of controls.
Definition: lduMatrix.H:106
void print(Ostream &os) const
Print summary of solver performance to the given stream.
OSstream & masterStream(const label communicator)
Convert to OSstream.
Definition: messageStream.C:59
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
fvModels source(alpha1, mixture.thermo1().rho())
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
Input from memory buffer stream.
Definition: IStringStream.H:49
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
static autoPtr< smoother > New(const word &fieldName, const lduMatrix &matrix, const FieldField< Field, scalar > &interfaceBouCoeffs, const FieldField< Field, scalar > &interfaceIntCoeffs, const lduInterfaceFieldPtrsList &interfaces, const dictionary &solverControls)
Return a new smoother.
const labelField & restrictAddressing(const label leveli) const
Return cell restrict addressing of given level.
messageStream Info
scalarField & diag()
Definition: lduMatrix.C:186
Preconditioned conjugate gradient solver for symmetric lduMatrices using a run-time selectable precon...
Definition: PCG.H:49
Preconditioned bi-conjugate gradient stabilised solver for asymmetric lduMatrices using a run-time se...
Definition: PBiCGStab.H:64
label nCells(const label leveli) const
Return number of coarse cells (before processor agglomeration)
scalar relTol_
Convergence tolerance relative to the initial.
Definition: lduMatrix.H:121