src/foam/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "GaussSeidelSmoother.H"
  27
  28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  29
  30 namespace Foam
  31 {
  32     defineTypeNameAndDebug(GaussSeidelSmoother, 0);
  33
  34     lduSmoother::addsymMatrixConstructorToTable<GaussSeidelSmoother>
  35         addGaussSeidelSmootherSymMatrixConstructorToTable_;
  36
  37     lduSmoother::addasymMatrixConstructorToTable<GaussSeidelSmoother>
  38         addGaussSeidelSmootherAsymMatrixConstructorToTable_;
  39 }
  40
  41
  42 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  43
  44 Foam::GaussSeidelSmoother::GaussSeidelSmoother
  45 (
  46     const lduMatrix& matrix,
  47     const FieldField<Field, scalar>& coupleBouCoeffs,
  48     const FieldField<Field, scalar>& coupleIntCoeffs,
  49     const lduInterfaceFieldPtrsList& interfaces
  50 )
  51 :
  52     lduSmoother
  53     (
  54         matrix,
  55         coupleBouCoeffs,
  56         coupleIntCoeffs,
  57         interfaces
  58     )
  59 {}
  60
  61
  62 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  63
  64 void Foam::GaussSeidelSmoother::smooth
  65 (
  66     scalarField& x,
  67     const lduMatrix& matrix_,
  68     const scalarField& b,
  69     const FieldField<Field, scalar>& coupleBouCoeffs_,
  70     const lduInterfaceFieldPtrsList& interfaces_,
  71     const direction cmpt,
  72     const label nSweeps
  73 )
  74 {
  75     register scalar* __restrict__ xPtr = x.begin();
  76
  77     register const label nCells = x.size();
  78
  79     scalarField bPrime(nCells);
  80     register scalar* __restrict__ bPrimePtr = bPrime.begin();
  81
  82     register const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
  83     register const scalar* const __restrict__ upperPtr =
  84         matrix_.upper().begin();
  85     register const scalar* const __restrict__ lowerPtr =
  86         matrix_.lower().begin();
  87
  88     register const label* const __restrict__ uPtr =
  89         matrix_.lduAddr().upperAddr().begin();
  90
  91     register const label* const __restrict__ ownStartPtr =
  92         matrix_.lduAddr().ownerStartAddr().begin();
  93
  94
  95     // Coupled boundary initialisation.  The coupled boundary is treated
  96     // as an effective jacobi interface in the boundary.
  97     // Note: there is a change of sign in the coupled
  98     // interface update.  The reason for this is that the
  99     // internal coefficients are all located at the l.h.s. of
 100     // the matrix whereas the "implicit" coefficients on the
 101     // coupled boundaries are all created as if the
 102     // coefficient contribution is of a b-kind (i.e. they
 103     // have a sign as if they are on the r.h.s. of the matrix.
 104     // To compensate for this, it is necessary to turn the
 105     // sign of the contribution.
 106
 107     // Handled by LHS switch on initMatrixInterfaces and updateMatrixInterfaces
 108     // HJ, 22/May/2013
 109
 110     for (label sweep = 0; sweep < nSweeps; sweep++)
 111     {
 112         bPrime = b;
 113
 114         // Update from lhs
 115         matrix_.initMatrixInterfaces
 116         (
 117             coupleBouCoeffs_,
 118             interfaces_,
 119             x,
 120             bPrime,
 121             cmpt,
 122             true         // switch to lhs
 123         );
 124
 125         // Update from lhs
 126         matrix_.updateMatrixInterfaces
 127         (
 128             coupleBouCoeffs_,
 129             interfaces_,
 130             x,
 131             bPrime,
 132             cmpt,
 133             true         // switch to lhs
 134         );
 135
 136         register scalar curX;
 137         register label fStart;
 138         register label fEnd = ownStartPtr[0];
 139
 140         for (register label cellI = 0; cellI < nCells; cellI++)
 141         {
 142             // Start and end of this row
 143             fStart = fEnd;
 144             fEnd = ownStartPtr[cellI + 1];
 145
 146             // Get the accumulated neighbour side
 147             curX = bPrimePtr[cellI];
 148
 149             // Accumulate the owner product side
 150             for (register label curFace = fStart; curFace < fEnd; curFace++)
 151             {
 152                 curX -= upperPtr[curFace]*xPtr[uPtr[curFace]];
 153             }
 154
 155             // Finish current x
 156             curX /= diagPtr[cellI];
 157
 158             // Distribute the neighbour side using current x
 159             for (register label curFace = fStart; curFace < fEnd; curFace++)
 160             {
 161                 bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curX;
 162             }
 163
 164             xPtr[cellI] = curX;
 165         }
 166     }
 167 }
 168
 169
 170 void Foam::GaussSeidelSmoother::smooth
 171 (
 172     scalarField& x,
 173     const scalarField& b,
 174     const direction cmpt,
 175     const label nSweeps
 176 ) const
 177 {
 178     smooth
 179     (
 180         x,
 181         matrix_,
 182         b,
 183         coupleBouCoeffs_,
 184         interfaces_,
 185         cmpt,
 186         nSweeps
 187     );
 188 }
 189
 190
 191 // ************************************************************************* //