src/lduSolvers/lduPrecon/ILUCp/ILUCp.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   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     ILUCp
  26
  27 Description
  28     ILU preconditioning with arbitrary level of fill in (p), based on Crout
  29     algorithm.
  30
  31     Reference: Saad, Y.: Iterative Methods for Sparse Linear Systems (2nd
  32     Edition), SIAM, 2003.
  33
  34 Author
  35     Vuko Vukcevic, FMENA Zagreb. All rights reserved
  36
  37 \*---------------------------------------------------------------------------*/
  38
  39 #include "ILUCp.H"
  40 #include "addToRunTimeSelectionTable.H"
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44 namespace Foam
  45 {
  46     defineTypeNameAndDebug(ILUCp, 0);
  47
  48     lduPreconditioner::
  49         addasymMatrixConstructorToTable<ILUCp>
  50         addILUCpPreconditionerAsymMatrixConstructorToTable_;
  51
  52     // Add to symmetric constructor table as well until we implement Choleskyp
  53     // preconditioner. VV, 27/Jun/2015.
  54     lduPreconditioner::
  55         addsymMatrixConstructorToTable<ILUCp>
  56         addILUCpPreconditionerSymMatrixConstructorToTable_;
  57 }
  58
  59
  60 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  61
  62
  63 void Foam::ILUCp::calcFactorization()
  64 {
  65     if (!matrix_.diagonal())
  66     {
  67         // Get necessary const access to extended ldu addressing
  68         const extendedLduAddressing& addr = extMatrix_.extendedLduAddr();
  69
  70         // Get upper/lower extended addressing
  71         const label* const __restrict__ uPtr = addr.extendedUpperAddr().begin();
  72         const label* const __restrict__ lPtr = addr.extendedLowerAddr().begin();
  73
  74         // Get extended owner start addressing
  75         const label* const __restrict__ ownStartPtr =
  76             addr.extendedOwnerStartAddr().begin();
  77
  78         // Get extended losort and losort start addressing
  79         const label* const __restrict__ lsrPtr =
  80             addr.extendedLosortAddr().begin();
  81         const label* const __restrict__ lsrStartPtr =
  82             addr.extendedLosortStartAddr().begin();
  83
  84         // Get access to factored matrix entries
  85         scalar* __restrict__ diagPtr = preconDiag_.begin();
  86         scalar* __restrict__ upperPtr = extMatrix_.extendedUpper().begin();
  87         scalar* __restrict__ lowerPtr = extMatrix_.extendedLower().begin();
  88
  89         // Get access to working fields
  90         scalar* __restrict__ zPtr = z_.begin();
  91         scalar* __restrict__ wPtr = w_.begin();
  92
  93         // Get number of rows
  94         const label nRows = preconDiag_.size();
  95
  96         // Define start and end face ("virtual" face when extended addressing is
  97         // used) of this row/column, and number of non zero off diagonal entries
  98         register label fStart, fEnd, fLsrStart, fLsrEnd;
  99
 100         // Crout LU factorization
 101
 102         // Row by row loop (k - loop).
 103         for (register label rowI = 0; rowI < nRows; ++rowI)
 104         {
 105             // Start and end of k-th row (upper) and k-th column (lower)
 106             fStart = ownStartPtr[rowI];
 107             fEnd = ownStartPtr[rowI + 1];
 108
 109             // Initialize temporary working diagonal
 110             zDiag_ = diagPtr[rowI];
 111
 112             // Initialize temporary working row field
 113             for (register label faceI = fStart; faceI < fEnd; ++faceI)
 114             {
 115                 // Note: z addressed by neighbour of face (column index for
 116                 // upper), w addressed by neighbour of face (row index for
 117                 // lower)
 118                 zPtr[uPtr[faceI]] = upperPtr[faceI];
 119                 wPtr[uPtr[faceI]] = lowerPtr[faceI];
 120             }
 121
 122             // Start and end of k-th row (lower) and k-th column (upper)
 123             fLsrStart = lsrStartPtr[rowI];
 124             fLsrEnd = lsrStartPtr[rowI + 1];
 125
 126             // Lower coeff loop (first i - loop)
 127             for
 128             (
 129                 register label faceLsrI = fLsrStart;
 130                 faceLsrI < fLsrEnd;
 131                 ++faceLsrI
 132             )
 133             {
 134                 // Get losort coefficient for this face
 135                 const register label losortCoeff = lsrPtr[faceLsrI];
 136
 137                 // Get corresponding row index for upper (i label)
 138                 const label i = lPtr[losortCoeff];
 139
 140                 // Update diagonal
 141                 zDiag_ -= lowerPtr[losortCoeff]*upperPtr[losortCoeff];
 142
 143                 // Get end of row for cell i
 144                 const register label fEndRowi = ownStartPtr[i + 1];
 145
 146                 // Upper coeff loop (additional loop to avoid checking the
 147                 // existence of certain upper coeffs)
 148                 for
 149                 (
 150                     // Diagonal is already updated (losortCoeff + 1 = start)
 151                     register label faceI = losortCoeff + 1;
 152                     faceI < fEndRowi;
 153                     ++faceI
 154                 )
 155                 {
 156                     zPtr[uPtr[faceI]] -= lowerPtr[losortCoeff]*upperPtr[faceI];
 157                     wPtr[uPtr[faceI]] -= upperPtr[losortCoeff]*lowerPtr[faceI];
 158                 }
 159             }
 160
 161             // Update diagonal entry, inverting it for future use
 162             scalar& diagRowI = diagPtr[rowI];
 163             diagRowI = 1.0/zDiag_;
 164
 165             // Index for updating L and U
 166             register label zwIndex;
 167
 168             // Update upper and lower coeffs
 169             for (register label faceI = fStart; faceI < fEnd; ++faceI)
 170             {
 171                 // Get index for current face
 172                 zwIndex = uPtr[faceI];
 173
 174                 // Update L and U decomposition for this row (column)
 175                 upperPtr[faceI] = zPtr[zwIndex];
 176                 lowerPtr[faceI] = wPtr[zwIndex]*diagRowI;
 177             }
 178
 179             // Reset temporary working fields
 180             zDiag_ = 0;
 181
 182             // Only reset parts of the working fields that have been updated in
 183             // this step (for this row and column)
 184             for
 185             (
 186                 register label faceLsrI = fLsrStart;
 187                 faceLsrI < fLsrEnd;
 188                 ++faceLsrI
 189             )
 190             {
 191                 // Get losort coefficient for this face
 192                 const register label losortCoeff = lsrPtr[faceLsrI];
 193
 194                 // Get corresponding row index for upper (i label)
 195                 const label i = lPtr[losortCoeff];
 196
 197                 // Get end of row for cell i
 198                 const register label fEndRowi = ownStartPtr[i + 1];
 199
 200                 for
 201                 (
 202                     register label faceI = losortCoeff + 1;
 203                     faceI < fEndRowi;
 204                     ++faceI
 205                 )
 206                 {
 207                     zPtr[uPtr[faceI]] = 0.0;
 208                     wPtr[uPtr[faceI]] = 0.0;
 209                 }
 210             }
 211         }
 212     }
 213     else
 214     {
 215         forAll (preconDiag_, i)
 216         {
 217             preconDiag_[i] = 1.0/preconDiag_[i];
 218         }
 219     }
 220 }
 221
 222
 223 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 224
 225 Foam::ILUCp::ILUCp
 226 (
 227     const lduMatrix& matrix,
 228     const FieldField<Field, scalar>& coupleBouCoeffs,
 229     const FieldField<Field, scalar>& coupleIntCoeffs,
 230     const lduInterfaceFieldPtrsList& interfaces,
 231     const dictionary& dict
 232 )
 233 :
 234     lduPreconditioner
 235     (
 236         matrix,
 237         coupleBouCoeffs,
 238         coupleIntCoeffs,
 239         interfaces
 240     ),
 241     preconDiag_(matrix_.diag()),
 242     p_(readLabel(dict.lookup("fillInLevel"))),
 243     extMatrix_
 244     (
 245         matrix,
 246         p_,
 247         matrix.mesh().thisDb().parent().lookupObject
 248         <
 249             polyMesh
 250         >(polyMesh::defaultRegion)
 251     ),
 252     zDiag_(0),
 253     z_(preconDiag_.size(), 0),
 254     w_(preconDiag_.size(), 0)
 255 {
 256     calcFactorization();
 257 }
 258
 259
 260 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 261
 262 Foam::ILUCp::~ILUCp()
 263 {}
 264
 265
 266 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 267
 268 void Foam::ILUCp::precondition
 269 (
 270     scalarField& x,
 271     const scalarField& b,
 272     const direction
 273 ) const
 274 {
 275     if (!matrix_.diagonal())
 276     {
 277         // Get matrix addressing
 278         const extendedLduAddressing& addr = extMatrix_.extendedLduAddr();
 279         const unallocLabelList& upperAddr = addr.extendedUpperAddr();
 280         const unallocLabelList& lowerAddr = addr.extendedLowerAddr();
 281         const unallocLabelList& losortAddr = addr.extendedLosortAddr();
 282
 283         // Get upper and lower matrix factors
 284         const scalarField& lower = extMatrix_.extendedLower();
 285         const scalarField& upper = extMatrix_.extendedUpper();
 286
 287         // Solve Lz = b with forward substitution. lower is chosen to be unit
 288         // triangular. z does not need to be stored
 289
 290         // Initialize x field
 291         x = b;
 292
 293         register label losortCoeffI;
 294         register label rowI;
 295
 296         // Forward substitution loop
 297         forAll (lower, coeffI)
 298         {
 299             // Get current losortCoeff to ensure row by row access
 300             losortCoeffI = losortAddr[coeffI];
 301
 302             // Subtract already updated lower part from the solution
 303             x[upperAddr[losortCoeffI]] -=
 304                 lower[losortCoeffI]*x[lowerAddr[losortCoeffI]];
 305         }
 306
 307         // Solve Ux = b with back substitution. U is chosen to be upper
 308         // triangular with diagonal entries corresponding to preconDiag_
 309
 310         // Multiply with inverse diagonal
 311         x *= preconDiag_;
 312
 313         // Back substitution loop
 314         forAllReverse (upper, coeffI)
 315         {
 316             // Get row index
 317             rowI = lowerAddr[coeffI];
 318
 319             // Subtract already updated upper part from the solution
 320             x[rowI] -= upper[coeffI]*x[upperAddr[coeffI]]*preconDiag_[rowI];
 321         }
 322     }
 323     else
 324     {
 325         WarningIn
 326         (
 327             "void ILUCp::precondition"
 328             "(scalarField& x, const scalarField& b, const direction cmpt)"
 329         )   << "Unnecessary use of ILUCp preconditioner for diagonal matrix. "
 330             << nl
 331             << "Use diagonal preconditioner instead."
 332             << endl;
 333
 334         // Diagonal preconditioning
 335         forAll(x, i)
 336         {
 337             x[i] = b[i]*preconDiag_[i];
 338         }
 339     }
 340 }
 341
 342
 343 void Foam::ILUCp::preconditionT
 344 (
 345     scalarField& x,
 346     const scalarField& b,
 347     const direction cmpt
 348 ) const
 349 {
 350     if (!matrix_.diagonal())
 351     {
 352         // Get matrix addressing
 353         const extendedLduAddressing& addr = extMatrix_.extendedLduAddr();
 354         const unallocLabelList& upperAddr = addr.extendedUpperAddr();
 355         const unallocLabelList& lowerAddr = addr.extendedLowerAddr();
 356         const unallocLabelList& losortAddr = addr.extendedLosortAddr();
 357
 358         // Get upper and lower matrix factors
 359         const scalarField& lower = extMatrix_.extendedLower();
 360         const scalarField& upper = extMatrix_.extendedUpper();
 361
 362         // Solve U^T z = b with forward substitution. lower is chosen to
 363         // be unit triangular - U^T (transpose U) "contains" diagonal entries. z
 364         // does not need to be stored.
 365
 366         // Initialize x field
 367         forAll(x, i)
 368         {
 369             x[i] = b[i]*preconDiag_[i];
 370         }
 371
 372         register label losortCoeffI;
 373         register label rowI;
 374
 375         // Forward substitution loop
 376         forAll (upper, coeffI)
 377         {
 378             // Get current losortCoeff to ensure row by row access
 379             losortCoeffI = losortAddr[coeffI];
 380
 381             // Get row index
 382             rowI = upperAddr[losortCoeffI];
 383
 384             // Subtract already updated lower (upper transpose) part from the
 385             // solution
 386             x[rowI] -= upper[losortCoeffI]*x[lowerAddr[losortCoeffI]]*
 387                 preconDiag_[rowI];
 388         }
 389
 390         // Solve L^T x = z with back substitution. L^T is unit upper triangular
 391
 392         // Back substitution loop
 393         forAllReverse (lower, coeffI)
 394         {
 395             // Subtract already updated upper part from the solution
 396             x[lowerAddr[coeffI]] -= lower[coeffI]*x[upperAddr[coeffI]];
 397         }
 398     }
 399     else
 400     {
 401         WarningIn
 402         (
 403             "void ILUCp::preconditionT"
 404             "(scalarField& x, const scalarField& b, const direction cmpt)"
 405         )   << "Unnecessary use of ILUCp preconditioner for diagonal matrix. "
 406             << nl
 407             << "Use diagonal preconditioner instead."
 408             << endl;
 409
 410         // Diagonal preconditioning
 411         forAll(x, i)
 412         {
 413             x[i] = b[i]*preconDiag_[i];
 414         }
 415     }
 416 }
 417
 418
 419 // ************************************************************************* //