src/coupledMatrix/coupledLduPrecon/GaussSeidelPrecon/coupledGaussSeidelPrecon.C

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