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