src/OpenFOAM/matrices/lduMatrix/preconditioners/DICPreconditioner/DICPreconditioner.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 "DICPreconditioner.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33     defineTypeNameAndDebug(DICPreconditioner, 0);
  34
  35     lduMatrix::preconditioner::
  36         addsymMatrixConstructorToTable<DICPreconditioner>
  37         addDICPreconditionerSymMatrixConstructorToTable_;
  38 }
  39
  40
  41 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  42
  43 Foam::DICPreconditioner::DICPreconditioner
  44 (
  45     const lduMatrix::solver& sol,
  46     Istream&
  47 )
  48 :
  49     lduMatrix::preconditioner(sol),
  50     rD_(sol.matrix().diag())
  51 {
  52     calcReciprocalD(rD_, sol.matrix());
  53 }
  54
  55
  56 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  57
  58 void Foam::DICPreconditioner::calcReciprocalD
  59 (
  60     scalarField& rD,
  61     const lduMatrix& matrix
  62 )
  63 {
  64     scalar* __restrict__ rDPtr = rD.begin();
  65
  66     const label* const __restrict__ uPtr = matrix.lduAddr().upperAddr().begin();
  67     const label* const __restrict__ lPtr = matrix.lduAddr().lowerAddr().begin();
  68     const scalar* const __restrict__ upperPtr = matrix.upper().begin();
  69
  70     // Calculate the DIC diagonal
  71     register const label nFaces = matrix.upper().size();
  72     for (register label face=0; face<nFaces; face++)
  73     {
  74         #ifdef ICC_IA64_PREFETCH
  75         __builtin_prefetch (&uPtr[face+96],0,0);
  76         __builtin_prefetch (&lPtr[face+96],0,0);
  77         __builtin_prefetch (&upperPtr[face+96],0,1);
  78         __builtin_prefetch (&rDPtr[lPtr[face+24]],0,1);
  79         __builtin_prefetch (&rDPtr[uPtr[face+24]],1,1);
  80         #endif
  81
  82         rDPtr[uPtr[face]] -= upperPtr[face]*upperPtr[face]/rDPtr[lPtr[face]];
  83     }
  84
  85
  86     // Calculate the reciprocal of the preconditioned diagonal
  87     register const label nCells = rD.size();
  88     #ifdef ICC_IA64_PREFETCH
  89     #pragma ivdep
  90     #endif
  91     for (register label cell=0; cell<nCells; cell++)
  92     {
  93         #ifdef ICC_IA64_PREFETCH
  94         __builtin_prefetch (&rDPtr[cell+96],0,1);
  95         #endif
  96
  97         rDPtr[cell] = 1.0/rDPtr[cell];
  98     }
  99 }
 100
 101
 102 void Foam::DICPreconditioner::precondition
 103 (
 104     scalarField& wA,
 105     const scalarField& rA,
 106     const direction
 107 ) const
 108 {
 109     scalar* __restrict__ wAPtr = wA.begin();
 110     const scalar* __restrict__ rAPtr = rA.begin();
 111     const scalar* __restrict__ rDPtr = rD_.begin();
 112
 113     const label* const __restrict__ uPtr =
 114         solver_.matrix().lduAddr().upperAddr().begin();
 115     const label* const __restrict__ lPtr =
 116         solver_.matrix().lduAddr().lowerAddr().begin();
 117     const scalar* const __restrict__ upperPtr = solver_.matrix().upper().begin();
 118
 119     register label nCells = wA.size();
 120     register label nFaces = solver_.matrix().upper().size();
 121     register label nFacesM1 = nFaces - 1;
 122
 123     #ifdef ICC_IA64_PREFETCH
 124     #pragma ivdep
 125     #endif
 126
 127     for (register label cell=0; cell<nCells; cell++)
 128     {
 129         #ifdef ICC_IA64_PREFETCH
 130         __builtin_prefetch (&wAPtr[cell+96],0,1);
 131         __builtin_prefetch (&rDPtr[cell+96],0,1);
 132         __builtin_prefetch (&rAPtr[cell+96],0,1);
 133         #endif
 134
 135         wAPtr[cell] = rDPtr[cell]*rAPtr[cell];
 136     }
 137
 138     #ifdef ICC_IA64_PREFETCH
 139     #pragma noprefetch uPtr,lPtr,upperPtr,rDPtr,wAPtr
 140     #pragma nounroll
 141     #endif
 142
 143     for (register label face=0; face<nFaces; face++)
 144     {
 145         #ifdef ICC_IA64_PREFETCH
 146         __builtin_prefetch (&uPtr[face+96],0,0);
 147         __builtin_prefetch (&lPtr[face+96],0,0);
 148         __builtin_prefetch (&upperPtr[face+96],0,0);
 149         __builtin_prefetch (&rDPtr[uPtr[face+32]],0,1);
 150         __builtin_prefetch (&wAPtr[uPtr[face+32]],0,1);
 151         __builtin_prefetch (&wAPtr[lPtr[face+32]],0,1);
 152         #endif
 153
 154         wAPtr[uPtr[face]] -= rDPtr[uPtr[face]]*upperPtr[face]*wAPtr[lPtr[face]];
 155     }
 156
 157     #ifdef ICC_IA64_PREFETCH
 158     #pragma noprefetch uPtr,lPtr,rDPtr,wAPtr
 159     #pragma nounroll
 160     #endif
 161
 162     for (register label face=nFacesM1; face>=0; face--)
 163     {
 164         #ifdef ICC_IA64_PREFETCH
 165         __builtin_prefetch (&uPtr[face-95],0,0);
 166         __builtin_prefetch (&lPtr[face-95],0,0);
 167         __builtin_prefetch (&rDPtr[lPtr[face-16]],0,1);
 168         __builtin_prefetch (&wAPtr[lPtr[face-16]],0,1);
 169         __builtin_prefetch (&wAPtr[uPtr[face-16]],0,1);
 170         __builtin_prefetch (&rDPtr[lPtr[face-24]],0,1);
 171         __builtin_prefetch (&wAPtr[lPtr[face-24]],0,1);
 172         __builtin_prefetch (&wAPtr[uPtr[face-24]],0,1);
 173         __builtin_prefetch (&rDPtr[lPtr[face-32]],0,1);
 174         __builtin_prefetch (&wAPtr[lPtr[face-32]],0,1);
 175         __builtin_prefetch (&wAPtr[uPtr[face-32]],0,1);
 176         #endif
 177
 178         wAPtr[lPtr[face]] -= rDPtr[lPtr[face]]*upperPtr[face]*wAPtr[uPtr[face]];
 179     }
 180 }
 181
 182
 183 // ************************************************************************* //