src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 Description
  26     Multiply a given vector (second argument) by the matrix or its transpose
  27     and return the result in the first argument.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "lduMatrix.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 void Foam::lduMatrix::Amul
  36 (
  37     scalarField& Apsi,
  38     const tmp<scalarField>& tpsi,
  39     const FieldField<Field, scalar>& interfaceBouCoeffs,
  40     const lduInterfaceFieldPtrsList& interfaces,
  41     const direction cmpt
  42 ) const
  43 {
  44     scalar* __restrict__ ApsiPtr = Apsi.begin();
  45
  46     const scalarField& psi = tpsi();
  47     const scalar* const __restrict__ psiPtr = psi.begin();
  48
  49     const scalar* const __restrict__ diagPtr = diag().begin();
  50
  51     const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
  52     const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
  53
  54     const scalar* const __restrict__ upperPtr = upper().begin();
  55     const scalar* const __restrict__ lowerPtr = lower().begin();
  56
  57     // Initialise the update of interfaced interfaces
  58     initMatrixInterfaces
  59     (
  60         interfaceBouCoeffs,
  61         interfaces,
  62         psi,
  63         Apsi,
  64         cmpt
  65     );
  66
  67     register const label nCells = diag().size();
  68     for (register label cell=0; cell<nCells; cell++)
  69     {
  70         #ifdef ICC_IA64_PREFETCH
  71         __builtin_prefetch (&psiPtr[cell+96],0,1);
  72         __builtin_prefetch (&diagPtr[cell+96],0,1);
  73         __builtin_prefetch (&ApsiPtr[cell+96],1,1);
  74         #endif
  75
  76         ApsiPtr[cell] = diagPtr[cell]*psiPtr[cell];
  77     }
  78
  79
  80     register const label nFaces = upper().size();
  81     #ifdef ICC_IA64_PREFETCH
  82     #pragma swp
  83     #endif
  84     for (register label face=0; face<nFaces; face++)
  85     {
  86         #ifdef ICC_IA64_PREFETCH
  87         __builtin_prefetch (&uPtr[face+32],0,0);
  88         __builtin_prefetch (&lPtr[face+32],0,0);
  89         __builtin_prefetch (&lowerPtr[face+32],0,1);
  90         __builtin_prefetch (&psiPtr[lPtr[face+32]],0,1);
  91         __builtin_prefetch (&ApsiPtr[uPtr[face+32]],0,1);
  92         #endif
  93
  94         ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];
  95
  96         #ifdef ICC_IA64_PREFETCH
  97         __builtin_prefetch (&upperPtr[face+32],0,1);
  98         __builtin_prefetch (&psiPtr[uPtr[face+32]],0,1);
  99         __builtin_prefetch (&ApsiPtr[lPtr[face+32]],0,1);
 100         #endif
 101
 102         ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];
 103     }
 104
 105     // Update interface interfaces
 106     updateMatrixInterfaces
 107     (
 108         interfaceBouCoeffs,
 109         interfaces,
 110         psi,
 111         Apsi,
 112         cmpt
 113     );
 114
 115     tpsi.clear();
 116 }
 117
 118
 119 void Foam::lduMatrix::Tmul
 120 (
 121     scalarField& Tpsi,
 122     const tmp<scalarField>& tpsi,
 123     const FieldField<Field, scalar>& interfaceIntCoeffs,
 124     const lduInterfaceFieldPtrsList& interfaces,
 125     const direction cmpt
 126 ) const
 127 {
 128     scalar* __restrict__ TpsiPtr = Tpsi.begin();
 129
 130     const scalarField& psi = tpsi();
 131     const scalar* const __restrict__ psiPtr = psi.begin();
 132
 133     const scalar* const __restrict__ diagPtr = diag().begin();
 134
 135     const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
 136     const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
 137
 138     const scalar* const __restrict__ lowerPtr = lower().begin();
 139     const scalar* const __restrict__ upperPtr = upper().begin();
 140
 141     // Initialise the update of interfaced interfaces
 142     initMatrixInterfaces
 143     (
 144         interfaceIntCoeffs,
 145         interfaces,
 146         psi,
 147         Tpsi,
 148         cmpt
 149     );
 150
 151     register const label nCells = diag().size();
 152     for (register label cell=0; cell<nCells; cell++)
 153     {
 154         #ifdef ICC_IA64_PREFETCH
 155         __builtin_prefetch (&psiPtr[cell+96],0,1);
 156         __builtin_prefetch (&diagPtr[cell+96],0,1);
 157         __builtin_prefetch (&TpsiPtr[cell+96],1,1);
 158         #endif
 159
 160         TpsiPtr[cell] = diagPtr[cell]*psiPtr[cell];
 161     }
 162
 163     register const label nFaces = upper().size();
 164     for (register label face=0; face<nFaces; face++)
 165     {
 166         #ifdef ICC_IA64_PREFETCH
 167         __builtin_prefetch (&uPtr[face+32],0,0);
 168         __builtin_prefetch (&lPtr[face+32],0,0);
 169         __builtin_prefetch (&upperPtr[face+32],0,1);
 170         __builtin_prefetch (&psiPtr[lPtr[face+32]],0,1);
 171         __builtin_prefetch (&TpsiPtr[uPtr[face+32]],0,1);
 172         #endif
 173
 174         TpsiPtr[uPtr[face]] += upperPtr[face]*psiPtr[lPtr[face]];
 175
 176         #ifdef ICC_IA64_PREFETCH
 177         __builtin_prefetch (&lowerPtr[face+32],0,1);
 178         __builtin_prefetch (&psiPtr[uPtr[face+32]],0,1);
 179         __builtin_prefetch (&TpsiPtr[lPtr[face+32]],0,1);
 180         #endif
 181
 182         TpsiPtr[lPtr[face]] += lowerPtr[face]*psiPtr[uPtr[face]];
 183     }
 184
 185     // Update interface interfaces
 186     updateMatrixInterfaces
 187     (
 188         interfaceIntCoeffs,
 189         interfaces,
 190         psi,
 191         Tpsi,
 192         cmpt
 193     );
 194
 195     tpsi.clear();
 196 }
 197
 198
 199 void Foam::lduMatrix::sumA
 200 (
 201     scalarField& sumA,
 202     const FieldField<Field, scalar>& interfaceBouCoeffs,
 203     const lduInterfaceFieldPtrsList& interfaces
 204 ) const
 205 {
 206     scalar* __restrict__ sumAPtr = sumA.begin();
 207
 208     const scalar* __restrict__ diagPtr = diag().begin();
 209
 210     const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
 211     const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
 212
 213     const scalar* __restrict__ lowerPtr = lower().begin();
 214     const scalar* __restrict__ upperPtr = upper().begin();
 215
 216     register const label nCells = diag().size();
 217     register const label nFaces = upper().size();
 218
 219     for (register label cell=0; cell<nCells; cell++)
 220     {
 221         #ifdef ICC_IA64_PREFETCH
 222         __builtin_prefetch (&diagPtr[cell+96],0,1);
 223         __builtin_prefetch (&sumAPtr[cell+96],1,1);
 224         #endif
 225
 226         sumAPtr[cell] = diagPtr[cell];
 227     }
 228
 229     #ifdef ICC_IA64_PREFETCH
 230     #pragma swp
 231     #endif
 232
 233     for (register label face=0; face<nFaces; face++)
 234     {
 235         #ifdef ICC_IA64_PREFETCH
 236         __builtin_prefetch (&uPtr[face+32],0,0);
 237         __builtin_prefetch (&lPtr[face+32],0,0);
 238         __builtin_prefetch (&lowerPtr[face+32],0,1);
 239         __builtin_prefetch (&sumAPtr[uPtr[face+32]],0,1);
 240         #endif
 241
 242         sumAPtr[uPtr[face]] += lowerPtr[face];
 243
 244         #ifdef ICC_IA64_PREFETCH
 245         __builtin_prefetch (&upperPtr[face+32],0,1);
 246         __builtin_prefetch (&sumAPtr[lPtr[face+32]],0,1);
 247         #endif
 248
 249         sumAPtr[lPtr[face]] += upperPtr[face];
 250     }
 251
 252     // Add the interface internal coefficients to diagonal
 253     // and the interface boundary coefficients to the sum-off-diagonal
 254     forAll(interfaces, patchI)
 255     {
 256         if (interfaces.set(patchI))
 257         {
 258             const unallocLabelList& pa = lduAddr().patchAddr(patchI);
 259             const scalarField& pCoeffs = interfaceBouCoeffs[patchI];
 260
 261             forAll(pa, face)
 262             {
 263                 sumAPtr[pa[face]] -= pCoeffs[face];
 264             }
 265         }
 266     }
 267 }
 268
 269
 270 void Foam::lduMatrix::residual
 271 (
 272     scalarField& rA,
 273     const scalarField& psi,
 274     const scalarField& source,
 275     const FieldField<Field, scalar>& interfaceBouCoeffs,
 276     const lduInterfaceFieldPtrsList& interfaces,
 277     const direction cmpt
 278 ) const
 279 {
 280     scalar* __restrict__ rAPtr = rA.begin();
 281
 282     const scalar* const __restrict__ psiPtr = psi.begin();
 283     const scalar* const __restrict__ diagPtr = diag().begin();
 284     const scalar* const __restrict__ sourcePtr = source.begin();
 285
 286     const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
 287     const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
 288
 289     const scalar* const __restrict__ upperPtr = upper().begin();
 290     const scalar* const __restrict__ lowerPtr = lower().begin();
 291
 292     // Parallel boundary initialisation.
 293     // Note: there is a change of sign in the coupled
 294     // interface update.  The reason for this is that the
 295     // internal coefficients are all located at the l.h.s. of
 296     // the matrix whereas the "implicit" coefficients on the
 297     // coupled boundaries are all created as if the
 298     // coefficient contribution is of a source-kind (i.e. they
 299     // have a sign as if they are on the r.h.s. of the matrix.
 300     // To compensate for this, it is necessary to turn the
 301     // sign of the contribution.
 302
 303     FieldField<Field, scalar> mBouCoeffs(interfaceBouCoeffs.size());
 304
 305     forAll(mBouCoeffs, patchi)
 306     {
 307         if (interfaces.set(patchi))
 308         {
 309             mBouCoeffs.set(patchi, -interfaceBouCoeffs[patchi]);
 310         }
 311     }
 312
 313     // Initialise the update of interfaced interfaces
 314     initMatrixInterfaces
 315     (
 316         mBouCoeffs,
 317         interfaces,
 318         psi,
 319         rA,
 320         cmpt
 321     );
 322
 323     register const label nCells = diag().size();
 324     for (register label cell=0; cell<nCells; cell++)
 325     {
 326         #ifdef ICC_IA64_PREFETCH
 327         __builtin_prefetch (&psiPtr[cell+96],0,1);
 328         __builtin_prefetch (&diagPtr[cell+96],0,1);
 329         __builtin_prefetch (&sourcePtr[cell+96],0,1);
 330         __builtin_prefetch (&rAPtr[cell+96],1,1);
 331         #endif
 332
 333         rAPtr[cell] = sourcePtr[cell] - diagPtr[cell]*psiPtr[cell];
 334     }
 335
 336
 337     register const label nFaces = upper().size();
 338     #ifdef ICC_IA64_PREFETCH
 339     #pragma swp
 340     #endif
 341     for (register label face=0; face<nFaces; face++)
 342     {
 343         #ifdef ICC_IA64_PREFETCH
 344         __builtin_prefetch (&uPtr[face+32],0,0);
 345         __builtin_prefetch (&lPtr[face+32],0,0);
 346         __builtin_prefetch (&lowerPtr[face+32],0,1);
 347         __builtin_prefetch (&psiPtr[lPtr[face+32]],0,1);
 348         __builtin_prefetch (&rAPtr[uPtr[face+32]],0,1);
 349         #endif
 350
 351         rAPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
 352
 353         #ifdef ICC_IA64_PREFETCH
 354         __builtin_prefetch (&upperPtr[face+32],0,1);
 355         __builtin_prefetch (&psiPtr[uPtr[face+32]],0,1);
 356         __builtin_prefetch (&rAPtr[lPtr[face+32]],0,1);
 357         #endif
 358
 359         rAPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
 360     }
 361
 362     // Update interface interfaces
 363     updateMatrixInterfaces
 364     (
 365         mBouCoeffs,
 366         interfaces,
 367         psi,
 368         rA,
 369         cmpt
 370     );
 371 }
 372
 373
 374 Foam::tmp<Foam::scalarField> Foam::lduMatrix::residual
 375 (
 376     const scalarField& psi,
 377     const scalarField& source,
 378     const FieldField<Field, scalar>& interfaceBouCoeffs,
 379     const lduInterfaceFieldPtrsList& interfaces,
 380     const direction cmpt
 381 ) const
 382 {
 383     tmp<scalarField> trA(new scalarField(psi.size()));
 384     residual(trA(), psi, source, interfaceBouCoeffs, interfaces, cmpt);
 385     return trA;
 386 }
 387
 388
 389 // ************************************************************************* //