src/foam/matrices/lduMatrix/preconditioners/DILUPreconditioner/DILUPreconditioner.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 "DILUPreconditioner.H"
  27
  28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  29
  30 namespace Foam
  31 {
  32     defineTypeNameAndDebug(DILUPreconditioner, 0);
  33
  34     lduPreconditioner::
  35         addsymMatrixConstructorToTable<DILUPreconditioner>
  36         addDILUPreconditionerSymMatrixConstructorToTable_;
  37
  38     lduPreconditioner::
  39         addasymMatrixConstructorToTable<DILUPreconditioner>
  40         addDILUPreconditionerAsymMatrixConstructorToTable_;
  41 }
  42
  43
  44 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  45
  46 Foam::DILUPreconditioner::DILUPreconditioner
  47 (
  48     const lduMatrix& matrix,
  49     const FieldField<Field, scalar>& coupleBouCoeffs,
  50     const FieldField<Field, scalar>& coupleIntCoeffs,
  51     const lduInterfaceFieldPtrsList& interfaces,
  52     const dictionary&
  53 )
  54 :
  55     lduPreconditioner
  56     (
  57         matrix,
  58         coupleBouCoeffs,
  59         coupleIntCoeffs,
  60         interfaces
  61     ),
  62     rD_(matrix.diag())
  63 {
  64     calcReciprocalD(rD_, matrix);
  65 }
  66
  67
  68 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  69
  70 void Foam::DILUPreconditioner::calcReciprocalD
  71 (
  72     scalarField& rD,
  73     const lduMatrix& matrix
  74 )
  75 {
  76     scalar* __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     register label nFaces = matrix.upper().size();
  87     for (register label face=0; face<nFaces; face++)
  88     {
  89         rDPtr[uPtr[face]] -= upperPtr[face]*lowerPtr[face]/rDPtr[lPtr[face]];
  90     }
  91
  92
  93     // Calculate the reciprocal of the preconditioned diagonal
  94     register label nCells = rD.size();
  95
  96     for (register label cell=0; cell<nCells; cell++)
  97     {
  98         rDPtr[cell] = 1.0/rDPtr[cell];
  99     }
 100 }
 101
 102
 103 void Foam::DILUPreconditioner::precondition
 104 (
 105     scalarField& wA,
 106     const scalarField& rA,
 107     const direction
 108 ) const
 109 {
 110     scalar* __restrict__ wAPtr = wA.begin();
 111     const scalar* __restrict__ rAPtr = rA.begin();
 112     const scalar* __restrict__ rDPtr = rD_.begin();
 113
 114     const label* const __restrict__ uPtr =
 115         matrix_.lduAddr().upperAddr().begin();
 116     const label* const __restrict__ lPtr =
 117         matrix_.lduAddr().lowerAddr().begin();
 118     const label* const __restrict__ losortPtr =
 119         matrix_.lduAddr().losortAddr().begin();
 120
 121     const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
 122     const scalar* const __restrict__ lowerPtr = matrix_.lower().begin();
 123
 124     register label nCells = wA.size();
 125     register label nFaces = matrix_.upper().size();
 126     register label nFacesM1 = nFaces - 1;
 127
 128     for (register label cell=0; cell<nCells; cell++)
 129     {
 130         wAPtr[cell] = rDPtr[cell]*rAPtr[cell];
 131     }
 132
 133
 134     register label sface;
 135
 136     for (register label face=0; face<nFaces; face++)
 137     {
 138         sface = losortPtr[face];
 139         wAPtr[uPtr[sface]] -=
 140             rDPtr[uPtr[sface]]*lowerPtr[sface]*wAPtr[lPtr[sface]];
 141     }
 142
 143     for (register label face=nFacesM1; face>=0; face--)
 144     {
 145         wAPtr[lPtr[face]] -=
 146             rDPtr[lPtr[face]]*upperPtr[face]*wAPtr[uPtr[face]];
 147     }
 148 }
 149
 150
 151 void Foam::DILUPreconditioner::preconditionT
 152 (
 153     scalarField& wT,
 154     const scalarField& rT,
 155     const direction
 156 ) const
 157 {
 158     scalar* __restrict__ wTPtr = wT.begin();
 159     const scalar* __restrict__ rTPtr = rT.begin();
 160     const scalar* __restrict__ rDPtr = rD_.begin();
 161
 162     const label* const __restrict__ uPtr =
 163         matrix_.lduAddr().upperAddr().begin();
 164     const label* const __restrict__ lPtr =
 165         matrix_.lduAddr().lowerAddr().begin();
 166     const label* const __restrict__ losortPtr =
 167         matrix_.lduAddr().losortAddr().begin();
 168
 169     const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
 170     const scalar* const __restrict__ lowerPtr = matrix_.lower().begin();
 171
 172     register label nCells = wT.size();
 173     register label nFaces = matrix_.upper().size();
 174     register label nFacesM1 = nFaces - 1;
 175
 176     for (register label cell=0; cell<nCells; cell++)
 177     {
 178         wTPtr[cell] = rDPtr[cell]*rTPtr[cell];
 179     }
 180
 181     for (register label face=0; face<nFaces; face++)
 182     {
 183         wTPtr[uPtr[face]] -=
 184             rDPtr[uPtr[face]]*upperPtr[face]*wTPtr[lPtr[face]];
 185     }
 186
 187
 188     register label sface;
 189
 190     for (register label face=nFacesM1; face>=0; face--)
 191     {
 192         sface = losortPtr[face];
 193         wTPtr[lPtr[sface]] -=
 194             rDPtr[lPtr[sface]]*lowerPtr[sface]*wTPtr[uPtr[sface]];
 195     }
 196 }
 197
 198
 199 // ************************************************************************* //