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