src/OpenFOAM/matrices/lduMatrix/preconditioners/FDICPreconditioner/FDICPreconditioner.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 "FDICPreconditioner.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33     defineTypeNameAndDebug(FDICPreconditioner, 0);
  34
  35     lduMatrix::preconditioner::
  36         addsymMatrixConstructorToTable<FDICPreconditioner>
  37         addFDICPreconditionerSymMatrixConstructorToTable_;
  38 }
  39
  40
  41 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  42
  43 Foam::FDICPreconditioner::FDICPreconditioner
  44 (
  45     const lduMatrix::solver& sol,
  46     const dictionary&
  47 )
  48 :
  49     lduMatrix::preconditioner(sol),
  50     rD_(sol.matrix().diag()),
  51     rDuUpper_(sol.matrix().upper().size()),
  52     rDlUpper_(sol.matrix().upper().size())
  53 {
  54     scalar* __restrict__ rDPtr = rD_.begin();
  55     scalar* __restrict__ rDuUpperPtr = rDuUpper_.begin();
  56     scalar* __restrict__ rDlUpperPtr = rDlUpper_.begin();
  57
  58     const label* const __restrict__ uPtr =
  59         solver_.matrix().lduAddr().upperAddr().begin();
  60     const label* const __restrict__ lPtr =
  61         solver_.matrix().lduAddr().lowerAddr().begin();
  62     const scalar* const __restrict__ upperPtr = solver_.matrix().upper().begin();
  63
  64     register label nCells = rD_.size();
  65     register label nFaces = solver_.matrix().upper().size();
  66
  67     for (register label face=0; face<nFaces; face++)
  68     {
  69         rDPtr[uPtr[face]] -= sqr(upperPtr[face])/rDPtr[lPtr[face]];
  70     }
  71
  72     // Generate reciprocal FDIC
  73     for (register label cell=0; cell<nCells; cell++)
  74     {
  75         rDPtr[cell] = 1.0/rDPtr[cell];
  76     }
  77
  78     for (register label face=0; face<nFaces; face++)
  79     {
  80         rDuUpperPtr[face] = rDPtr[uPtr[face]]*upperPtr[face];
  81         rDlUpperPtr[face] = rDPtr[lPtr[face]]*upperPtr[face];
  82     }
  83 }
  84
  85
  86 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  87
  88 void Foam::FDICPreconditioner::precondition
  89 (
  90     scalarField& wA,
  91     const scalarField& rA,
  92     const direction
  93 ) const
  94 {
  95     scalar* __restrict__ wAPtr = wA.begin();
  96     const scalar* __restrict__ rAPtr = rA.begin();
  97     const scalar* __restrict__ rDPtr = rD_.begin();
  98
  99     const label* const __restrict__ uPtr =
 100         solver_.matrix().lduAddr().upperAddr().begin();
 101     const label* const __restrict__ lPtr =
 102         solver_.matrix().lduAddr().lowerAddr().begin();
 103
 104     const scalar* const __restrict__ rDuUpperPtr = rDuUpper_.begin();
 105     const scalar* const __restrict__ rDlUpperPtr = rDlUpper_.begin();
 106
 107     register label nCells = wA.size();
 108     register label nFaces = solver_.matrix().upper().size();
 109     register label nFacesM1 = nFaces - 1;
 110
 111     for (register label cell=0; cell<nCells; cell++)
 112     {
 113         wAPtr[cell] = rDPtr[cell]*rAPtr[cell];
 114     }
 115
 116     for (register label face=0; face<nFaces; face++)
 117     {
 118         wAPtr[uPtr[face]] -= rDuUpperPtr[face]*wAPtr[lPtr[face]];
 119     }
 120
 121     for (register label face=nFacesM1; face>=0; face--)
 122     {
 123         wAPtr[lPtr[face]] -= rDlUpperPtr[face]*wAPtr[uPtr[face]];
 124     }
 125 }
 126
 127
 128 // ************************************************************************* //