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

PCG.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 <OpenFOAM/PCG.H>
00027 
00028 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00029 
00030 namespace Foam
00031 {
00032     defineTypeNameAndDebug(PCG, 0);
00033 
00034     lduMatrix::solver::addsymMatrixConstructorToTable<PCG>
00035         addPCGSymMatrixConstructorToTable_;
00036 }
00037 
00038 
00039 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00040 
00041 Foam::PCG::PCG
00042 (
00043     const word& fieldName,
00044     const lduMatrix& matrix,
00045     const FieldField<Field, scalar>& interfaceBouCoeffs,
00046     const FieldField<Field, scalar>& interfaceIntCoeffs,
00047     const lduInterfaceFieldPtrsList& interfaces,
00048     const dictionary& solverControls
00049 )
00050 :
00051     lduMatrix::solver
00052     (
00053         fieldName,
00054         matrix,
00055         interfaceBouCoeffs,
00056         interfaceIntCoeffs,
00057         interfaces,
00058         solverControls
00059     )
00060 {}
00061 
00062 
00063 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00064 
00065 Foam::lduMatrix::solverPerformance Foam::PCG::solve
00066 (
00067     scalarField& psi,
00068     const scalarField& source,
00069     const direction cmpt
00070 ) const
00071 {
00072     // --- Setup class containing solver performance data
00073     lduMatrix::solverPerformance solverPerf
00074     (
00075         lduMatrix::preconditioner::getName(controlDict_) + typeName,
00076         fieldName_
00077     );
00078 
00079     register label nCells = psi.size();
00080 
00081     scalar* __restrict__ psiPtr = psi.begin();
00082 
00083     scalarField pA(nCells);
00084     scalar* __restrict__ pAPtr = pA.begin();
00085 
00086     scalarField wA(nCells);
00087     scalar* __restrict__ wAPtr = wA.begin();
00088 
00089     scalar wArA = matrix_.great_;
00090     scalar wArAold = wArA;
00091 
00092     // --- Calculate A.psi
00093     matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
00094 
00095     // --- Calculate initial residual field
00096     scalarField rA(source - wA);
00097     scalar* __restrict__ rAPtr = rA.begin();
00098 
00099     // --- Calculate normalisation factor
00100     scalar normFactor = this->normFactor(psi, source, wA, pA);
00101 
00102     if (lduMatrix::debug >= 2)
00103     {
00104         Info<< "   Normalisation factor = " << normFactor << endl;
00105     }
00106 
00107     // --- Calculate normalised residual norm
00108     solverPerf.initialResidual() = gSumMag(rA)/normFactor;
00109     solverPerf.finalResidual() = solverPerf.initialResidual();
00110 
00111     // --- Check convergence, solve if not converged
00112     if (!solverPerf.checkConvergence(tolerance_, relTol_))
00113     {
00114         // --- Select and construct the preconditioner
00115         autoPtr<lduMatrix::preconditioner> preconPtr =
00116         lduMatrix::preconditioner::New
00117         (
00118             *this,
00119             controlDict_
00120         );
00121 
00122         // --- Solver iteration
00123         do
00124         {
00125             // --- Store previous wArA
00126             wArAold = wArA;
00127 
00128             // --- Precondition residual
00129             preconPtr->precondition(wA, rA, cmpt);
00130 
00131             // --- Update search directions:
00132             wArA = gSumProd(wA, rA);
00133 
00134             if (solverPerf.nIterations() == 0)
00135             {
00136                 for (register label cell=0; cell<nCells; cell++)
00137                 {
00138                     pAPtr[cell] = wAPtr[cell];
00139                 }
00140             }
00141             else
00142             {
00143                 scalar beta = wArA/wArAold;
00144 
00145                 for (register label cell=0; cell<nCells; cell++)
00146                 {
00147                     pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
00148                 }
00149             }
00150 
00151 
00152             // --- Update preconditioned residual
00153             matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
00154 
00155             scalar wApA = gSumProd(wA, pA);
00156 
00157 
00158             // --- Test for singularity
00159             if (solverPerf.checkSingularity(mag(wApA)/normFactor)) break;
00160 
00161 
00162             // --- Update solution and residual:
00163 
00164             scalar alpha = wArA/wApA;
00165 
00166             for (register label cell=0; cell<nCells; cell++)
00167             {
00168                 psiPtr[cell] += alpha*pAPtr[cell];
00169                 rAPtr[cell] -= alpha*wAPtr[cell];
00170             }
00171 
00172             solverPerf.finalResidual() = gSumMag(rA)/normFactor;
00173 
00174         } while
00175         (
00176             solverPerf.nIterations()++ < maxIter_
00177         && !(solverPerf.checkConvergence(tolerance_, relTol_))
00178         );
00179     }
00180 
00181     return solverPerf;
00182 }
00183 
00184 
00185 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines