src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 "DILUSmoother.H"
  28 #include "DILUPreconditioner.H"
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34     defineTypeNameAndDebug(DILUSmoother, 0);
  35
  36     lduMatrix::smoother::addasymMatrixConstructorToTable<DILUSmoother>
  37         addDILUSmootherAsymMatrixConstructorToTable_;
  38 }
  39
  40
  41 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  42
  43 Foam::DILUSmoother::DILUSmoother
  44 (
  45     const word& fieldName,
  46     const lduMatrix& matrix,
  47     const FieldField<Field, scalar>& interfaceBouCoeffs,
  48     const FieldField<Field, scalar>& interfaceIntCoeffs,
  49     const lduInterfaceFieldPtrsList& interfaces
  50 )
  51 :
  52     lduMatrix::smoother
  53     (
  54         fieldName,
  55         matrix,
  56         interfaceBouCoeffs,
  57         interfaceIntCoeffs,
  58         interfaces
  59     ),
  60     rD_(matrix_.diag())
  61 {
  62     DILUPreconditioner::calcReciprocalD(rD_, matrix_);
  63 }
  64
  65
  66 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  67
  68 void Foam::DILUSmoother::smooth
  69 (
  70     scalarField& psi,
  71     const scalarField& source,
  72     const direction cmpt,
  73     const label nSweeps
  74 ) const
  75 {
  76     const scalar* const __restrict__ rDPtr = rD_.begin();
  77
  78     const label* const __restrict__ uPtr =
  79         matrix_.lduAddr().upperAddr().begin();
  80     const label* const __restrict__ lPtr =
  81         matrix_.lduAddr().lowerAddr().begin();
  82
  83     const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
  84     const scalar* const __restrict__ lowerPtr = matrix_.lower().begin();
  85
  86     // Temporary storage for the residual
  87     scalarField rA(rD_.size());
  88     scalar* __restrict__ rAPtr = rA.begin();
  89
  90     for (label sweep=0; sweep<nSweeps; sweep++)
  91     {
  92         matrix_.residual
  93         (
  94             rA,
  95             psi,
  96             source,
  97             interfaceBouCoeffs_,
  98             interfaces_,
  99             cmpt
 100         );
 101
 102         rA *= rD_;
 103
 104         register label nFaces = matrix_.upper().size();
 105         for (register label face=0; face<nFaces; face++)
 106         {
 107             register label u = uPtr[face];
 108             rAPtr[u] -= rDPtr[u]*lowerPtr[face]*rAPtr[lPtr[face]];
 109         }
 110
 111         register label nFacesM1 = nFaces - 1;
 112         for (register label face=nFacesM1; face>=0; face--)
 113         {
 114             register label l = lPtr[face];
 115             rAPtr[l] -= rDPtr[l]*upperPtr[face]*rAPtr[uPtr[face]];
 116         }
 117
 118         psi += rA;
 119     }
 120 }
 121
 122
 123 // ************************************************************************* //