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