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 bool switchToLhs
 200 ) const
 201 {
 202     const PtrList<lduMatrix>& matrices = *this;
 203
 204     // Note.  The comms design requires all non-processor interfaces
 205     // to be updated first, followed by the update of processor
 206     // interfaces.  The reason is that non-processor coupled
 207     // interfaces require a complex comms pattern involving more than
 208     // pairwise communications.
 209     // Under normal circumstances this is achieved naturally, since
 210     // processor interfaces come last on the list and other coupled
 211     // interfaces execute complex comms at init() level.
 212     // For coupled matrices, the update loop needs to be split over
 213     // all matrices by hand
 214     // Bug fix: Zeljko Tukovic, 7/Apr/2015
 215
 216     // Init update all non-processor coupled interfaces
 217     forAll (matrices, rowI)
 218     {
 219         if
 220         (
 221             Pstream::defaultCommsType() == Pstream::blocking
 222          || Pstream::defaultCommsType() == Pstream::nonBlocking
 223         )
 224         {
 225             forAll (interfaces[rowI], interfaceI)
 226             {
 227                 if (interfaces[rowI].set(interfaceI))
 228                 {
 229                     if
 230                     (
 231                        !isA<processorLduInterfaceField>
 232                         (
 233                             interfaces[rowI][interfaceI]
 234                         )
 235                     )
 236                     {
 237                         interfaces[rowI][interfaceI].initInterfaceMatrixUpdate
 238                         (
 239                             x[rowI],
 240                             result[rowI],
 241                             matrices[rowI],
 242                             coupleCoeffs[rowI][interfaceI],
 243                             cmpt,
 244                             static_cast<const Pstream::commsTypes>
 245                             (
 246                                 Pstream::defaultCommsType()
 247                             ),
 248                             switchToLhs
 249                         );
 250                     }
 251                 }
 252             }
 253         }
 254         else if (Pstream::defaultCommsType() == Pstream::scheduled)
 255         {
 256             // ERROR: Does not work with scheduled comms.
 257             // To investigate.  HJ, 11/Jun/2015
 258             forAll (matrices, rowI)
 259             {
 260                 const lduSchedule& patchSchedule =
 261                     matrices[rowI].patchSchedule();
 262
 263                 // Loop over the "global" patches are on the list of
 264                 // interfaces but beyond the end of the schedule
 265                 // which only handles "normal" patches
 266                 for
 267                 (
 268                     label interfaceI = patchSchedule.size()/2;
 269                     interfaceI < interfaces.size();
 270                     interfaceI++
 271                 )
 272                 {
 273                     if (interfaces[rowI].set(interfaceI))
 274                     {
 275                         if
 276                         (
 277                             !isA<processorLduInterfaceField>
 278                             (
 279                                 interfaces[rowI][interfaceI]
 280                             )
 281                         )
 282                         {
 283                             interfaces[rowI][interfaceI].
 284                                 initInterfaceMatrixUpdate
 285                                 (
 286                                     x[rowI],
 287                                     result[rowI],
 288                                     matrices[rowI],
 289                                     coupleCoeffs[rowI][interfaceI],
 290                                     cmpt,
 291                                     static_cast<const Pstream::commsTypes>
 292                                     (
 293                                         Pstream::defaultCommsType()
 294                                     ),
 295                                     switchToLhs
 296                                 );
 297                         }
 298                     }
 299                 }
 300             }
 301         }
 302         else
 303         {
 304             FatalErrorIn("void coupledLduMatrix::initMatrixInterfaces")
 305                 << "Unsuported communications type "
 306                 << Pstream::commsTypeNames[Pstream::defaultCommsType()]
 307                 << exit(FatalError);
 308         }
 309     }
 310
 311     // Init update for all processor interfaces
 312     forAll (matrices, rowI)
 313     {
 314         if
 315         (
 316             Pstream::defaultCommsType() == Pstream::blocking
 317          || Pstream::defaultCommsType() == Pstream::nonBlocking
 318         )
 319         {
 320             forAll (interfaces[rowI], interfaceI)
 321             {
 322                 if (interfaces[rowI].set(interfaceI))
 323                 {
 324                     if
 325                     (
 326                         isA<processorLduInterfaceField>
 327                         (
 328                             interfaces[rowI][interfaceI]
 329                         )
 330                     )
 331                     {
 332                         interfaces[rowI][interfaceI].initInterfaceMatrixUpdate
 333                         (
 334                             x[rowI],
 335                             result[rowI],
 336                             matrices[rowI],
 337                             coupleCoeffs[rowI][interfaceI],
 338                             cmpt,
 339                             static_cast<const Pstream::commsTypes>
 340                             (
 341                                 Pstream::defaultCommsType()
 342                             ),
 343                             switchToLhs
 344                         );
 345                     }
 346                 }
 347             }
 348         }
 349         else if (Pstream::defaultCommsType() == Pstream::scheduled)
 350         {
 351             // ERROR: Does not work with scheduled comms.
 352             // To investigate.  HJ, 11/Jun/2015
 353             forAll (matrices, rowI)
 354             {
 355                 const lduSchedule& patchSchedule =
 356                     matrices[rowI].patchSchedule();
 357
 358                 // Loop over the "global" patches are on the list of
 359                 // interfaces but beyond the end of the schedule
 360                 // which only handles "normal" patches
 361                 for
 362                 (
 363                     label interfaceI = patchSchedule.size()/2;
 364                     interfaceI < interfaces.size();
 365                     interfaceI++
 366                 )
 367                 {
 368                     if (interfaces[rowI].set(interfaceI))
 369                     {
 370                         if
 371                         (
 372                             isA<processorLduInterfaceField>
 373                             (
 374                                 interfaces[rowI][interfaceI]
 375                             )
 376                         )
 377                         {
 378                             interfaces[rowI][interfaceI].
 379                                 initInterfaceMatrixUpdate
 380                                 (
 381                                     x[rowI],
 382                                     result[rowI],
 383                                     matrices[rowI],
 384                                     coupleCoeffs[rowI][interfaceI],
 385                                     cmpt,
 386                                     static_cast<const Pstream::commsTypes>
 387                                     (
 388                                         Pstream::defaultCommsType()
 389                                     ),
 390                                     switchToLhs
 391                                 );
 392                         }
 393                     }
 394                 }
 395             }
 396         }
 397         else
 398         {
 399             FatalErrorIn("void coupledLduMatrix::initMatrixInterfaces")
 400                 << "Unsuported communications type "
 401                 << Pstream::commsTypeNames[Pstream::defaultCommsType()]
 402                 << exit(FatalError);
 403         }
 404     }
 405 }
 406
 407
 408 void Foam::coupledLduMatrix::updateMatrixInterfaces
 409 (
 410     const PtrList<FieldField<Field, scalar> >& coupleCoeffs,
 411     const lduInterfaceFieldPtrsListList& interfaces,
 412     const FieldField<Field, scalar>& x,
 413     FieldField<Field, scalar>& result,
 414     const direction cmpt,
 415     const bool switchToLhs
 416 ) const
 417 {
 418     const PtrList<lduMatrix>& matrices = *this;
 419
 420     forAll (matrices, rowI)
 421     {
 422         matrices[rowI].updateMatrixInterfaces
 423         (
 424             coupleCoeffs[rowI],
 425             interfaces[rowI],
 426             x[rowI],
 427             result[rowI],
 428             cmpt,
 429             switchToLhs
 430         );
 431     }
 432 }
 433
 434
 435 // ************************************************************************* //