Merge branch 'master' of ssh://git.code.sf.net/p/foam-extend/foam-extend-3.2
[foam-extend-3.2.git] / src / lduSolvers / lduPrecon / CholeskyPrecon / CholeskyPrecon.C
blob79de904e36985283e10537e6c6e6f8be496fa96d
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | foam-extend: Open Source CFD
4    \\    /   O peration     | Version:     3.2
5     \\  /    A nd           | Web:         http://www.foam-extend.org
6      \\/     M anipulation  | For copyright notice see file Copyright
7 -------------------------------------------------------------------------------
8 License
9     This file is part of foam-extend.
11     foam-extend is free software: you can redistribute it and/or modify it
12     under the terms of the GNU General Public License as published by the
13     Free Software Foundation, either version 3 of the License, or (at your
14     option) any later version.
16     foam-extend is distributed in the hope that it will be useful, but
17     WITHOUT ANY WARRANTY; without even the implied warranty of
18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19     General Public License for more details.
21     You should have received a copy of the GNU General Public License
22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
24 Class
25     CholeskyPrecon
27 Description
28     Incmplete Cholesky preconditioning with no fill-in
30 Author
31     Hrvoje Jasak, Wikki Ltd.  All rights reserved
33 \*---------------------------------------------------------------------------*/
35 #include "CholeskyPrecon.H"
36 #include "addToRunTimeSelectionTable.H"
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 namespace Foam
42     defineTypeNameAndDebug(CholeskyPrecon, 0);
44     lduPreconditioner::
45         addsymMatrixConstructorToTable<CholeskyPrecon>
46         addCholeskyPreconditionerSymMatrixConstructorToTable_;
50 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
52 void Foam::CholeskyPrecon::calcPreconDiag()
54     // Precondition the diagonal
55     if (matrix_.symmetric())
56     {
57         const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
58         const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
60         // Get off-diagonal matrix coefficients
61         const scalarField& upper = matrix_.upper();
63         forAll (upper, coeffI)
64         {
65             preconDiag_[upperAddr[coeffI]] -=
66                 sqr(upper[coeffI])/preconDiag_[lowerAddr[coeffI]];
67         }
68     }
70     // Invert the diagonal for future use
71     forAll (preconDiag_, i)
72     {
73         preconDiag_[i] = 1.0/preconDiag_[i];
74     }
78 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
80 Foam::CholeskyPrecon::CholeskyPrecon
82     const lduMatrix& matrix,
83     const FieldField<Field, scalar>& coupleBouCoeffs,
84     const FieldField<Field, scalar>& coupleIntCoeffs,
85     const lduInterfaceFieldPtrsList& interfaces,
86     const dictionary& dict
89     lduPreconditioner
90     (
91         matrix,
92         coupleBouCoeffs,
93         coupleIntCoeffs,
94         interfaces
95     ),
96     preconDiag_(matrix_.diag())
98     calcPreconDiag();
102 Foam::CholeskyPrecon::CholeskyPrecon
104     const lduMatrix& matrix,
105     const FieldField<Field, scalar>& coupleBouCoeffs,
106     const FieldField<Field, scalar>& coupleIntCoeffs,
107     const lduInterfaceFieldPtrsList& interfaces
110     lduPreconditioner
111     (
112         matrix,
113         coupleBouCoeffs,
114         coupleIntCoeffs,
115         interfaces
116     ),
117     preconDiag_(matrix_.diag())
119     calcPreconDiag();
123 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
125 Foam::CholeskyPrecon::~CholeskyPrecon()
129 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
131 void Foam::CholeskyPrecon::precondition
133     scalarField& x,
134     const scalarField& b,
135     const direction
136 ) const
138     forAll(x, i)
139     {
140         x[i] = b[i]*preconDiag_[i];
141     }
143     if (matrix_.symmetric())
144     {
145         const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
146         const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
148         // Get off-diagonal matrix coefficients
149         const scalarField& upper = matrix_.upper();
151         forAll (upper, coeffI)
152         {
153             x[upperAddr[coeffI]] -=
154                 preconDiag_[upperAddr[coeffI]]*
155                 upper[coeffI]*x[lowerAddr[coeffI]];
156         }
158         forAllReverse (upper, coeffI)
159         {
160             x[lowerAddr[coeffI]] -=
161                 preconDiag_[lowerAddr[coeffI]]*
162                 upper[coeffI]*x[upperAddr[coeffI]];
163         }
164     }
168 // ************************************************************************* //