src/coupledMatrix/coupledLduPrecon/GaussSeidelPrecon/coupledGaussSeidelPrecon.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     coupledGaussSeidelPrecon
  26
  27 Description
  28     GaussSeidel preconditioning
  29
  30 Author
  31     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "coupledGaussSeidelPrecon.H"
  36 #include "addToRunTimeSelectionTable.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42     defineTypeNameAndDebug(coupledGaussSeidelPrecon, 0);
  43
  44     addToRunTimeSelectionTable
  45     (
  46         coupledLduPrecon,
  47         coupledGaussSeidelPrecon,
  48         dictionary
  49     );
  50 }
  51
  52
  53 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  54
  55 void Foam::coupledGaussSeidelPrecon::forwardSweep
  56 (
  57     const lduMatrix& matrix,
  58     scalarField& x,
  59     scalarField& bPrime
  60 ) const
  61 {
  62     const scalarField& diag = matrix.diag();
  63     const scalarField& lower = matrix.lower();
  64     const scalarField& upper = matrix.upper();
  65
  66     const labelList& upperAddr = matrix.lduAddr().upperAddr();
  67     const labelList& ownStartAddr = matrix.lduAddr().ownerStartAddr();
  68
  69     const label nRows = x.size();
  70     label fStart, fEnd;
  71
  72     for (register label rowI = 0; rowI < nRows; rowI++)
  73     {
  74         // lRow is equal to rowI
  75         scalar& curX = x[rowI];
  76
  77         // Grab the accumulated neighbour side
  78         curX = bPrime[rowI];
  79
  80         // Start and end of this row
  81         fStart = ownStartAddr[rowI];
  82         fEnd = ownStartAddr[rowI + 1];
  83
  84         // Accumulate the owner product side
  85         for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
  86         {
  87             curX -= upper[curCoeff]*x[upperAddr[curCoeff]];
  88         }
  89
  90         // Finish current x
  91         curX /= diag[rowI];
  92
  93         // Distribute the neighbour side using current x
  94         for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
  95         {
  96             bPrime[upperAddr[curCoeff]] -= lower[curCoeff]*curX;
  97         }
  98     }
  99 }
 100
 101
 102 void Foam::coupledGaussSeidelPrecon::reverseSweep
 103 (
 104     const lduMatrix& matrix,
 105     scalarField& x,
 106     scalarField& bPrime
 107 ) const
 108 {
 109     const scalarField& diag = matrix.diag();
 110     const scalarField& lower = matrix.lower();
 111     const scalarField& upper = matrix.upper();
 112
 113     const labelList& upperAddr = matrix.lduAddr().upperAddr();
 114     const labelList& ownStartAddr = matrix.lduAddr().ownerStartAddr();
 115
 116     const label nRows = x.size();
 117     label fStart, fEnd;
 118
 119     for (register label rowI = nRows - 1; rowI >= 0; rowI--)
 120     {
 121         // lRow is equal to rowI
 122         scalar& curX = x[rowI];
 123
 124         // Grab the accumulated neighbour side
 125         curX = bPrime[rowI];
 126
 127         // Start and end of this row
 128         fStart = ownStartAddr[rowI];
 129         fEnd = ownStartAddr[rowI + 1];
 130
 131         // Accumulate the owner product side
 132         for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 133         {
 134             curX -= upper[curCoeff]*x[upperAddr[curCoeff]];
 135         }
 136
 137         // Finish current x
 138         curX /= diag[rowI];
 139
 140         // Distribute the neighbour side using current x
 141         for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 142         {
 143             bPrime[upperAddr[curCoeff]] -= lower[curCoeff]*curX;
 144         }
 145     }
 146 }
 147
 148
 149 void Foam::coupledGaussSeidelPrecon::forwardSweepTranspose
 150 (
 151     const lduMatrix& matrix,
 152     scalarField& x,
 153     scalarField& bPrime
 154 ) const
 155 {
 156     const scalarField& diag = matrix.diag();
 157     const scalarField& lower = matrix.lower();
 158     const scalarField& upper = matrix.upper();
 159
 160     const labelList& upperAddr = matrix.lduAddr().upperAddr();
 161     const labelList& ownStartAddr = matrix.lduAddr().ownerStartAddr();
 162
 163     const label nRows = x.size();
 164     label fStart, fEnd;
 165
 166     for (register label rowI = 0; rowI < nRows; rowI++)
 167     {
 168         // lRow is equal to rowI
 169         scalar& curX = x[rowI];
 170
 171         // Grab the accumulated neighbour side
 172         curX = bPrime[rowI];
 173
 174         // Start and end of this row
 175         fStart = ownStartAddr[rowI];
 176         fEnd = ownStartAddr[rowI + 1];
 177
 178         // Accumulate the owner product side
 179         for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 180         {
 181             // Transpose multiplication.  HJ, 19/Jan/2009
 182             curX -= lower[curCoeff]*x[upperAddr[curCoeff]];
 183         }
 184
 185         // Finish current x
 186         curX /= diag[rowI];
 187
 188         // Distribute the neighbour side using current x
 189         for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 190         {
 191             // Transpose multiplication.  HJ, 19/Jan/2009
 192             bPrime[upperAddr[curCoeff]] -= upper[curCoeff]*curX;
 193         }
 194     }
 195 }
 196
 197
 198 void Foam::coupledGaussSeidelPrecon::reverseSweepTranspose
 199 (
 200     const lduMatrix& matrix,
 201     scalarField& x,
 202     scalarField& bPrime
 203 ) const
 204 {
 205     const scalarField& diag = matrix.diag();
 206     const scalarField& lower = matrix.lower();
 207     const scalarField& upper = matrix.upper();
 208
 209     const labelList& upperAddr = matrix.lduAddr().upperAddr();
 210     const labelList& ownStartAddr = matrix.lduAddr().ownerStartAddr();
 211
 212     const label nRows = x.size();
 213     label fStart, fEnd;
 214
 215     for (register label rowI = nRows - 1; rowI >= 0; rowI--)
 216     {
 217         // lRow is equal to rowI
 218         scalar& curX = x[rowI];
 219
 220         // Grab the accumulated neighbour side
 221         curX = bPrime[rowI];
 222
 223         // Start and end of this row
 224         fStart = ownStartAddr[rowI];
 225         fEnd = ownStartAddr[rowI + 1];
 226
 227         // Accumulate the owner product side
 228         for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 229         {
 230             // Transpose multiplication.  HJ, 19/Jan/2009
 231             curX -= lower[curCoeff]*x[upperAddr[curCoeff]];
 232         }
 233
 234         // Finish current x
 235         curX /= diag[rowI];
 236
 237         // Distribute the neighbour side using current x
 238         for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 239         {
 240             // Transpose multiplication.  HJ, 19/Jan/2009
 241             bPrime[upperAddr[curCoeff]] -= upper[curCoeff]*curX;
 242         }
 243     }
 244 }
 245
 246
 247 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 248
 249 Foam::coupledGaussSeidelPrecon::coupledGaussSeidelPrecon
 250 (
 251     const coupledLduMatrix& matrix,
 252     const PtrList<FieldField<Field, scalar> >& bouCoeffs,
 253     const PtrList<FieldField<Field, scalar> >& intCoeffs,
 254     const lduInterfaceFieldPtrsListList& interfaces
 255 )
 256 :
 257     coupledLduPrecon
 258     (
 259         matrix,
 260         bouCoeffs,
 261         intCoeffs,
 262         interfaces
 263     ),
 264     mBouCoeffs_(bouCoeffs),
 265     bPrime_(matrix.size())
 266 {
 267     // Invert boundary coefficients
 268     forAll (mBouCoeffs_, rowI)
 269     {
 270         mBouCoeffs_[rowI] *= -1;
 271     }
 272
 273     // Hook bPrime components
 274     forAll (matrix_, rowI)
 275     {
 276         bPrime_.set(rowI, new scalarField(matrix_[rowI].lduAddr().size(), 0));
 277     }
 278 }
 279
 280
 281 Foam::coupledGaussSeidelPrecon::coupledGaussSeidelPrecon
 282 (
 283     const coupledLduMatrix& matrix,
 284     const PtrList<FieldField<Field, scalar> >& bouCoeffs,
 285     const PtrList<FieldField<Field, scalar> >& intCoeffs,
 286     const lduInterfaceFieldPtrsListList& interfaces,
 287     const dictionary& dict
 288 )
 289 :
 290     coupledLduPrecon
 291     (
 292         matrix,
 293         bouCoeffs,
 294         intCoeffs,
 295         interfaces
 296     ),
 297     mBouCoeffs_(bouCoeffs),
 298     bPrime_(matrix.size())
 299 {
 300     // Invert boundary coefficients
 301     forAll (mBouCoeffs_, rowI)
 302     {
 303         mBouCoeffs_[rowI] *= -1;
 304     }
 305
 306     // Hook bPrime components
 307     forAll (matrix_, rowI)
 308     {
 309         bPrime_.set(rowI, new scalarField(matrix_[rowI].lduAddr().size(), 0));
 310     }
 311 }
 312
 313
 314 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 315
 316 void Foam::coupledGaussSeidelPrecon::precondition
 317 (
 318     FieldField<Field, scalar>& x,
 319     const FieldField<Field, scalar>& b,
 320     const direction cmpt
 321 ) const
 322 {
 323     // Execute preconditioning
 324     if (matrix_.diagonal())
 325     {
 326         forAll (matrix_, rowI)
 327         {
 328             x[rowI] = b[rowI]/matrix_[rowI].diag();
 329         }
 330     }
 331     else if (matrix_.symmetric() || matrix_.asymmetric())
 332     {
 333         bPrime_ = b;
 334
 335         // Parallel boundary update
 336         {
 337             matrix_.initMatrixInterfaces
 338             (
 339                 mBouCoeffs_,
 340                 interfaces_,
 341                 x,
 342                 bPrime_,
 343                 cmpt,
 344                 true         // switch to lhs
 345             );
 346
 347             matrix_.updateMatrixInterfaces
 348             (
 349                 mBouCoeffs_,
 350                 interfaces_,
 351                 x,
 352                 bPrime_,
 353                 cmpt,
 354                 true         // switch to lhs
 355             );
 356         }
 357
 358         // Forward sweep
 359         forAll (matrix_, rowI)
 360         {
 361             forwardSweep(matrix_[rowI], x[rowI], bPrime_[rowI]);
 362         }
 363
 364         // Reverse sweep
 365         forAllReverse (matrix_, rowI)
 366         {
 367             reverseSweep(matrix_[rowI], x[rowI], bPrime_[rowI]);
 368         }
 369     }
 370 }
 371
 372
 373 void Foam::coupledGaussSeidelPrecon::preconditionT
 374 (
 375     FieldField<Field, scalar>& x,
 376     const FieldField<Field, scalar>& b,
 377     const direction cmpt
 378 ) const
 379 {
 380     // Execute preconditioning
 381     if (matrix_.diagonal())
 382     {
 383         forAll (matrix_, rowI)
 384         {
 385             x[rowI] = b[rowI]/matrix_[rowI].diag();
 386         }
 387     }
 388     else if (matrix_.symmetric() || matrix_.asymmetric())
 389     {
 390         bPrime_ = b;
 391
 392         // Parallel boundary update
 393         {
 394             matrix_.initMatrixInterfaces
 395             (
 396                 mBouCoeffs_,
 397                 interfaces_,
 398                 x,
 399                 bPrime_,
 400                 cmpt
 401             );
 402
 403             matrix_.updateMatrixInterfaces
 404             (
 405                 mBouCoeffs_,
 406                 interfaces_,
 407                 x,
 408                 bPrime_,
 409                 cmpt
 410             );
 411         }
 412
 413         // Forward sweep
 414         forAll (matrix_, rowI)
 415         {
 416             forwardSweepTranspose(matrix_[rowI], x[rowI], bPrime_[rowI]);
 417         }
 418
 419         // Parallel boundary update
 420         {
 421             matrix_.initMatrixInterfaces
 422             (
 423                 mBouCoeffs_,
 424                 interfaces_,
 425                 x,
 426                 bPrime_,
 427                 cmpt
 428             );
 429
 430             matrix_.updateMatrixInterfaces
 431             (
 432                 mBouCoeffs_,
 433                 interfaces_,
 434                 x,
 435                 bPrime_,
 436                 cmpt
 437             );
 438         }
 439
 440         // Reverse sweep
 441         forAllReverse (matrix_, rowI)
 442         {
 443             reverseSweepTranspose(matrix_[rowI], x[rowI], bPrime_[rowI]);
 444         }
 445     }
 446 }
 447
 448
 449 // ************************************************************************* //