src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.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 "LUscalarMatrix.H"
  28 #include "lduMatrix.H"
  29 #include "procLduMatrix.H"
  30 #include "procLduInterface.H"
  31
  32 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  33
  34 Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix)
  35 :
  36     scalarSquareMatrix(matrix),
  37     pivotIndices_(n())
  38 {
  39     LUDecompose(*this, pivotIndices_);
  40 }
  41
  42
  43 Foam::LUscalarMatrix::LUscalarMatrix
  44 (
  45     const lduMatrix& ldum,
  46     const FieldField<Field, scalar>& interfaceCoeffs,
  47     const lduInterfaceFieldPtrsList& interfaces
  48 )
  49 {
  50     if (Pstream::parRun())
  51     {
  52         PtrList<procLduMatrix> lduMatrices(Pstream::nProcs());
  53
  54         label lduMatrixi = 0;
  55
  56         lduMatrices.set
  57         (
  58             lduMatrixi++,
  59             new procLduMatrix
  60             (
  61                 ldum,
  62                 interfaceCoeffs,
  63                 interfaces
  64             )
  65         );
  66
  67         if (Pstream::master())
  68         {
  69             for
  70             (
  71                 int slave=Pstream::firstSlave();
  72                 slave<=Pstream::lastSlave();
  73                 slave++
  74             )
  75             {
  76                 lduMatrices.set
  77                 (
  78                     lduMatrixi++,
  79                     new procLduMatrix(IPstream(Pstream::scheduled, slave)())
  80                 );
  81             }
  82         }
  83         else
  84         {
  85             OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
  86             procLduMatrix cldum
  87             (
  88                 ldum,
  89                 interfaceCoeffs,
  90                 interfaces
  91             );
  92             toMaster<< cldum;
  93
  94         }
  95
  96         if (Pstream::master())
  97         {
  98             label nCells = 0;
  99             forAll(lduMatrices, i)
 100             {
 101                 nCells += lduMatrices[i].size();
 102             }
 103
 104             scalarSquareMatrix m(nCells, nCells, 0.0);
 105             transfer(m);
 106             convert(lduMatrices);
 107         }
 108     }
 109     else
 110     {
 111         label nCells = ldum.lduAddr().size();
 112         scalarSquareMatrix m(nCells, nCells, 0.0);
 113         transfer(m);
 114         convert(ldum, interfaceCoeffs, interfaces);
 115     }
 116
 117     if (Pstream::master())
 118     {
 119         pivotIndices_.setSize(n());
 120         LUDecompose(*this, pivotIndices_);
 121     }
 122 }
 123
 124
 125 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 126
 127 void Foam::LUscalarMatrix::convert
 128 (
 129     const lduMatrix& ldum,
 130     const FieldField<Field, scalar>& interfaceCoeffs,
 131     const lduInterfaceFieldPtrsList& interfaces
 132 )
 133 {
 134     const label* __restrict__ uPtr = ldum.lduAddr().upperAddr().begin();
 135     const label* __restrict__ lPtr = ldum.lduAddr().lowerAddr().begin();
 136
 137     const scalar* __restrict__ diagPtr = ldum.diag().begin();
 138     const scalar* __restrict__ upperPtr = ldum.upper().begin();
 139     const scalar* __restrict__ lowerPtr = ldum.lower().begin();
 140
 141     register const label nCells = ldum.diag().size();
 142     register const label nFaces = ldum.upper().size();
 143
 144     for (register label cell=0; cell<nCells; cell++)
 145     {
 146         operator[](cell)[cell] = diagPtr[cell];
 147     }
 148
 149     for (register label face=0; face<nFaces; face++)
 150     {
 151         label uCell = uPtr[face];
 152         label lCell = lPtr[face];
 153
 154         operator[](uCell)[lCell] = lowerPtr[face];
 155         operator[](lCell)[uCell] = upperPtr[face];
 156     }
 157
 158     forAll(interfaces, inti)
 159     {
 160         if (interfaces.set(inti))
 161         {
 162             const lduInterface& interface = interfaces[inti].interface();
 163
 164             const label* __restrict__ ulPtr = interface.faceCells().begin();
 165             const scalar* __restrict__ upperLowerPtr =
 166                 interfaceCoeffs[inti].begin();
 167
 168             register label inFaces = interface.faceCells().size()/2;
 169
 170             for (register label face=0; face<inFaces; face++)
 171             {
 172                 label uCell = ulPtr[face];
 173                 label lCell = ulPtr[face + inFaces];
 174
 175                 operator[](uCell)[lCell] -= upperLowerPtr[face + inFaces];
 176                 operator[](lCell)[uCell] -= upperLowerPtr[face];
 177             }
 178         }
 179     }
 180
 181     //printDiagonalDominance();
 182 }
 183
 184
 185 void Foam::LUscalarMatrix::convert
 186 (
 187     const PtrList<procLduMatrix>& lduMatrices
 188 )
 189 {
 190     procOffsets_.setSize(lduMatrices.size() + 1);
 191     procOffsets_[0] = 0;
 192
 193     forAll(lduMatrices, ldumi)
 194     {
 195         procOffsets_[ldumi+1] = procOffsets_[ldumi] + lduMatrices[ldumi].size();
 196     }
 197
 198     forAll(lduMatrices, ldumi)
 199     {
 200         const procLduMatrix& lduMatrixi = lduMatrices[ldumi];
 201         label offset = procOffsets_[ldumi];
 202
 203         const label* __restrict__ uPtr = lduMatrixi.upperAddr_.begin();
 204         const label* __restrict__ lPtr = lduMatrixi.lowerAddr_.begin();
 205
 206         const scalar* __restrict__ diagPtr = lduMatrixi.diag_.begin();
 207         const scalar* __restrict__ upperPtr = lduMatrixi.upper_.begin();
 208         const scalar* __restrict__ lowerPtr = lduMatrixi.lower_.begin();
 209
 210         register const label nCells = lduMatrixi.size();
 211         register const label nFaces = lduMatrixi.upper_.size();
 212
 213         for (register label cell=0; cell<nCells; cell++)
 214         {
 215             label globalCell = cell + offset;
 216             operator[](globalCell)[globalCell] = diagPtr[cell];
 217         }
 218
 219         for (register label face=0; face<nFaces; face++)
 220         {
 221             label uCell = uPtr[face] + offset;
 222             label lCell = lPtr[face] + offset;
 223
 224             operator[](uCell)[lCell] = lowerPtr[face];
 225             operator[](lCell)[uCell] = upperPtr[face];
 226         }
 227
 228         const PtrList<procLduInterface>& interfaces =
 229             lduMatrixi.interfaces_;
 230
 231         forAll(interfaces, inti)
 232         {
 233             const procLduInterface& interface = interfaces[inti];
 234
 235             if (interface.myProcNo_ == interface.neighbProcNo_)
 236             {
 237                 const label* __restrict__ ulPtr = interface.faceCells_.begin();
 238
 239                 const scalar* __restrict__ upperLowerPtr =
 240                     interface.coeffs_.begin();
 241
 242                 register label inFaces = interface.faceCells_.size()/2;
 243
 244                 for (register label face=0; face<inFaces; face++)
 245                 {
 246                     label uCell = ulPtr[face] + offset;
 247                     label lCell = ulPtr[face + inFaces] + offset;
 248
 249                     operator[](uCell)[lCell] -= upperLowerPtr[face + inFaces];
 250                     operator[](lCell)[uCell] -= upperLowerPtr[face];
 251                 }
 252             }
 253             else if (interface.myProcNo_ < interface.neighbProcNo_)
 254             {
 255                 const PtrList<procLduInterface>& neiInterfaces =
 256                     lduMatrices[interface.neighbProcNo_].interfaces_;
 257
 258                 label neiInterfacei = -1;
 259
 260                 forAll(neiInterfaces, ninti)
 261                 {
 262                     if
 263                     (
 264                         neiInterfaces[ninti].neighbProcNo_
 265                      == interface.myProcNo_
 266                     )
 267                     {
 268                         neiInterfacei = ninti;
 269                         break;
 270                     }
 271                 }
 272
 273                 if (neiInterfacei == -1)
 274                 {
 275                     FatalErrorIn("LUscalarMatrix::convert") << exit(FatalError);
 276                 }
 277
 278                 const procLduInterface& neiInterface =
 279                     neiInterfaces[neiInterfacei];
 280
 281                 const label* __restrict__ uPtr = interface.faceCells_.begin();
 282                 const label* __restrict__ lPtr = neiInterface.faceCells_.begin();
 283
 284                 const scalar* __restrict__ upperPtr = interface.coeffs_.begin();
 285                 const scalar* __restrict__ lowerPtr = neiInterface.coeffs_.begin();
 286
 287                 register label inFaces = interface.faceCells_.size();
 288                 label neiOffset = procOffsets_[interface.neighbProcNo_];
 289
 290                 for (register label face=0; face<inFaces; face++)
 291                 {
 292                     label uCell = uPtr[face] + offset;
 293                     label lCell = lPtr[face] + neiOffset;
 294
 295                     operator[](uCell)[lCell] -= lowerPtr[face];
 296                     operator[](lCell)[uCell] -= upperPtr[face];
 297                 }
 298             }
 299         }
 300     }
 301
 302     //printDiagonalDominance();
 303 }
 304
 305
 306 void Foam::LUscalarMatrix::printDiagonalDominance() const
 307 {
 308     for (label i=0; i<n(); i++)
 309     {
 310         scalar sum = 0.0;
 311         for (label j=0; j<n(); j++)
 312         {
 313             if (i != j)
 314             {
 315                 sum += operator[](i)[j];
 316             }
 317         }
 318         Info<< mag(sum)/mag(operator[](i)[i]) << endl;
 319     }
 320 }
 321
 322
 323 // ************************************************************************* //