FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

GAMGSolverSolve.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "GAMGSolver.H"
00027 #include <OpenFOAM/ICCG.H>
00028 #include <OpenFOAM/BICCG.H>
00029 #include <OpenFOAM/SubField.H>
00030 
00031 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00032 
00033 Foam::lduMatrix::solverPerformance Foam::GAMGSolver::solve
00034 (
00035     scalarField& psi,
00036     const scalarField& source,
00037     const direction cmpt
00038 ) const
00039 {
00040     // Setup class containing solver performance data
00041     lduMatrix::solverPerformance solverPerf(typeName, fieldName_);
00042 
00043     // Calculate A.psi used to calculate the initial residual
00044     scalarField Apsi(psi.size());
00045     matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
00046 
00047     // Create the storage for the finestCorrection which may be used as a
00048     // temporary in normFactor
00049     scalarField finestCorrection(psi.size());
00050 
00051     // Calculate normalisation factor
00052     scalar normFactor = this->normFactor(psi, source, Apsi, finestCorrection);
00053 
00054     if (debug >= 2)
00055     {
00056         Pout<< "   Normalisation factor = " << normFactor << endl;
00057     }
00058 
00059     // Calculate initial finest-grid residual field
00060     scalarField finestResidual(source - Apsi);
00061 
00062     // Calculate normalised residual for convergence test
00063     solverPerf.initialResidual() = gSumMag(finestResidual)/normFactor;
00064     solverPerf.finalResidual() = solverPerf.initialResidual();
00065 
00066 
00067     // Check convergence, solve if not converged
00068     if (!solverPerf.checkConvergence(tolerance_, relTol_))
00069     {
00070         // Create coarse grid correction fields
00071         PtrList<scalarField> coarseCorrFields;
00072 
00073         // Create coarse grid sources
00074         PtrList<scalarField> coarseSources;
00075 
00076         // Create the smoothers for all levels
00077         PtrList<lduMatrix::smoother> smoothers;
00078 
00079         // Initialise the above data structures
00080         initVcycle(coarseCorrFields, coarseSources, smoothers);
00081 
00082         do
00083         {
00084             Vcycle
00085             (
00086                 smoothers,
00087                 psi,
00088                 source,
00089                 Apsi,
00090                 finestCorrection,
00091                 finestResidual,
00092                 coarseCorrFields,
00093                 coarseSources,
00094                 cmpt
00095             );
00096 
00097             // Calculate finest level residual field
00098             matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
00099             finestResidual = source;
00100             finestResidual -= Apsi;
00101 
00102             solverPerf.finalResidual() = gSumMag(finestResidual)/normFactor;
00103 
00104             if (debug >= 2)
00105             {
00106                 solverPerf.print();
00107             }
00108         } while
00109         (
00110             ++solverPerf.nIterations() < maxIter_
00111          && !(solverPerf.checkConvergence(tolerance_, relTol_))
00112         );
00113     }
00114 
00115     return solverPerf;
00116 }
00117 
00118 
00119 void Foam::GAMGSolver::Vcycle
00120 (
00121     const PtrList<lduMatrix::smoother>& smoothers,
00122     scalarField& psi,
00123     const scalarField& source,
00124     scalarField& Apsi,
00125     scalarField& finestCorrection,
00126     scalarField& finestResidual,
00127     PtrList<scalarField>& coarseCorrFields,
00128     PtrList<scalarField>& coarseSources,
00129     const direction cmpt
00130 ) const
00131 {
00132     //debug = 2;
00133 
00134     const label coarsestLevel = matrixLevels_.size() - 1;
00135 
00136     // Restrict finest grid residual for the next level up
00137     agglomeration_.restrictField(coarseSources[0], finestResidual, 0);
00138 
00139     if (debug >= 2 && nPreSweeps_)
00140     {
00141         Pout<< "Pre-smoothing scaling factors: ";
00142     }
00143 
00144 
00145     // Residual restriction (going to coarser levels)
00146     for (label leveli = 0; leveli < coarsestLevel; leveli++)
00147     {
00148         // If the optional pre-smoothing sweeps are selected
00149         // smooth the coarse-grid field for the restriced source
00150         if (nPreSweeps_)
00151         {
00152             coarseCorrFields[leveli] = 0.0;
00153 
00154             smoothers[leveli + 1].smooth
00155             (
00156                 coarseCorrFields[leveli],
00157                 coarseSources[leveli],
00158                 cmpt,
00159                 nPreSweeps_ + leveli
00160             );
00161 
00162             scalarField::subField ACf
00163             (
00164                 Apsi,
00165                 coarseCorrFields[leveli].size()
00166             );
00167 
00168             // Scale coarse-grid correction field
00169             // but not on the coarsest level because it evaluates to 1
00170             if (scaleCorrection_ && leveli < coarsestLevel - 1)
00171             {
00172                 scalar sf = scalingFactor
00173                 (
00174                     const_cast<scalarField&>(ACf.operator const scalarField&()),
00175                     matrixLevels_[leveli],
00176                     coarseCorrFields[leveli],
00177                     interfaceLevelsBouCoeffs_[leveli],
00178                     interfaceLevels_[leveli],
00179                     coarseSources[leveli],
00180                     cmpt
00181                 );
00182 
00183                 if (debug >= 2)
00184                 {
00185                     Pout<< sf << " ";
00186                 }
00187 
00188                 coarseCorrFields[leveli] *= sf;
00189             }
00190 
00191             // Correct the residual with the new solution
00192             matrixLevels_[leveli].Amul
00193             (
00194                 const_cast<scalarField&>(ACf.operator const scalarField&()),
00195                 coarseCorrFields[leveli],
00196                 interfaceLevelsBouCoeffs_[leveli],
00197                 interfaceLevels_[leveli],
00198                 cmpt
00199             );
00200 
00201             coarseSources[leveli] -= ACf;
00202         }
00203 
00204         // Residual is equal to source
00205         agglomeration_.restrictField
00206         (
00207             coarseSources[leveli + 1],
00208             coarseSources[leveli],
00209             leveli + 1
00210         );
00211     }
00212 
00213     if (debug >= 2 && nPreSweeps_)
00214     {
00215         Pout<< endl;
00216     }
00217 
00218 
00219     // Solve Coarsest level with either an iterative or direct solver
00220     solveCoarsestLevel
00221     (
00222         coarseCorrFields[coarsestLevel],
00223         coarseSources[coarsestLevel]
00224     );
00225 
00226 
00227     if (debug >= 2)
00228     {
00229         Pout<< "Post-smoothing scaling factors: ";
00230     }
00231 
00232     // Smoothing and prolongation of the coarse correction fields
00233     // (going to finer levels)
00234     for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
00235     {
00236         // Create a field for the pre-smoothed correction field
00237         // as a sub-field of the finestCorrection which is not
00238         // currently being used
00239         scalarField::subField preSmoothedCoarseCorrField
00240         (
00241             finestCorrection,
00242             coarseCorrFields[leveli].size()
00243         );
00244 
00245         // Only store the preSmoothedCoarseCorrField is pre-smoothing is used
00246         if (nPreSweeps_)
00247         {
00248             preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
00249         }
00250 
00251         agglomeration_.prolongField
00252         (
00253             coarseCorrFields[leveli],
00254             coarseCorrFields[leveli + 1],
00255             leveli + 1
00256         );
00257 
00258         // Scale coarse-grid correction field
00259         // but not on the coarsest level because it evaluates to 1
00260         if (scaleCorrection_ && leveli < coarsestLevel - 1)
00261         {
00262             // Create A.psi for this coarse level as a sub-field of Apsi
00263             scalarField::subField ACf
00264             (
00265                 Apsi,
00266                 coarseCorrFields[leveli].size()
00267             );
00268 
00269             scalar sf = scalingFactor
00270             (
00271                 const_cast<scalarField&>(ACf.operator const scalarField&()),
00272                 matrixLevels_[leveli],
00273                 coarseCorrFields[leveli],
00274                 interfaceLevelsBouCoeffs_[leveli],
00275                 interfaceLevels_[leveli],
00276                 coarseSources[leveli],
00277                 cmpt
00278             );
00279 
00280 
00281             if (debug >= 2)
00282             {
00283                 Pout<< sf << " ";
00284             }
00285 
00286             coarseCorrFields[leveli] *= sf;
00287         }
00288 
00289         // Only add the preSmoothedCoarseCorrField is pre-smoothing is used
00290         if (nPreSweeps_)
00291         {
00292             coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
00293         }
00294 
00295         smoothers[leveli + 1].smooth
00296         (
00297             coarseCorrFields[leveli],
00298             coarseSources[leveli],
00299             cmpt,
00300             nPostSweeps_ + leveli
00301         );
00302     }
00303 
00304     // Prolong the finest level correction
00305     agglomeration_.prolongField
00306     (
00307         finestCorrection,
00308         coarseCorrFields[0],
00309         0
00310     );
00311 
00312     if (scaleCorrection_)
00313     {
00314         // Calculate finest level scaling factor
00315         scalar fsf = scalingFactor
00316         (
00317             Apsi,
00318             matrix_,
00319             finestCorrection,
00320             interfaceBouCoeffs_,
00321             interfaces_,
00322             finestResidual,
00323             cmpt
00324         );
00325 
00326         if (debug >= 2)
00327         {
00328             Pout<< fsf << endl;
00329         }
00330 
00331         forAll(psi, i)
00332         {
00333             psi[i] += fsf*finestCorrection[i];
00334         }
00335     }
00336     else
00337     {
00338         forAll(psi, i)
00339         {
00340             psi[i] += finestCorrection[i];
00341         }
00342     }
00343 
00344     smoothers[0].smooth
00345     (
00346         psi,
00347         source,
00348         cmpt,
00349         nFinestSweeps_
00350     );
00351 }
00352 
00353 
00354 void Foam::GAMGSolver::initVcycle
00355 (
00356     PtrList<scalarField>& coarseCorrFields,
00357     PtrList<scalarField>& coarseSources,
00358     PtrList<lduMatrix::smoother>& smoothers
00359 ) const
00360 {
00361     coarseCorrFields.setSize(matrixLevels_.size());
00362     coarseSources.setSize(matrixLevels_.size());
00363     smoothers.setSize(matrixLevels_.size() + 1);
00364 
00365     // Create the smoother for the finest level
00366     smoothers.set
00367     (
00368         0,
00369         lduMatrix::smoother::New
00370         (
00371             fieldName_,
00372             matrix_,
00373             interfaceBouCoeffs_,
00374             interfaceIntCoeffs_,
00375             interfaces_,
00376             controlDict_
00377         )
00378     );
00379 
00380     forAll (matrixLevels_, leveli)
00381     {
00382         coarseCorrFields.set
00383         (
00384             leveli,
00385             new scalarField
00386             (
00387                 agglomeration_.meshLevel(leveli + 1).lduAddr().size()
00388             )
00389         );
00390 
00391         coarseSources.set
00392         (
00393             leveli,
00394             new scalarField
00395             (
00396                 agglomeration_.meshLevel(leveli + 1).lduAddr().size()
00397             )
00398         );
00399 
00400         smoothers.set
00401         (
00402             leveli + 1,
00403             lduMatrix::smoother::New
00404             (
00405                 fieldName_,
00406                 matrixLevels_[leveli],
00407                 interfaceLevelsBouCoeffs_[leveli],
00408                 interfaceLevelsIntCoeffs_[leveli],
00409                 interfaceLevels_[leveli],
00410                 controlDict_
00411             )
00412         );
00413     }
00414 }
00415 
00416 
00417 void Foam::GAMGSolver::solveCoarsestLevel
00418 (
00419     scalarField& coarsestCorrField,
00420     const scalarField& coarsestSource
00421 ) const
00422 {
00423     if (directSolveCoarsest_)
00424     {
00425         coarsestCorrField = coarsestSource;
00426         coarsestLUMatrixPtr_->solve(coarsestCorrField);
00427     }
00428     else
00429     {
00430         const label coarsestLevel = matrixLevels_.size() - 1;
00431         coarsestCorrField = 0;
00432         lduMatrix::solverPerformance coarseSolverPerf;
00433 
00434         if (matrixLevels_[coarsestLevel].asymmetric())
00435         {
00436             coarseSolverPerf = BICCG
00437             (
00438                 "coarsestLevelCorr",
00439                 matrixLevels_[coarsestLevel],
00440                 interfaceLevelsBouCoeffs_[coarsestLevel],
00441                 interfaceLevelsIntCoeffs_[coarsestLevel],
00442                 interfaceLevels_[coarsestLevel],
00443                 tolerance_,
00444                 relTol_
00445             ).solve
00446             (
00447                 coarsestCorrField,
00448                 coarsestSource
00449             );
00450         }
00451         else
00452         {
00453             coarseSolverPerf = ICCG
00454             (
00455                 "coarsestLevelCorr",
00456                 matrixLevels_[coarsestLevel],
00457                 interfaceLevelsBouCoeffs_[coarsestLevel],
00458                 interfaceLevelsIntCoeffs_[coarsestLevel],
00459                 interfaceLevels_[coarsestLevel],
00460                 tolerance_,
00461                 relTol_
00462             ).solve
00463             (
00464                 coarsestCorrField,
00465                 coarsestSource
00466             );
00467         }
00468 
00469         if (debug >= 2)
00470         {
00471             coarseSolverPerf.print();
00472         }
00473     }
00474 }
00475 
00476 
00477 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines