GAMGSolverSolve.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-2016 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 "PBiCG.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 restriced 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::PBiCGsolverDict
537 (
538  const scalar tol,
539  const scalar relTol
540 ) const
541 {
542  dictionary dict(IStringStream("solver PBiCG; 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  label oldWarn = UPstream::warnComm;
560  UPstream::warnComm = coarseComm;
561 
562  if (directSolveCoarsest_)
563  {
564  coarsestLUMatrixPtr_->solve(coarsestCorrField, coarsestSource);
565  }
566  //else if
567  //(
568  // agglomeration_.processorAgglomerate()
569  // && procMatrixLevels_.set(coarsestLevel)
570  //)
571  //{
572  // //const labelList& agglomProcIDs = agglomeration_.agglomProcIDs
573  // //(
574  // // coarsestLevel
575  // //);
576  // //
577  // //scalarField allSource;
578  // //
579  // //globalIndex cellOffsets;
580  // //if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
581  // //{
582  // // cellOffsets.offsets() =
583  // // agglomeration_.cellOffsets(coarsestLevel);
584  // //}
585  // //
586  // //cellOffsets.gather
587  // //(
588  // // coarseComm,
589  // // agglomProcIDs,
590  // // coarsestSource,
591  // // allSource
592  // //);
593  // //
594  // //scalarField allCorrField;
595  // //solverPerformance coarseSolverPerf;
596  //
597  // label solveComm = agglomeration_.procCommunicator(coarsestLevel);
598  // label oldWarn = UPstream::warnComm;
599  // UPstream::warnComm = solveComm;
600  //
601  //
602  // coarsestCorrField = 0;
603  // solverPerformance coarseSolverPerf;
604  //
605  // if (Pstream::myProcNo(solveComm) != -1)
606  // {
607  // const lduMatrix& allMatrix = procMatrixLevels_[coarsestLevel];
608  //
609  // {
610  // Pout<< "** Master:Solving on comm:" << solveComm
611  // << " with procs:" << UPstream::procID(solveComm) << endl;
612  //
613  // if (allMatrix.asymmetric())
614  // {
615  // coarseSolverPerf = PBiCG
616  // (
617  // "coarsestLevelCorr",
618  // allMatrix,
619  // procInterfaceLevelsBouCoeffs_[coarsestLevel],
620  // procInterfaceLevelsIntCoeffs_[coarsestLevel],
621  // procInterfaceLevels_[coarsestLevel],
622  // PBiCGsolverDict(tolerance_, relTol_)
623  // ).solve
624  // (
625  // coarsestCorrField,
626  // coarsestSource
627  // );
628  // }
629  // else
630  // {
631  // coarseSolverPerf = PCG
632  // (
633  // "coarsestLevelCorr",
634  // allMatrix,
635  // procInterfaceLevelsBouCoeffs_[coarsestLevel],
636  // procInterfaceLevelsIntCoeffs_[coarsestLevel],
637  // procInterfaceLevels_[coarsestLevel],
638  // PCGsolverDict(tolerance_, relTol_)
639  // ).solve
640  // (
641  // coarsestCorrField,
642  // coarsestSource
643  // );
644  // }
645  // }
646  // }
647  //
648  // UPstream::warnComm = oldWarn;
649  // Pout<< "done master solve." << endl;
650  //
651  // //// Scatter to all processors
652  // //coarsestCorrField.setSize(coarsestSource.size());
653  // //cellOffsets.scatter
654  // //(
655  // // coarseComm,
656  // // agglomProcIDs,
657  // // allCorrField,
658  // // coarsestCorrField
659  // //);
660  //
661  // if (debug >= 2)
662  // {
663  // coarseSolverPerf.print(Info.masterStream(coarseComm));
664  // }
665  //
666  // Pout<< "procAgglom: coarsestSource :" << coarsestSource << endl;
667  // Pout<< "procAgglom: coarsestCorrField:" << coarsestCorrField << endl;
668  //}
669  else
670  {
671  coarsestCorrField = 0;
672  solverPerformance coarseSolverPerf;
673 
674  if (matrixLevels_[coarsestLevel].asymmetric())
675  {
676  coarseSolverPerf = PBiCG
677  (
678  "coarsestLevelCorr",
679  matrixLevels_[coarsestLevel],
680  interfaceLevelsBouCoeffs_[coarsestLevel],
681  interfaceLevelsIntCoeffs_[coarsestLevel],
682  interfaceLevels_[coarsestLevel],
683  PBiCGsolverDict(tolerance_, relTol_)
684  ).solve
685  (
686  coarsestCorrField,
687  coarsestSource
688  );
689  }
690  else
691  {
692  coarseSolverPerf = PCG
693  (
694  "coarsestLevelCorr",
695  matrixLevels_[coarsestLevel],
696  interfaceLevelsBouCoeffs_[coarsestLevel],
697  interfaceLevelsIntCoeffs_[coarsestLevel],
698  interfaceLevels_[coarsestLevel],
699  PCGsolverDict(tolerance_, relTol_)
700  ).solve
701  (
702  coarsestCorrField,
703  coarsestSource
704  );
705  }
706 
707  if (debug >= 2)
708  {
709  coarseSolverPerf.print(Info.masterStream(coarseComm));
710  }
711  }
712 
713  UPstream::warnComm = oldWarn;
714 }
715 
716 
717 // ************************************************************************* //
scalar gSumMag(const FieldField< Field, Type > &f)
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
uint8_t direction
Definition: direction.H:46
const FieldField< Field, scalar > & interfaceIntCoeffs_
Definition: lduMatrix.H:102
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const FieldField< Field, scalar > & interfaceBouCoeffs_
Definition: lduMatrix.H:101
const Type & finalResidual() const
Return final residual.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
tmp< surfaceScalarField > interpolate(const RhoType &rho)
Pre-declare related SubField type.
Definition: Field.H:61
bool add(entry *, bool mergeEntry=false)
Add a new entry.
Definition: dictionary.C:737
dynamicFvMesh & mesh
scalar tolerance_
Final convergence tolerance.
Definition: lduMatrix.H:115
label nCells(const label leveli) const
Return number of coarse cells (before processor agglomeration)
void print(Ostream &os) const
Print summary of solver performance to the given stream.
bool checkConvergence(const Type &tolerance, const Type &relTolerance)
Check, store and return convergence.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics...
label nIterations() const
Return number of iterations.
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PCG.C:66
const labelField & restrictAddressing(const label leveli) const
Return cell restrict addressing of given level.
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PBiCG.C:66
static label warnComm
Debugging: warn for use of any communicator differing from warnComm.
Definition: UPstream.H:277
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
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices using a run-time selectable pr...
Definition: PBiCG.H:49
lduInterfaceFieldPtrsList interfaces_
Definition: lduMatrix.H:103
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dictionary controlDict_
Dictionary of controls.
Definition: lduMatrix.H:106
OSstream & masterStream(const label communicator)
Convert to OSstream.
Definition: messageStream.C:59
void setSize(const label)
Reset size of List.
Definition: List.C:295
lduMatrix is a general matrix class in which the coefficients are stored as three arrays...
Definition: lduMatrix.H:79
void prolongField(Field< Type > &ff, const Field< Type > &cf, const label coarseLevelIndex, const bool procAgglom) const
Prolong (interpolate by injection) cell field.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:62
Input from memory buffer stream.
Definition: IStringStream.H:49
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 Type & initialResidual() const
Return initial residual.
messageStream Info
void restrictField(Field< Type > &cf, const Field< Type > &ff, const label fineLevelIndex, const bool procAgglom) const
Restrict (integrate by summation) cell field.
scalarField & diag()
Definition: lduMatrix.C:183
Preconditioned conjugate gradient solver for symmetric lduMatrices using a run-time selectable precon...
Definition: PCG.H:49
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve.
scalar relTol_
Convergence tolerance relative to the initial.
Definition: lduMatrix.H:118