src/coupledMatrix/coupledLduMatrix/coupledLduMatrix.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 Class
  25     coupledLduMatrix
  26
  27 Description
  28     Collection of lduMatrices solved together as a block system
  29
  30 Author
  31     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "coupledLduMatrix.H"
  36 #include "processorLduInterfaceField.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42     defineTypeNameAndDebug(coupledLduMatrix, 1);
  43 }
  44
  45
  46 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  47
  48 // Construct given size
  49 Foam::coupledLduMatrix::coupledLduMatrix(const label size)
  50 :
  51     PtrList<lduMatrix>(size)
  52 {}
  53
  54
  55 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  56
  57 Foam::coupledLduMatrix::~coupledLduMatrix()
  58 {}
  59
  60
  61 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  62
  63 bool Foam::coupledLduMatrix::diagonal() const
  64 {
  65     const PtrList<lduMatrix>& matrices = *this;
  66
  67     bool diag = true;
  68
  69     forAll (matrices, matrixI)
  70     {
  71         diag = diag && matrices[matrixI].diagonal();
  72     }
  73
  74     return diag;
  75 }
  76
  77
  78 bool Foam::coupledLduMatrix::symmetric() const
  79 {
  80     const PtrList<lduMatrix>& matrices = *this;
  81
  82     bool sym = true;
  83
  84     forAll (matrices, matrixI)
  85     {
  86         sym =
  87             (sym && matrices[matrixI].diagonal())
  88          || (sym && matrices[matrixI].symmetric());
  89     }
  90
  91     return sym;
  92 }
  93
  94
  95 bool Foam::coupledLduMatrix::asymmetric() const
  96 {
  97     const PtrList<lduMatrix>& matrices = *this;
  98
  99     bool asym = false;
 100
 101     forAll (matrices, matrixI)
 102     {
 103         asym = (asym || matrices[matrixI].asymmetric());
 104     }
 105
 106     return asym;
 107 }
 108
 109
 110 void Foam::coupledLduMatrix::Amul
 111 (
 112     FieldField<Field, scalar>& result,
 113     const FieldField<Field, scalar>& x,
 114     const PtrList<FieldField<Field, scalar> >& bouCoeffs,
 115     const lduInterfaceFieldPtrsListList& interfaces,
 116     const direction cmpt
 117 ) const
 118 {
 119     const PtrList<lduMatrix>& matrices = *this;
 120
 121     // Reset product to zero
 122     result = 0;
 123
 124     // Initialise the update of coupled interfaces
 125     initMatrixInterfaces
 126     (
 127         bouCoeffs,
 128         interfaces,
 129         x,
 130         result,
 131         cmpt
 132     );
 133
 134     forAll (matrices, rowI)
 135     {
 136         matrices[rowI].AmulCore(result[rowI], x[rowI]);
 137     }
 138
 139     // Update couple interfaces
 140     updateMatrixInterfaces
 141     (
 142         bouCoeffs,
 143         interfaces,
 144         x,
 145         result,
 146         cmpt
 147     );
 148 }
 149
 150
 151 void Foam::coupledLduMatrix::Tmul
 152 (
 153     FieldField<Field, scalar>& result,
 154     const FieldField<Field, scalar>& x,
 155     const PtrList<FieldField<Field, scalar> >& intCoeffs,
 156     const lduInterfaceFieldPtrsListList& interfaces,
 157     const direction cmpt
 158 ) const
 159 {
 160     const PtrList<lduMatrix>& matrices = *this;
 161
 162     // Reset product to zero
 163     result = 0;
 164
 165     // Initialise the update of coupled interfaces
 166     initMatrixInterfaces
 167     (
 168         intCoeffs,
 169         interfaces,
 170         x,
 171         result,
 172         cmpt
 173     );
 174
 175     forAll (matrices, rowI)
 176     {
 177         matrices[rowI].TmulCore(result[rowI], x[rowI]);
 178     }
 179
 180     // Update couple interfaces
 181     updateMatrixInterfaces
 182     (
 183         intCoeffs,
 184         interfaces,
 185         x,
 186         result,
 187         cmpt
 188     );
 189 }
 190
 191
 192 void Foam::coupledLduMatrix::initMatrixInterfaces
 193 (
 194     const PtrList<FieldField<Field, scalar> >& coupleCoeffs,
 195     const lduInterfaceFieldPtrsListList& interfaces,
 196     const FieldField<Field, scalar>& x,
 197     FieldField<Field, scalar>& result,
 198     const direction cmpt
 199 ) const
 200 {
 201     const PtrList<lduMatrix>& matrices = *this;
 202
 203     // Note.  The comms design requires all non-processor interfaces
 204     // to be updated first, followed by the update of processor
 205     // interfaces.  The reason is that non-processor coupled
 206     // interfaces require a complex comms pattern involving more than
 207     // pairwise communications.
 208     // Under normal circumstances this is achieved naturall, since
 209     // processor interfaces come last on the list and other coupled
 210     // interfaces execute complex comms at init() level.
 211     // For coupled matrices, the update loop needs to be split over
 212     // all matrices by hand
 213     // Bug fix: Zeljko Tukovic, 7/Apr/2015
 214
 215     // Init update all non-processor coupled interfaces
 216     forAll (matrices, rowI)
 217     {
 218         if
 219         (
 220             Pstream::defaultCommsType() == Pstream::blocking
 221          || Pstream::defaultCommsType() == Pstream::nonBlocking
 222         )
 223         {
 224             forAll (interfaces[rowI], interfaceI)
 225             {
 226                 if (interfaces[rowI].set(interfaceI))
 227                 {
 228                     if
 229                     (
 230                        !isA<processorLduInterfaceField>
 231                         (
 232                             interfaces[rowI][interfaceI]
 233                         )
 234                     )
 235                     {
 236                         interfaces[rowI][interfaceI].initInterfaceMatrixUpdate
 237                         (
 238                             x[rowI],
 239                             result[rowI],
 240                             matrices[rowI],
 241                             coupleCoeffs[rowI][interfaceI],
 242                             cmpt,
 243                             static_cast<const Pstream::commsTypes>
 244                             (
 245                                 Pstream::defaultCommsType()
 246                             ),
 247                             false
 248                         );
 249                     }
 250                 }
 251             }
 252         }
 253         else
 254         {
 255             matrices[rowI].initMatrixInterfaces
 256             (
 257                 coupleCoeffs[rowI],
 258                 interfaces[rowI],
 259                 x[rowI],
 260                 result[rowI],
 261                 cmpt
 262             );
 263         }
 264     }
 265
 266     // Init update for all processor interfaces
 267     forAll (matrices, rowI)
 268     {
 269         if
 270         (
 271             Pstream::defaultCommsType() == Pstream::blocking
 272          || Pstream::defaultCommsType() == Pstream::nonBlocking
 273         )
 274         {
 275             forAll (interfaces[rowI], interfaceI)
 276             {
 277                 if (interfaces[rowI].set(interfaceI))
 278                 {
 279                     if
 280                     (
 281                         isA<processorLduInterfaceField>
 282                         (
 283                             interfaces[rowI][interfaceI]
 284                         )
 285                     )
 286                     {
 287                         interfaces[rowI][interfaceI].initInterfaceMatrixUpdate
 288                         (
 289                             x[rowI],
 290                             result[rowI],
 291                             matrices[rowI],
 292                             coupleCoeffs[rowI][interfaceI],
 293                             cmpt,
 294                             static_cast<const Pstream::commsTypes>
 295                             (
 296                                 Pstream::defaultCommsType()
 297                             ),
 298                             false
 299                         );
 300                     }
 301                 }
 302             }
 303         }
 304     }
 305 }
 306
 307
 308 void Foam::coupledLduMatrix::updateMatrixInterfaces
 309 (
 310     const PtrList<FieldField<Field, scalar> >& coupleCoeffs,
 311     const lduInterfaceFieldPtrsListList& interfaces,
 312     const FieldField<Field, scalar>& x,
 313     FieldField<Field, scalar>& result,
 314     const direction cmpt
 315 ) const
 316 {
 317     const PtrList<lduMatrix>& matrices = *this;
 318
 319     forAll (matrices, rowI)
 320     {
 321         matrices[rowI].updateMatrixInterfaces
 322         (
 323             coupleCoeffs[rowI],
 324             interfaces[rowI],
 325             x[rowI],
 326             result[rowI],
 327             cmpt
 328         );
 329     }
 330 }
 331
 332
 333 // ************************************************************************* //