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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2016-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "GAMGSolver.H"
30#include "SubField.H"
31#include "PrecisionAdaptor.H"
32
33// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
34
36(
37 scalarField& psi_s,
38 const scalarField& source,
39 const direction cmpt
40) const
41{
43 solveScalarField& psi = tpsi.ref();
44
46
47 // Setup class containing solver performance data
48 solverPerformance solverPerf(typeName, fieldName_);
49
50 // Calculate A.psi used to calculate the initial residual
53
54 // Create the storage for the finestCorrection which may be used as a
55 // temporary in normFactor
56 solveScalarField finestCorrection(psi.size());
57
58 // Calculate normalisation factor
59 solveScalar normFactor =
60 this->normFactor(psi, tsource(), Apsi, finestCorrection);
61
62 if ((log_ >= 2) || (debug >= 2))
63 {
64 Pout<< " Normalisation factor = " << normFactor << endl;
65 }
66
67 // Calculate initial finest-grid residual field
68 solveScalarField finestResidual(tsource() - Apsi);
69
71 (
74 true
75 );
76
77 // Calculate normalised residual for convergence test
78 solverPerf.initialResidual() = gSumMag
79 (
80 finestResidual,
81 matrix().mesh().comm()
82 )/normFactor;
83 solverPerf.finalResidual() = solverPerf.initialResidual();
84
85
86 // Check convergence, solve if not converged
87 if
88 (
89 minIter_ > 0
90 || !solverPerf.checkConvergence(tolerance_, relTol_, log_)
91 )
92 {
93 // Create coarse grid correction fields
94 PtrList<solveScalarField> coarseCorrFields;
95
96 // Create coarse grid sources
97 PtrList<solveScalarField> coarseSources;
98
99 // Create the smoothers for all levels
101
102 // Scratch fields if processor-agglomerated coarse level meshes
103 // are bigger than original. Usually not needed
104 solveScalarField scratch1;
105 solveScalarField scratch2;
106
107 // Initialise the above data structures
108 initVcycle
109 (
110 coarseCorrFields,
111 coarseSources,
112 smoothers,
113 scratch1,
114 scratch2
115 );
116
117 do
118 {
119 Vcycle
120 (
121 smoothers,
122 psi,
123 source,
124 Apsi,
125 finestCorrection,
126 finestResidual,
127
128 (scratch1.size() ? scratch1 : Apsi),
129 (scratch2.size() ? scratch2 : finestCorrection),
130
131 coarseCorrFields,
132 coarseSources,
133 cmpt
134 );
135
136 // Calculate finest level residual field
138 finestResidual = tsource();
139 finestResidual -= Apsi;
140
141 solverPerf.finalResidual() = gSumMag
142 (
143 finestResidual,
144 matrix().mesh().comm()
145 )/normFactor;
146
147 if ((log_ >= 2) || (debug >= 2))
148 {
149 solverPerf.print(Info.masterStream(matrix().mesh().comm()));
150 }
151 } while
152 (
153 (
154 ++solverPerf.nIterations() < maxIter_
155 && !solverPerf.checkConvergence(tolerance_, relTol_, log_)
156 )
157 || solverPerf.nIterations() < minIter_
158 );
159 }
160
162 (
165 false
166 );
167
168 return solverPerf;
169}
170
171
172void Foam::GAMGSolver::Vcycle
173(
174 const PtrList<lduMatrix::smoother>& smoothers,
176 const scalarField& source,
177 solveScalarField& Apsi,
178 solveScalarField& finestCorrection,
179 solveScalarField& finestResidual,
180
181 solveScalarField& scratch1,
182 solveScalarField& scratch2,
183
184 PtrList<solveScalarField>& coarseCorrFields,
185 PtrList<solveScalarField>& coarseSources,
186 const direction cmpt
187) const
188{
189 //debug = 2;
190
191 const label coarsestLevel = matrixLevels_.size() - 1;
192
193 // Restrict finest grid residual for the next level up.
194 agglomeration_.restrictField(coarseSources[0], finestResidual, 0, true);
195
196 if (nPreSweeps_ && ((log_ >= 2) || (debug >= 2)))
197 {
198 Pout<< "Pre-smoothing scaling factors: ";
199 }
200
201
202 // Residual restriction (going to coarser levels)
203 for (label leveli = 0; leveli < coarsestLevel; leveli++)
204 {
205 if (coarseSources.set(leveli + 1))
206 {
207 // If the optional pre-smoothing sweeps are selected
208 // smooth the coarse-grid field for the restricted source
209 if (nPreSweeps_)
210 {
211 coarseCorrFields[leveli] = 0.0;
212
213 smoothers[leveli + 1].scalarSmooth
214 (
215 coarseCorrFields[leveli],
216 coarseSources[leveli], //coarseSource,
217 cmpt,
218 min
219 (
220 nPreSweeps_ + preSweepsLevelMultiplier_*leveli,
221 maxPreSweeps_
222 )
223 );
224
226 (
227 scratch1,
228 coarseCorrFields[leveli].size()
229 );
230
231 // Scale coarse-grid correction field
232 // but not on the coarsest level because it evaluates to 1
233 if (scaleCorrection_ && leveli < coarsestLevel - 1)
234 {
235 scale
236 (
237 coarseCorrFields[leveli],
238 const_cast<solveScalarField&>
239 (
240 ACf.operator const solveScalarField&()
241 ),
242 matrixLevels_[leveli],
243 interfaceLevelsBouCoeffs_[leveli],
244 interfaceLevels_[leveli],
245 coarseSources[leveli],
246 cmpt
247 );
248 }
249
250 // Correct the residual with the new solution
251 matrixLevels_[leveli].Amul
252 (
253 const_cast<solveScalarField&>
254 (
255 ACf.operator const solveScalarField&()
256 ),
257 coarseCorrFields[leveli],
258 interfaceLevelsBouCoeffs_[leveli],
259 interfaceLevels_[leveli],
260 cmpt
261 );
262
263 coarseSources[leveli] -= ACf;
264 }
265
266 // Residual is equal to source
267 agglomeration_.restrictField
268 (
269 coarseSources[leveli + 1],
270 coarseSources[leveli],
271 leveli + 1,
272 true
273 );
274 }
275 }
276
277 if (nPreSweeps_ && ((log_ >= 2) || (debug >= 2)))
278 {
279 Pout<< endl;
280 }
281
282
283 // Solve Coarsest level with either an iterative or direct solver
284 if (coarseCorrFields.set(coarsestLevel))
285 {
286 solveCoarsestLevel
287 (
288 coarseCorrFields[coarsestLevel],
289 coarseSources[coarsestLevel]
290 );
291 }
292
293 if ((log_ >= 2) || (debug >= 2))
294 {
295 Pout<< "Post-smoothing scaling factors: ";
296 }
297
298 // Smoothing and prolongation of the coarse correction fields
299 // (going to finer levels)
300
301 solveScalarField dummyField(0);
302
303 for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
304 {
305 if (coarseCorrFields.set(leveli))
306 {
307 // Create a field for the pre-smoothed correction field
308 // as a sub-field of the finestCorrection which is not
309 // currently being used
310 solveScalarField::subField preSmoothedCoarseCorrField
311 (
312 scratch2,
313 coarseCorrFields[leveli].size()
314 );
315
316 // Only store the preSmoothedCoarseCorrField if pre-smoothing is
317 // used
318 if (nPreSweeps_)
319 {
320 preSmoothedCoarseCorrField = coarseCorrFields[leveli];
321 }
322
323 agglomeration_.prolongField
324 (
325 coarseCorrFields[leveli],
326 (
327 coarseCorrFields.set(leveli + 1)
328 ? coarseCorrFields[leveli + 1]
329 : dummyField // dummy value
330 ),
331 leveli + 1,
332 true
333 );
334
335
336 // Create A.psi for this coarse level as a sub-field of Apsi
338 (
339 scratch1,
340 coarseCorrFields[leveli].size()
341 );
342 solveScalarField& ACfRef =
343 const_cast
344 <
346 >(ACf.operator const solveScalarField&());
347
348 if (interpolateCorrection_) //&& leveli < coarsestLevel - 2)
349 {
350 if (coarseCorrFields.set(leveli+1))
351 {
353 (
354 coarseCorrFields[leveli],
355 ACfRef,
356 matrixLevels_[leveli],
357 interfaceLevelsBouCoeffs_[leveli],
358 interfaceLevels_[leveli],
359 agglomeration_.restrictAddressing(leveli + 1),
360 coarseCorrFields[leveli + 1],
361 cmpt
362 );
363 }
364 else
365 {
367 (
368 coarseCorrFields[leveli],
369 ACfRef,
370 matrixLevels_[leveli],
371 interfaceLevelsBouCoeffs_[leveli],
372 interfaceLevels_[leveli],
373 cmpt
374 );
375 }
376 }
377
378 // Scale coarse-grid correction field
379 // but not on the coarsest level because it evaluates to 1
380 if
381 (
382 scaleCorrection_
383 && (interpolateCorrection_ || leveli < coarsestLevel - 1)
384 )
385 {
386 scale
387 (
388 coarseCorrFields[leveli],
389 ACfRef,
390 matrixLevels_[leveli],
391 interfaceLevelsBouCoeffs_[leveli],
392 interfaceLevels_[leveli],
393 coarseSources[leveli],
394 cmpt
395 );
396 }
397
398 // Only add the preSmoothedCoarseCorrField if pre-smoothing is
399 // used
400 if (nPreSweeps_)
401 {
402 coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
403 }
404
405 smoothers[leveli + 1].scalarSmooth
406 (
407 coarseCorrFields[leveli],
408 coarseSources[leveli], //coarseSource,
409 cmpt,
410 min
411 (
412 nPostSweeps_ + postSweepsLevelMultiplier_*leveli,
413 maxPostSweeps_
414 )
415 );
416 }
417 }
418
419 // Prolong the finest level correction
420 agglomeration_.prolongField
421 (
422 finestCorrection,
423 coarseCorrFields[0],
424 0,
425 true
426 );
427
428 if (interpolateCorrection_)
429 {
431 (
432 finestCorrection,
433 Apsi,
434 matrix_,
435 interfaceBouCoeffs_,
436 interfaces_,
437 agglomeration_.restrictAddressing(0),
438 coarseCorrFields[0],
439 cmpt
440 );
441 }
442
443 if (scaleCorrection_)
444 {
445 // Scale the finest level correction
446 scale
447 (
448 finestCorrection,
449 Apsi,
450 matrix_,
451 interfaceBouCoeffs_,
452 interfaces_,
453 finestResidual,
454 cmpt
455 );
456 }
457
458 forAll(psi, i)
459 {
460 psi[i] += finestCorrection[i];
461 }
462
463 smoothers[0].smooth
464 (
465 psi,
466 source,
467 cmpt,
468 nFinestSweeps_
469 );
470}
471
472
473void Foam::GAMGSolver::initVcycle
474(
475 PtrList<solveScalarField>& coarseCorrFields,
476 PtrList<solveScalarField>& coarseSources,
477 PtrList<lduMatrix::smoother>& smoothers,
478 solveScalarField& scratch1,
479 solveScalarField& scratch2
480) const
481{
482 label maxSize = matrix_.diag().size();
483
484 coarseCorrFields.setSize(matrixLevels_.size());
485 coarseSources.setSize(matrixLevels_.size());
486 smoothers.setSize(matrixLevels_.size() + 1);
487
488 // Create the smoother for the finest level
489 smoothers.set
490 (
491 0,
493 (
494 fieldName_,
495 matrix_,
496 interfaceBouCoeffs_,
497 interfaceIntCoeffs_,
498 interfaces_,
499 controlDict_
500 )
501 );
502
503 forAll(matrixLevels_, leveli)
504 {
505 if (agglomeration_.nCells(leveli) >= 0)
506 {
507 label nCoarseCells = agglomeration_.nCells(leveli);
508
509 coarseSources.set(leveli, new solveScalarField(nCoarseCells));
510 }
511
512 if (matrixLevels_.set(leveli))
513 {
514 const lduMatrix& mat = matrixLevels_[leveli];
515
516 label nCoarseCells = mat.diag().size();
517
518 maxSize = max(maxSize, nCoarseCells);
519
520 coarseCorrFields.set(leveli, new solveScalarField(nCoarseCells));
521
522 smoothers.set
523 (
524 leveli + 1,
526 (
527 fieldName_,
528 matrixLevels_[leveli],
529 interfaceLevelsBouCoeffs_[leveli],
530 interfaceLevelsIntCoeffs_[leveli],
531 interfaceLevels_[leveli],
532 controlDict_
533 )
534 );
535 }
536 }
537
538 if (maxSize > matrix_.diag().size())
539 {
540 // Allocate some scratch storage
541 scratch1.setSize(maxSize);
542 scratch2.setSize(maxSize);
543 }
544}
545
546
547Foam::dictionary Foam::GAMGSolver::PCGsolverDict
548(
549 const scalar tol,
550 const scalar relTol
551) const
552{
553 dictionary dict(IStringStream("solver PCG; preconditioner DIC;")());
554 dict.add("tolerance", tol);
555 dict.add("relTol", relTol);
556
557 return dict;
558}
559
560
561Foam::dictionary Foam::GAMGSolver::PBiCGStabSolverDict
562(
563 const scalar tol,
564 const scalar relTol
565) const
566{
567 dictionary dict(IStringStream("solver PBiCGStab; preconditioner DILU;")());
568 dict.add("tolerance", tol);
569 dict.add("relTol", relTol);
570
571 return dict;
572}
573
574
575void Foam::GAMGSolver::solveCoarsestLevel
576(
577 solveScalarField& coarsestCorrField,
578 const solveScalarField& coarsestSource
579) const
580{
581 const label coarsestLevel = matrixLevels_.size() - 1;
582
583 const label coarseComm = matrixLevels_[coarsestLevel].mesh().comm();
584
585 if (directSolveCoarsest_)
586 {
587 PrecisionAdaptor<scalar, solveScalar> tcorrField(coarsestCorrField);
588
589 coarsestLUMatrixPtr_->solve
590 (
591 tcorrField.ref(),
592 ConstPrecisionAdaptor<scalar, solveScalar>(coarsestSource)()
593 );
594 }
595 //else if
596 //(
597 // agglomeration_.processorAgglomerate()
598 // && procMatrixLevels_.set(coarsestLevel)
599 //)
600 //{
601 // //const labelList& agglomProcIDs = agglomeration_.agglomProcIDs
602 // //(
603 // // coarsestLevel
604 // //);
605 // //
606 // //scalarField allSource;
607 // //
608 // //globalIndex cellOffsets;
609 // //if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
610 // //{
611 // // cellOffsets.offsets() =
612 // // agglomeration_.cellOffsets(coarsestLevel);
613 // //}
614 // //
615 // //cellOffsets.gather
616 // //(
617 // // coarseComm,
618 // // agglomProcIDs,
619 // // coarsestSource,
620 // // allSource
621 // //);
622 // //
623 // //scalarField allCorrField;
624 // //solverPerformance coarseSolverPerf;
625 //
626 // label solveComm = agglomeration_.procCommunicator(coarsestLevel);
627 //
628 // coarsestCorrField = 0;
629 // solverPerformance coarseSolverPerf;
630 //
631 // if (Pstream::myProcNo(solveComm) != -1)
632 // {
633 // const lduMatrix& allMatrix = procMatrixLevels_[coarsestLevel];
634 //
635 // {
636 // Pout<< "** Master:Solving on comm:" << solveComm
637 // << " with procs:" << UPstream::procID(solveComm) << endl;
638 //
639 // if (allMatrix.asymmetric())
640 // {
641 // coarseSolverPerf = PBiCGStab
642 // (
643 // "coarsestLevelCorr",
644 // allMatrix,
645 // procInterfaceLevelsBouCoeffs_[coarsestLevel],
646 // procInterfaceLevelsIntCoeffs_[coarsestLevel],
647 // procInterfaceLevels_[coarsestLevel],
648 // PBiCGStabSolverDict(tolerance_, relTol_)
649 // ).solve
650 // (
651 // coarsestCorrField,
652 // coarsestSource
653 // );
654 // }
655 // else
656 // {
657 // coarseSolverPerf = PCG
658 // (
659 // "coarsestLevelCorr",
660 // allMatrix,
661 // procInterfaceLevelsBouCoeffs_[coarsestLevel],
662 // procInterfaceLevelsIntCoeffs_[coarsestLevel],
663 // procInterfaceLevels_[coarsestLevel],
664 // PCGsolverDict(tolerance_, relTol_)
665 // ).solve
666 // (
667 // coarsestCorrField,
668 // coarsestSource
669 // );
670 // }
671 // }
672 // }
673 //
674 // Pout<< "done master solve." << endl;
675 //
676 // //// Scatter to all processors
677 // //coarsestCorrField.setSize(coarsestSource.size());
678 // //cellOffsets.scatter
679 // //(
680 // // coarseComm,
681 // // agglomProcIDs,
682 // // allCorrField,
683 // // coarsestCorrField
684 // //);
685 //
686 // if (debug >= 2)
687 // {
688 // coarseSolverPerf.print(Info.masterStream(coarseComm));
689 // }
690 //
691 // Pout<< "procAgglom: coarsestSource :" << coarsestSource << endl;
692 // Pout<< "procAgglom: coarsestCorrField:" << coarsestCorrField << endl;
693 //}
694 else
695 {
696 coarsestCorrField = 0;
697 const solverPerformance coarseSolverPerf
698 (
699 coarsestSolverPtr_->scalarSolve
700 (
701 coarsestCorrField,
702 coarsestSource
703 )
704 );
705
706 if ((log_ >= 2) || debug)
707 {
708 coarseSolverPerf.print(Info.masterStream(coarseComm));
709 }
710 }
711}
712
713
714// ************************************************************************* //
A const Field/List wrapper with possible data conversion.
SubField< solveScalar > subField
Declare type of subField.
Definition: Field.H:89
A non-const Field/List wrapper with possible data conversion.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
const T * set(const label i) const
Definition: PtrList.H:138
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics.
const Type & finalResidual() const noexcept
Return final residual.
const labelType & nIterations() const noexcept
Return number of iterations.
void print(Ostream &os) const
Print summary of solver performance to the given stream.
const Type & initialResidual() const noexcept
Return initial residual.
bool checkConvergence(const Type &tolerance, const Type &relTolerance, const int logLevel=0)
Check, store and return convergence.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
entry * add(entry *entryPtr, bool mergeEntry=false)
Add a new entry.
Definition: dictionary.C:640
const lduMatrix & matrix_
Definition: lduMatrix.H:108
label maxIter_
Maximum number of iterations in the solver.
Definition: lduMatrix.H:123
scalar tolerance_
Final convergence tolerance.
Definition: lduMatrix.H:126
lduInterfaceFieldPtrsList interfaces_
Definition: lduMatrix.H:111
const FieldField< Field, scalar > & interfaceBouCoeffs_
Definition: lduMatrix.H:109
label minIter_
Minimum number of iterations in the solver.
Definition: lduMatrix.H:120
const lduMatrix & matrix() const noexcept
Definition: lduMatrix.H:233
solveScalarField::cmptType normFactor(const solveScalarField &psi, const solveScalarField &source, const solveScalarField &Apsi, solveScalarField &tmpField) const
int log_
Level of verbosity in the solver output statements.
Definition: lduMatrix.H:117
scalar relTol_
Convergence tolerance relative to the initial.
Definition: lduMatrix.H:129
void setResidualField(const scalarField &residual, const word &fieldName, const bool initial) const
Definition: lduMatrix.C:322
void Amul(solveScalarField &, const tmp< solveScalarField > &, const FieldField< Field, scalar > &, const lduInterfaceFieldPtrsList &, const direction cmpt) const
Matrix multiplication with updated interfaces.
OSstream & masterStream(const label communicator)
T & ref() const
Definition: refPtrI.H:203
const volScalarField & psi
dynamicFvMesh & mesh
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
messageStream Info
Information stream (stdout output on master, null elsewhere)
SolverPerformance< scalar > solverPerformance
SolverPerformance instantiated for a scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Field< solveScalar > solveScalarField
uint8_t direction
Definition: direction.H:56
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
Definition: curveTools.C:75
typeOfMag< Type >::type gSumMag(const FieldField< Field, Type > &f)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dictionary dict
CEqn solve()
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333