src/lduSolvers/lduPrecon/symGaussSeidelPrecon/symGaussSeidelPrecon.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     symGaussSeidelPrecon
  26
  27 Description
  28     Symmetric Gauss-Seidel preconditioning
  29
  30 Author
  31     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "symGaussSeidelPrecon.H"
  36 #include "addToRunTimeSelectionTable.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42     defineTypeNameAndDebug(symGaussSeidelPrecon, 0);
  43
  44     lduPreconditioner::
  45         addsymMatrixConstructorToTable<symGaussSeidelPrecon>
  46         addsymGaussSeidelPreconditionerSymMatrixConstructorToTable_;
  47
  48     lduPreconditioner::
  49         addasymMatrixConstructorToTable<symGaussSeidelPrecon>
  50         addsymGaussSeidelPreconditionerAsymMatrixConstructorToTable_;
  51 }
  52
  53
  54 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  55
  56 Foam::symGaussSeidelPrecon::symGaussSeidelPrecon
  57 (
  58     const lduMatrix& matrix,
  59     const FieldField<Field, scalar>& coupleBouCoeffs,
  60     const FieldField<Field, scalar>& coupleIntCoeffs,
  61     const lduInterfaceFieldPtrsList& interfaces
  62 )
  63 :
  64     lduPreconditioner
  65     (
  66         matrix,
  67         coupleBouCoeffs,
  68         coupleIntCoeffs,
  69         interfaces
  70     ),
  71     bPrime_(matrix.lduAddr().size())
  72 {}
  73
  74
  75 Foam::symGaussSeidelPrecon::symGaussSeidelPrecon
  76 (
  77     const lduMatrix& matrix,
  78     const FieldField<Field, scalar>& coupleBouCoeffs,
  79     const FieldField<Field, scalar>& coupleIntCoeffs,
  80     const lduInterfaceFieldPtrsList& interfaces,
  81     const dictionary& dict
  82 )
  83 :
  84     lduPreconditioner
  85     (
  86         matrix,
  87         coupleBouCoeffs,
  88         coupleIntCoeffs,
  89         interfaces
  90     ),
  91     bPrime_(matrix.lduAddr().size())
  92 {}
  93
  94
  95 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  96
  97 void Foam::symGaussSeidelPrecon::precondition
  98 (
  99     scalarField& x,
 100     const scalarField& b,
 101     const direction cmpt
 102 ) const
 103 {
 104     // Execute preconditioning
 105     if (matrix_.diagonal())
 106     {
 107         x = b/matrix_.diag();
 108     }
 109     else if (matrix_.symmetric() || matrix_.asymmetric())
 110     {
 111         scalar* __restrict__ xPtr = x.begin();
 112
 113         const scalar* __restrict__ diagPtr = matrix_.diag().begin();
 114
 115         scalar* __restrict__ bPrimePtr = bPrime_.begin();
 116
 117         const label* const __restrict__ uPtr =
 118             matrix_.lduAddr().upperAddr().begin();
 119
 120         const label* const __restrict__ ownStartPtr =
 121             matrix_.lduAddr().ownerStartAddr().begin();
 122
 123         const scalar* const __restrict__ lowerPtr = matrix_.lower().begin();
 124         const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
 125
 126         const label nRows = x.size();
 127
 128         label fStart, fEnd;
 129
 130         bPrime_ = b;
 131
 132         // Coupled boundary update
 133         {
 134             matrix_.initMatrixInterfaces
 135             (
 136                 coupleBouCoeffs_,
 137                 interfaces_,
 138                 x,
 139                 bPrime_,
 140                 cmpt,
 141                 true             // switch to lhs of system
 142             );
 143
 144             matrix_.updateMatrixInterfaces
 145             (
 146                 coupleBouCoeffs_,
 147                 interfaces_,
 148                 x,
 149                 bPrime_,
 150                 cmpt,
 151                 true             // switch to lhs of system
 152             );
 153         }
 154
 155         // Forward sweep
 156         for (register label rowI = 0; rowI < nRows; rowI++)
 157         {
 158             // lRow is equal to rowI
 159             scalar& curX = xPtr[rowI];
 160
 161             // Grab the accumulated neighbour side
 162             curX = bPrimePtr[rowI];
 163
 164             // Start and end of this row
 165             fStart = ownStartPtr[rowI];
 166             fEnd = ownStartPtr[rowI + 1];
 167
 168             // Accumulate the owner product side
 169             for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 170             {
 171                 curX -= upperPtr[curCoeff]*xPtr[uPtr[curCoeff]];
 172             }
 173
 174             // Finish current x
 175             curX /= diagPtr[rowI];
 176
 177             // Distribute the neighbour side using current x
 178             for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 179             {
 180                 bPrimePtr[uPtr[curCoeff]] -= lowerPtr[curCoeff]*curX;
 181             }
 182         }
 183
 184         // Reverse sweep
 185         for (register label rowI = nRows - 1; rowI >= 0; rowI--)
 186         {
 187             // lRow is equal to rowI
 188             scalar& curX = xPtr[rowI];
 189
 190             // Grab the accumulated neighbour side
 191             curX = bPrimePtr[rowI];
 192
 193             // Start and end of this row
 194             fStart = ownStartPtr[rowI];
 195             fEnd = ownStartPtr[rowI + 1];
 196
 197             // Accumulate the owner product side
 198             for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 199             {
 200                 curX -= upperPtr[curCoeff]*xPtr[uPtr[curCoeff]];
 201             }
 202
 203             // Finish current x
 204             curX /= diagPtr[rowI];
 205
 206             // No need to update bPrime on reverse sweep. VV, 20/May/2015.
 207         }
 208     }
 209 }
 210
 211
 212 void Foam::symGaussSeidelPrecon::preconditionT
 213 (
 214     scalarField& x,
 215     const scalarField& b,
 216     const direction cmpt
 217 ) const
 218 {
 219     // Execute preconditioning
 220     if (matrix_.diagonal())
 221     {
 222         x = b/matrix_.diag();
 223     }
 224     else if (matrix_.symmetric() || matrix_.asymmetric())
 225     {
 226         scalar* __restrict__ xPtr = x.begin();
 227
 228         const scalar* __restrict__ diagPtr = matrix_.diag().begin();
 229
 230         scalar* __restrict__ bPrimePtr = bPrime_.begin();
 231
 232         const label* const __restrict__ uPtr =
 233             matrix_.lduAddr().upperAddr().begin();
 234
 235         const label* const __restrict__ ownStartPtr =
 236             matrix_.lduAddr().ownerStartAddr().begin();
 237
 238         const scalar* const __restrict__ lowerPtr = matrix_.lower().begin();
 239         const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
 240
 241         const label nRows = x.size();
 242
 243         label fStart, fEnd;
 244
 245         bPrime_ = b;
 246
 247         // Coupled boundary update
 248         {
 249             matrix_.initMatrixInterfaces
 250             (
 251                 coupleBouCoeffs_,
 252                 interfaces_,
 253                 x,
 254                 bPrime_,
 255                 cmpt,
 256                 true             // switch to lhs of system
 257             );
 258
 259             matrix_.updateMatrixInterfaces
 260             (
 261                 coupleBouCoeffs_,
 262                 interfaces_,
 263                 x,
 264                 bPrime_,
 265                 cmpt,
 266                 true             // switch to lhs of system
 267             );
 268         }
 269
 270         // Forward sweep
 271         for (register label rowI = 0; rowI < nRows; rowI++)
 272         {
 273             // lRow is equal to rowI
 274             scalar& curX = xPtr[rowI];
 275
 276             // Grab the accumulated neighbour side
 277             curX = bPrimePtr[rowI];
 278
 279             // Start and end of this row
 280             fStart = ownStartPtr[rowI];
 281             fEnd = ownStartPtr[rowI + 1];
 282
 283             // Accumulate the owner product side
 284             for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 285             {
 286                 // Transpose multiplication.  HJ, 10/Jul/2007
 287                 curX -= lowerPtr[curCoeff]*xPtr[uPtr[curCoeff]];
 288             }
 289
 290             // Finish current x
 291             curX /= diagPtr[rowI];
 292
 293             // Distribute the neighbour side using current x
 294             for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 295             {
 296                 // Transpose multiplication.  HJ, 10/Jul/2007
 297                 bPrimePtr[uPtr[curCoeff]] -= upperPtr[curCoeff]*curX;
 298             }
 299         }
 300
 301         // Reverse sweep
 302         for (register label rowI = nRows - 1; rowI >= 0; rowI--)
 303         {
 304             // lRow is equal to rowI
 305             scalar& curX = xPtr[rowI];
 306
 307             // Grab the accumulated neighbour side
 308             curX = bPrimePtr[rowI];
 309
 310             // Start and end of this row
 311             fStart = ownStartPtr[rowI];
 312             fEnd = ownStartPtr[rowI + 1];
 313
 314             // Accumulate the owner product side
 315             for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
 316             {
 317                 // Transpose multiplication.  HJ, 10/Jul/2007
 318                 curX -= lowerPtr[curCoeff]*xPtr[uPtr[curCoeff]];
 319             }
 320
 321             // Finish current x
 322             curX /= diagPtr[rowI];
 323
 324             // No need to update bPrime on reverse sweep. VV, 20/May/2015.
 325         }
 326     }
 327 }
 328
 329
 330 // ************************************************************************* //