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-2014 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 "ICCG.H"
28 #include "BICCG.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.assign(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 void Foam::GAMGSolver::solveCoarsestLevel
523 (
524  scalarField& coarsestCorrField,
525  const scalarField& coarsestSource
526 ) const
527 {
528  const label coarsestLevel = matrixLevels_.size() - 1;
529 
530  label coarseComm = matrixLevels_[coarsestLevel].mesh().comm();
531  label oldWarn = UPstream::warnComm;
532  UPstream::warnComm = coarseComm;
533 
534  if (directSolveCoarsest_)
535  {
536  coarsestCorrField = coarsestSource;
537  coarsestLUMatrixPtr_->solve(coarsestCorrField);
538  }
539  //else if
540  //(
541  // agglomeration_.processorAgglomerate()
542  // && procMatrixLevels_.set(coarsestLevel)
543  //)
544  //{
545  // //const labelList& agglomProcIDs = agglomeration_.agglomProcIDs
546  // //(
547  // // coarsestLevel
548  // //);
549  // //
550  // //scalarField allSource;
551  // //
552  // //globalIndex cellOffsets;
553  // //if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
554  // //{
555  // // cellOffsets.offsets() =
556  // // agglomeration_.cellOffsets(coarsestLevel);
557  // //}
558  // //
559  // //cellOffsets.gather
560  // //(
561  // // coarseComm,
562  // // agglomProcIDs,
563  // // coarsestSource,
564  // // allSource
565  // //);
566  // //
567  // //scalarField allCorrField;
568  // //solverPerformance coarseSolverPerf;
569  //
570  // label solveComm = agglomeration_.procCommunicator(coarsestLevel);
571  // label oldWarn = UPstream::warnComm;
572  // UPstream::warnComm = solveComm;
573  //
574  //
575  // coarsestCorrField = 0;
576  // solverPerformance coarseSolverPerf;
577  //
578  // if (Pstream::myProcNo(solveComm) != -1)
579  // {
580  // const lduMatrix& allMatrix = procMatrixLevels_[coarsestLevel];
581  //
582  // {
583  // Pout<< "** Master:Solving on comm:" << solveComm
584  // << " with procs:" << UPstream::procID(solveComm) << endl;
585  //
586  // if (allMatrix.asymmetric())
587  // {
588  // coarseSolverPerf = BICCG
589  // (
590  // "coarsestLevelCorr",
591  // allMatrix,
592  // procInterfaceLevelsBouCoeffs_[coarsestLevel],
593  // procInterfaceLevelsIntCoeffs_[coarsestLevel],
594  // procInterfaceLevels_[coarsestLevel],
595  // tolerance_,
596  // relTol_
597  // ).solve
598  // (
599  // coarsestCorrField,
600  // coarsestSource
601  // );
602  // }
603  // else
604  // {
605  // coarseSolverPerf = ICCG
606  // (
607  // "coarsestLevelCorr",
608  // allMatrix,
609  // procInterfaceLevelsBouCoeffs_[coarsestLevel],
610  // procInterfaceLevelsIntCoeffs_[coarsestLevel],
611  // procInterfaceLevels_[coarsestLevel],
612  // tolerance_,
613  // relTol_
614  // ).solve
615  // (
616  // coarsestCorrField,
617  // coarsestSource
618  // );
619  // }
620  // }
621  // }
622  //
623  // UPstream::warnComm = oldWarn;
624  // Pout<< "done master solve." << endl;
625  //
626  // //// Scatter to all processors
627  // //coarsestCorrField.setSize(coarsestSource.size());
628  // //cellOffsets.scatter
629  // //(
630  // // coarseComm,
631  // // agglomProcIDs,
632  // // allCorrField,
633  // // coarsestCorrField
634  // //);
635  //
636  // if (debug >= 2)
637  // {
638  // coarseSolverPerf.print(Info.masterStream(coarseComm));
639  // }
640  //
641  // Pout<< "procAgglom: coarsestSource :" << coarsestSource << endl;
642  // Pout<< "procAgglom: coarsestCorrField:" << coarsestCorrField << endl;
643  //}
644  else
645  {
646  coarsestCorrField = 0;
647  solverPerformance coarseSolverPerf;
648 
649  if (matrixLevels_[coarsestLevel].asymmetric())
650  {
651  coarseSolverPerf = BICCG
652  (
653  "coarsestLevelCorr",
654  matrixLevels_[coarsestLevel],
655  interfaceLevelsBouCoeffs_[coarsestLevel],
656  interfaceLevelsIntCoeffs_[coarsestLevel],
657  interfaceLevels_[coarsestLevel],
658  tolerance_,
659  relTol_
660  ).solve
661  (
662  coarsestCorrField,
663  coarsestSource
664  );
665  }
666  else
667  {
668  coarseSolverPerf = ICCG
669  (
670  "coarsestLevelCorr",
671  matrixLevels_[coarsestLevel],
672  interfaceLevelsBouCoeffs_[coarsestLevel],
673  interfaceLevelsIntCoeffs_[coarsestLevel],
674  interfaceLevels_[coarsestLevel],
675  tolerance_,
676  relTol_
677  ).solve
678  (
679  coarsestCorrField,
680  coarsestSource
681  );
682  }
683 
684  if (debug >= 2)
685  {
686  coarseSolverPerf.print(Info.masterStream(coarseComm));
687  }
688  }
689 
690  UPstream::warnComm = oldWarn;
691 }
692 
693 
694 // ************************************************************************* //
const lduMatrix & matrix_
Definition: lduMatrix.H:98
Incomplete Cholesky preconditioned CG solver derived from the general preconditioned CG solver PCG bu...
Definition: ICCG.H:54
lduInterfaceFieldPtrsList interfaces_
Definition: lduMatrix.H:101
unsigned char direction
Definition: direction.H:43
Pre-declare related SubField type.
Definition: Field.H:61
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:142
Diagonal incomplete LU preconditioned BiCG solver derived from the general preconditioned BiCG solver...
Definition: BICCG.H:54
bool set(const label) const
Is element set.
Definition: PtrListI.H:107
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 solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PBiCG.C:66
void assign(const UList< T > &)
Assign elements to those from UList.
Definition: UList.C:37
const Type & initialResidual() const
Return initial residual.
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics...
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
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.
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PCG.C:66
messageStream Info
tmp< surfaceScalarField > interpolate(const RhoType &rho)
bool checkConvergence(const Type &tolerance, const Type &relTolerance)
Check, store and return convergence.
dynamicFvMesh & mesh
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
scalar relTol_
Convergence tolerance relative to the initial.
Definition: lduMatrix.H:116
const FieldField< Field, scalar > & interfaceIntCoeffs_
Definition: lduMatrix.H:100
void setSize(const label)
Reset size of List.
Definition: List.C:318
const Type & finalResidual() const
Return final residual.
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
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
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
scalarField & diag()
Definition: lduMatrix.C:183
scalar tolerance_
Final convergence tolerance.
Definition: lduMatrix.H:113
const labelField & restrictAddressing(const label leveli) const
Return cell restrict addressing of given level.
dictionary controlDict_
Dictionary of controls.
Definition: lduMatrix.H:104
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve.
label nCells(const label leveli) const
Return number of coarse cells (before processor agglomeration)
scalar gSumMag(const FieldField< Field, Type > &f)
const FieldField< Field, scalar > & interfaceBouCoeffs_
Definition: lduMatrix.H:99
OSstream & masterStream(const label communicator)
Convert to OSstream.
Definition: messageStream.C:60
label nIterations() const
Return number of iterations.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void restrictField(Field< Type > &cf, const Field< Type > &ff, const label fineLevelIndex, const bool procAgglom) const
Restrict (integrate by summation) cell field.
void prolongField(Field< Type > &ff, const Field< Type > &cf, const label coarseLevelIndex, const bool procAgglom) const
Prolong (interpolate by injection) cell field.
void print(Ostream &os) const
Print summary of solver performance to the given stream.
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53