src/OpenFOAM/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, 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     lduMatrix::smoother::addsymMatrixConstructorToTable<GaussSeidelSmoother>
  35         addGaussSeidelSmootherSymMatrixConstructorToTable_;
  36
  37     lduMatrix::smoother::addasymMatrixConstructorToTable<GaussSeidelSmoother>
  38         addGaussSeidelSmootherAsymMatrixConstructorToTable_;
  39 }
  40
  41
  42 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  43
  44 Foam::GaussSeidelSmoother::GaussSeidelSmoother
  45 (
  46     const word& fieldName,
  47     const lduMatrix& matrix,
  48     const FieldField<Field, scalar>& interfaceBouCoeffs,
  49     const FieldField<Field, scalar>& interfaceIntCoeffs,
  50     const lduInterfaceFieldPtrsList& interfaces
  51 )
  52 :
  53     lduMatrix::smoother
  54     (
  55         fieldName,
  56         matrix,
  57         interfaceBouCoeffs,
  58         interfaceIntCoeffs,
  59         interfaces
  60     )
  61 {}
  62
  63
  64 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  65
  66 void Foam::GaussSeidelSmoother::smooth
  67 (
  68     const word& fieldName_,
  69     scalarField& psi,
  70     const lduMatrix& matrix_,
  71     const scalarField& source,
  72     const FieldField<Field, scalar>& interfaceBouCoeffs_,
  73     const lduInterfaceFieldPtrsList& interfaces_,
  74     const direction cmpt,
  75     const label nSweeps
  76 )
  77 {
  78     register scalar* __restrict__ psiPtr = psi.begin();
  79
  80     register const label nCells = psi.size();
  81
  82     scalarField bPrime(nCells);
  83     register scalar* __restrict__ bPrimePtr = bPrime.begin();
  84
  85     register const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
  86     register const scalar* const __restrict__ upperPtr =
  87         matrix_.upper().begin();
  88     register const scalar* const __restrict__ lowerPtr =
  89         matrix_.lower().begin();
  90
  91     register const label* const __restrict__ uPtr =
  92         matrix_.lduAddr().upperAddr().begin();
  93
  94     register const label* const __restrict__ ownStartPtr =
  95         matrix_.lduAddr().ownerStartAddr().begin();
  96
  97
  98     // Parallel boundary initialisation.  The parallel boundary is treated
  99     // as an effective jacobi interface in the boundary.
 100     // Note: there is a change of sign in the coupled
 101     // interface update.  The reason for this is that the
 102     // internal coefficients are all located at the l.h.s. of
 103     // the matrix whereas the "implicit" coefficients on the
 104     // coupled boundaries are all created as if the
 105     // coefficient contribution is of a source-kind (i.e. they
 106     // have a sign as if they are on the r.h.s. of the matrix.
 107     // To compensate for this, it is necessary to turn the
 108     // sign of the contribution.
 109
 110     FieldField<Field, scalar> mBouCoeffs(interfaceBouCoeffs_.size());
 111
 112     forAll(mBouCoeffs, patchi)
 113     {
 114         if (interfaces_.set(patchi))
 115         {
 116             mBouCoeffs.set(patchi, -interfaceBouCoeffs_[patchi]);
 117         }
 118     }
 119
 120     for (label sweep=0; sweep<nSweeps; sweep++)
 121     {
 122         bPrime = source;
 123
 124         matrix_.initMatrixInterfaces
 125         (
 126             mBouCoeffs,
 127             interfaces_,
 128             psi,
 129             bPrime,
 130             cmpt
 131         );
 132
 133         matrix_.updateMatrixInterfaces
 134         (
 135             mBouCoeffs,
 136             interfaces_,
 137             psi,
 138             bPrime,
 139             cmpt
 140         );
 141
 142         register scalar curPsi;
 143         register label fStart;
 144         register label fEnd = ownStartPtr[0];
 145
 146         for (register label cellI=0; cellI<nCells; cellI++)
 147         {
 148             // Start and end of this row
 149             fStart = fEnd;
 150             fEnd = ownStartPtr[cellI + 1];
 151
 152             // Get the accumulated neighbour side
 153             curPsi = bPrimePtr[cellI];
 154
 155             // Accumulate the owner product side
 156             for (register label curFace=fStart; curFace<fEnd; curFace++)
 157             {
 158                 curPsi -= upperPtr[curFace]*psiPtr[uPtr[curFace]];
 159             }
 160
 161             // Finish current psi
 162             curPsi /= diagPtr[cellI];
 163
 164             // Distribute the neighbour side using current psi
 165             for (register label curFace=fStart; curFace<fEnd; curFace++)
 166             {
 167                 bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curPsi;
 168             }
 169
 170             psiPtr[cellI] = curPsi;
 171         }
 172     }
 173 }
 174
 175
 176 void Foam::GaussSeidelSmoother::smooth
 177 (
 178     scalarField& psi,
 179     const scalarField& source,
 180     const direction cmpt,
 181     const label nSweeps
 182 ) const
 183 {
 184     smooth
 185     (
 186         fieldName_,
 187         psi,
 188         matrix_,
 189         source,
 190         interfaceBouCoeffs_,
 191         interfaces_,
 192         cmpt,
 193         nSweeps
 194     );
 195 }
 196
 197
 198 // ************************************************************************* //