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