src/coupledMatrix/coupledLduPrecon/CholeskyPrecon/coupledCholeskyPrecon.C

   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.
  10
  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.
  15
  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.
  20
  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/>.
  23
  24 Class
  25     coupledCholeskyPrecon
  26
  27 Description
  28     Incomplete Cholesky preconditioning with no fill-in for coupled matrices
  29
  30 Author
  31     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
  32
  33 \*----------------------------------------------------------------------------*/
  34
  35 #include "coupledCholeskyPrecon.H"
  36 #include "addToRunTimeSelectionTable.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42     defineTypeNameAndDebug(coupledCholeskyPrecon, 0);
  43
  44     addToRunTimeSelectionTable
  45     (
  46         coupledLduPrecon,
  47         coupledCholeskyPrecon,
  48         dictionary
  49     );
  50 }
  51
  52
  53 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  54
  55 void Foam::coupledCholeskyPrecon::calcPreconDiag()
  56 {
  57     // Precondition the diagonal
  58     forAll (matrix_, rowI)
  59     {
  60         const lduMatrix& rowMatrix = matrix_[rowI];
  61
  62         preconDiag_.set(rowI, new scalarField(rowMatrix.diag()));
  63         scalarField& rowPreconDiag = preconDiag_[rowI];
  64
  65         if (rowMatrix.symmetric())
  66         {
  67             const unallocLabelList& upperAddr = rowMatrix.lduAddr().upperAddr();
  68             const unallocLabelList& lowerAddr = rowMatrix.lduAddr().lowerAddr();
  69
  70             // Get off-diagonal matrix coefficients
  71             const scalarField& upper = rowMatrix.upper();
  72
  73             forAll (upper, coeffI)
  74             {
  75                 rowPreconDiag[upperAddr[coeffI]] -=
  76                     sqr(upper[coeffI])/rowPreconDiag[lowerAddr[coeffI]];
  77             }
  78         }
  79         else if (rowMatrix.asymmetric())
  80         {
  81             const unallocLabelList& upperAddr = rowMatrix.lduAddr().upperAddr();
  82             const unallocLabelList& lowerAddr = rowMatrix.lduAddr().lowerAddr();
  83
  84             // Get off-diagonal matrix coefficients
  85             const scalarField& upper = rowMatrix.upper();
  86             const scalarField& lower = rowMatrix.lower();
  87
  88             forAll (upper, coeffI)
  89             {
  90                 rowPreconDiag[upperAddr[coeffI]] -=
  91                     upper[coeffI]*lower[coeffI]/
  92                     rowPreconDiag[lowerAddr[coeffI]];
  93             }
  94         }
  95
  96         // Invert the diagonal for future use
  97         forAll (rowPreconDiag, i)
  98         {
  99             rowPreconDiag[i] = 1.0/rowPreconDiag[i];
 100         }
 101     }
 102 }
 103
 104
 105 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 106
 107 Foam::coupledCholeskyPrecon::coupledCholeskyPrecon
 108 (
 109     const coupledLduMatrix& matrix,
 110     const PtrList<FieldField<Field, scalar> >& bouCoeffs,
 111     const PtrList<FieldField<Field, scalar> >& intCoeffs,
 112     const lduInterfaceFieldPtrsListList& interfaces
 113 )
 114 :
 115     coupledLduPrecon
 116     (
 117         matrix,
 118         bouCoeffs,
 119         intCoeffs,
 120         interfaces
 121     ),
 122     preconDiag_(matrix.size())
 123 {
 124     calcPreconDiag();
 125 }
 126
 127
 128 Foam::coupledCholeskyPrecon::coupledCholeskyPrecon
 129 (
 130     const coupledLduMatrix& matrix,
 131     const PtrList<FieldField<Field, scalar> >& bouCoeffs,
 132     const PtrList<FieldField<Field, scalar> >& intCoeffs,
 133     const lduInterfaceFieldPtrsListList& interfaces,
 134     const dictionary& dict
 135 )
 136 :
 137     coupledLduPrecon
 138     (
 139         matrix,
 140         bouCoeffs,
 141         intCoeffs,
 142         interfaces
 143     ),
 144     preconDiag_(matrix.size())
 145 {
 146     calcPreconDiag();
 147 }
 148
 149
 150 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 151
 152 void Foam::coupledCholeskyPrecon::precondition
 153 (
 154     FieldField<Field, scalar>& x,
 155     const FieldField<Field, scalar>& b,
 156     const direction cmpt
 157 ) const
 158 {
 159     // Cholesky precondition all matrices
 160     forAll (matrix_, rowI)
 161     {
 162         scalarField& rowX = x[rowI];
 163         const scalarField& rowB = b[rowI];
 164
 165         const lduMatrix& rowMatrix = matrix_[rowI];
 166         const scalarField& rowPreconDiag = preconDiag_[rowI];
 167
 168         forAll(rowX, i)
 169         {
 170             rowX[i] = rowB[i]*rowPreconDiag[i];
 171         }
 172
 173         if (rowMatrix.symmetric())
 174         {
 175             const unallocLabelList& upperAddr = rowMatrix.lduAddr().upperAddr();
 176             const unallocLabelList& lowerAddr = rowMatrix.lduAddr().lowerAddr();
 177
 178             // Get off-diagonal matrix coefficients
 179             const scalarField& upper = rowMatrix.upper();
 180
 181             forAll (upper, coeffI)
 182             {
 183                 rowX[upperAddr[coeffI]] -=
 184                     rowPreconDiag[upperAddr[coeffI]]*
 185                     upper[coeffI]*rowX[lowerAddr[coeffI]];
 186             }
 187
 188             forAllReverse (upper, coeffI)
 189             {
 190                 rowX[lowerAddr[coeffI]] -=
 191                     rowPreconDiag[lowerAddr[coeffI]]*
 192                     upper[coeffI]*rowX[upperAddr[coeffI]];
 193             }
 194         }
 195         else if (rowMatrix.asymmetric())
 196         {
 197             const unallocLabelList& upperAddr = rowMatrix.lduAddr().upperAddr();
 198             const unallocLabelList& lowerAddr = rowMatrix.lduAddr().lowerAddr();
 199             const unallocLabelList& losortAddr =
 200                 rowMatrix.lduAddr().losortAddr();
 201
 202             // Get off-diagonal matrix coefficients
 203             const scalarField& upper = rowMatrix.upper();
 204             const scalarField& lower = rowMatrix.lower();
 205
 206             label losortCoeff;
 207
 208             forAll (lower, coeffI)
 209             {
 210                 losortCoeff = losortAddr[coeffI];
 211
 212                 rowX[upperAddr[losortCoeff]] -=
 213                     rowPreconDiag[upperAddr[losortCoeff]]*
 214                     lower[losortCoeff]*rowX[lowerAddr[losortCoeff]];
 215             }
 216
 217             forAllReverse (upper, coeffI)
 218             {
 219                 rowX[lowerAddr[coeffI]] -=
 220                     rowPreconDiag[lowerAddr[coeffI]]*
 221                     upper[coeffI]*rowX[upperAddr[coeffI]];
 222             }
 223         }
 224     }
 225 }
 226
 227
 228 void Foam::coupledCholeskyPrecon::preconditionT
 229 (
 230     FieldField<Field, scalar>& x,
 231     const FieldField<Field, scalar>& b,
 232     const direction cmpt
 233 ) const
 234 {
 235     // Cholesky precondition all matrices
 236     forAll (matrix_, rowI)
 237     {
 238         scalarField& rowX = x[rowI];
 239         const scalarField& rowB = b[rowI];
 240
 241         const lduMatrix& rowMatrix = matrix_[rowI];
 242         const scalarField& rowPreconDiag = preconDiag_[rowI];
 243
 244         forAll(rowX, i)
 245         {
 246             rowX[i] = rowB[i]*rowPreconDiag[i];
 247         }
 248
 249         if (rowMatrix.symmetric())
 250         {
 251             const unallocLabelList& upperAddr = rowMatrix.lduAddr().upperAddr();
 252             const unallocLabelList& lowerAddr = rowMatrix.lduAddr().lowerAddr();
 253
 254             // Get off-diagonal matrix coefficients
 255             const scalarField& upper = rowMatrix.upper();
 256
 257             forAll (upper, coeffI)
 258             {
 259                 // For symmetric matrix, there is no change.  HJ, 19/Jan/2009
 260                 rowX[upperAddr[coeffI]] -=
 261                     rowPreconDiag[upperAddr[coeffI]]*
 262                     upper[coeffI]*rowX[lowerAddr[coeffI]];
 263             }
 264
 265             forAllReverse (upper, coeffI)
 266             {
 267                 // For symmetric matrix, there is no change.  HJ, 19/Jan/2009
 268                 rowX[lowerAddr[coeffI]] -=
 269                     rowPreconDiag[lowerAddr[coeffI]]*
 270                     upper[coeffI]*rowX[upperAddr[coeffI]];
 271             }
 272         }
 273         else if (rowMatrix.asymmetric())
 274         {
 275             const unallocLabelList& upperAddr = rowMatrix.lduAddr().upperAddr();
 276             const unallocLabelList& lowerAddr = rowMatrix.lduAddr().lowerAddr();
 277             const unallocLabelList& losortAddr =
 278                 rowMatrix.lduAddr().losortAddr();
 279
 280             // Get off-diagonal matrix coefficients
 281             const scalarField& upper = rowMatrix.upper();
 282             const scalarField& lower = rowMatrix.lower();
 283
 284             label losortCoeff;
 285
 286             forAll (lower, coeffI)
 287             {
 288                 // Transpose multiplication.  HJ, 19/Jan/2009
 289                 rowX[upperAddr[coeffI]] -=
 290                     rowPreconDiag[upperAddr[coeffI]]*
 291                     upper[coeffI]*rowX[lowerAddr[coeffI]];
 292             }
 293
 294             forAllReverse (upper, coeffI)
 295             {
 296                 losortCoeff = losortAddr[coeffI];
 297
 298                 // Transpose multiplication.  HJ, 19/Jan/2009
 299                 rowX[lowerAddr[losortCoeff]] -=
 300                     rowPreconDiag[lowerAddr[losortCoeff]]*
 301                     lower[coeffI]*rowX[upperAddr[losortCoeff]];
 302             }
 303         }
 304     }
 305 }
 306
 307
 308 // ************************************************************************* //