src/lduSolvers/amg/amgPolicy/pamgPolicy.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     pamgPolicy
  26
  27 Description
  28     Pairwise agglomerative AMG policy.  Legacy implementation
  29
  30 Author
  31     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "pamgPolicy.H"
  36 #include "amgMatrix.H"
  37 #include "boolList.H"
  38 #include "addToRunTimeSelectionTable.H"
  39 #include "GAMGInterfaceField.H"
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45     defineTypeNameAndDebug(pamgPolicy, 0);
  46
  47     addToRunTimeSelectionTable(amgPolicy, pamgPolicy, matrix);
  48
  49 } // End namespace Foam
  50
  51
  52 const Foam::debug::tolerancesSwitch
  53 Foam::pamgPolicy::diagFactor_
  54 (
  55     "pamgDiagFactor",
  56     1e-8
  57 );
  58
  59
  60 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  61
  62 void Foam::pamgPolicy::calcChild()
  63 {
  64     // Algorithm:
  65     // 1) Create temporary equation addressing using a double-pass algorithm.
  66     //    to create the offset table.
  67     // 2) Loop through all equations and for each equation find the best fit
  68     //    neighbour.  If all neighbours are grouped, add equation to best group
  69
  70     // Get addressing
  71     const label nEqns = matrix_.lduAddr().size();
  72
  73     const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
  74     const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
  75
  76     // Get off-diagonal matrix coefficients
  77     const scalarField& upper = matrix_.upper();
  78
  79     // For each equation calculate coeffs
  80     labelList cols(upperAddr.size() + lowerAddr.size());
  81     labelList rowOffsets(nEqns + 1);
  82
  83     // Memory management
  84     {
  85         labelList nNbrs(nEqns, 0);
  86
  87         forAll (upperAddr, coeffI)
  88         {
  89             nNbrs[upperAddr[coeffI]]++;
  90         }
  91
  92         forAll (lowerAddr, coeffI)
  93         {
  94             nNbrs[lowerAddr[coeffI]]++;
  95         }
  96
  97         rowOffsets[0] = 0;
  98
  99         forAll (nNbrs, eqnI)
 100         {
 101             rowOffsets[eqnI + 1] = rowOffsets[eqnI] + nNbrs[eqnI];
 102         }
 103
 104         // Reset the list to use as counter
 105         nNbrs = 0;
 106
 107         forAll (upperAddr, coeffI)
 108         {
 109             cols
 110             [
 111                 rowOffsets[upperAddr[coeffI]] + nNbrs[upperAddr[coeffI]]
 112             ] = coeffI;
 113
 114             nNbrs[upperAddr[coeffI]]++;
 115         }
 116
 117         forAll (lowerAddr, coeffI)
 118         {
 119             cols
 120             [
 121                 rowOffsets[lowerAddr[coeffI]] + nNbrs[lowerAddr[coeffI]]
 122             ] = coeffI;
 123
 124             nNbrs[lowerAddr[coeffI]]++;
 125         }
 126     }
 127
 128     // Reset child array
 129     child_.setSize(nEqns, -1);
 130
 131     // Calculate agglomeration
 132
 133     // Get matrix coefficients
 134     const scalarField& diag = matrix_.diag();
 135
 136     scalarField magOffDiag;
 137
 138     if (matrix_.asymmetric())
 139     {
 140         magOffDiag = Foam::max(mag(matrix_.upper()), mag(matrix_.lower()));
 141     }
 142     else if (matrix_.symmetric())
 143     {
 144         magOffDiag = mag(matrix_.upper());
 145     }
 146     else
 147     {
 148         // Diag only matrix.  Reset and return
 149         child_ = 0;
 150         nCoarseEqns_ = 1;
 151
 152         return;
 153     }
 154
 155     nCoarseEqns_ = 0;
 156
 157     // Gather disconnected and weakly connected equations into cluster zero
 158     // Weak connection is assumed to be the one where the off-diagonal
 159     // coefficient is smaller than diagFactor_*diag
 160     {
 161         // Algorithm
 162         // Mark all cells to belong to zero cluster
 163         // Go through all upper and lower coefficients and for the ones
 164         // larger than threshold mark the equations out of cluster zero
 165
 166         scalarField magScaledDiag = diagFactor_()*mag(diag);
 167
 168         boolList zeroCluster(diag.size(), true);
 169
 170         forAll (magOffDiag, coeffI)
 171         {
 172             if (magOffDiag[coeffI] > magScaledDiag[upperAddr[coeffI]])
 173             {
 174                 zeroCluster[upperAddr[coeffI]] = false;
 175             }
 176
 177             if (magOffDiag[coeffI] > magScaledDiag[lowerAddr[coeffI]])
 178             {
 179                 zeroCluster[lowerAddr[coeffI]] = false;
 180             }
 181         }
 182
 183         // Collect solo equations
 184         forAll (zeroCluster, eqnI)
 185         {
 186             if (zeroCluster[eqnI])
 187             {
 188                 // Found solo equation
 189                 nSolo_++;
 190
 191                 child_[eqnI] = nCoarseEqns_;
 192             }
 193         }
 194
 195         if (nSolo_ > 0)
 196         {
 197             // Found solo equations
 198             nCoarseEqns_++;
 199
 200             if (lduMatrix::debug >= 2)
 201             {
 202                 Pout<< "Found " << nSolo_ << " weakly connected equations."
 203                     << endl;
 204             }
 205         }
 206     }
 207
 208     // Go through the off-diagonal and create clusters, marking the child array
 209     for (label eqnI = 0; eqnI < nEqns; eqnI++)
 210     {
 211         if (child_[eqnI] < 0)
 212         {
 213             label matchCoeffNo = -1;
 214             scalar maxCoeff = -GREAT;
 215
 216             // Check row to find ungrouped neighbour with largest coefficient
 217             for
 218             (
 219                 label rowCoeffI = rowOffsets[eqnI];
 220                 rowCoeffI < rowOffsets[eqnI + 1];
 221                 rowCoeffI++
 222             )
 223             {
 224                 label coeffI = cols[rowCoeffI];
 225
 226                 // I don't know whether the current equation is owner
 227                 //  or neighbour.  Therefore I'll check both sides
 228                 if
 229                 (
 230                     child_[upperAddr[coeffI]] < 0
 231                  && child_[lowerAddr[coeffI]] < 0
 232                  && mag(upper[coeffI]) > maxCoeff
 233                 )
 234                 {
 235                     // Match found. Pick up all the necessary data
 236                     matchCoeffNo = coeffI;
 237                     maxCoeff = mag(upper[coeffI]);
 238                 }
 239             }
 240
 241             if (matchCoeffNo >= 0)
 242             {
 243                 // Make a new group
 244                 child_[upperAddr[matchCoeffNo]] = nCoarseEqns_;
 245                 child_[lowerAddr[matchCoeffNo]] = nCoarseEqns_;
 246                 nCoarseEqns_++;
 247             }
 248             else
 249             {
 250                 // No match. Find the best neighbouring cluster and
 251                 // put the equation there
 252                 label clusterMatchCoeffNo = -1;
 253                 scalar clusterMaxCoeff = -GREAT;
 254
 255                 for
 256                 (
 257                     label rowCoeffI = rowOffsets[eqnI];
 258                     rowCoeffI < rowOffsets[eqnI + 1];
 259                     rowCoeffI++
 260                 )
 261                 {
 262                     label coeffI = cols[rowCoeffI];
 263
 264                     if (mag(upper[coeffI]) > clusterMaxCoeff)
 265                     {
 266                         clusterMatchCoeffNo = coeffI;
 267                         clusterMaxCoeff = mag(upper[coeffI]);
 268                     }
 269                 }
 270
 271                 if (clusterMatchCoeffNo >= 0)
 272                 {
 273                     // Add the equation to the best cluster
 274                     child_[eqnI] = max
 275                     (
 276                         child_[upperAddr[clusterMatchCoeffNo]],
 277                         child_[lowerAddr[clusterMatchCoeffNo]]
 278                     );
 279                 }
 280                 else
 281                 {
 282                     // This is a Dirichlet point: no neighbours
 283                     // Put it into its own cluster
 284                     child_[eqnI] = nCoarseEqns_;
 285                     nCoarseEqns_++;
 286                 }
 287             }
 288         }
 289     }
 290
 291     // Reverse the map to facilitate later agglomeration.
 292     // Keep zero cluster at zero index
 293     forAll (child_, eqnI)
 294     {
 295         child_[eqnI] = nCoarseEqns_ - 1 - child_[eqnI];
 296     }
 297
 298     // The decision on parallel agglomeration needs to be made for the
 299     // whole gang of processes; otherwise I may end up with a different
 300     // number of agglomeration levels on different processors.
 301
 302     if (nCoarseEqns_ > minCoarseEqns() && 3*nCoarseEqns_ <= 2*nEqns)
 303     {
 304         coarsen_ = true;
 305     }
 306
 307     reduce(coarsen_, andOp<bool>());
 308
 309     if (lduMatrix::debug >= 2)
 310     {
 311         Pout << "Coarse level size: " << nCoarseEqns_;
 312
 313         if (coarsen_)
 314         {
 315             Pout << ".  Accepted" << endl;
 316         }
 317         else
 318         {
 319             Pout << ".  Rejected" << endl;
 320         }
 321     }
 322 }
 323
 324
 325 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 326
 327 Foam::pamgPolicy::pamgPolicy
 328 (
 329     const lduMatrix& matrix,
 330     const label groupSize,
 331     const label minCoarseEqns
 332 )
 333 :
 334     amgPolicy(groupSize, minCoarseEqns),
 335     matrix_(matrix),
 336     child_(),
 337     nSolo_(0),
 338     nCoarseEqns_(0),
 339     coarsen_(false)
 340 {
 341     calcChild();
 342 }
 343
 344 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 345
 346 Foam::pamgPolicy::~pamgPolicy()
 347 {}
 348
 349
 350 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 351
 352 Foam::autoPtr<Foam::amgMatrix> Foam::pamgPolicy::restrictMatrix
 353 (
 354     const FieldField<Field, scalar>& bouCoeffs,
 355     const FieldField<Field, scalar>& intCoeffs,
 356     const lduInterfaceFieldPtrsList& interfaceFields
 357 ) const
 358 {
 359     if (!coarsen_)
 360     {
 361         FatalErrorIn("autoPtr<amgMatrix> pamgPolicy::restrictMatrix() const")
 362             << "Requesting coarse matrix when it cannot be created"
 363             << abort(FatalError);
 364     }
 365
 366     // Construct the coarse matrix and ldu addressing for the next level
 367     // Algorithm:
 368     // 1) Loop through all fine coeffs. If the child labels on two sides are
 369     //    different, this creates a coarse coeff. Define owner and neighbour
 370     //    for this coeff based on cluster IDs.
 371     // 2) Check if the coeff has been seen before. If yes, add the coefficient
 372     //    to the appropriate field (stored with the equation). If no, create
 373     //    a new coeff with neighbour ID and add the coefficient
 374     // 3) Once all the coeffs have been created, loop through all clusters and
 375     //    insert the coeffs in the upper order. At the same time, collect the
 376     //    owner and neighbour addressing.
 377     // 4) Agglomerate the diagonal by summing up the fine diagonal
 378
 379     // Get addressing
 380     const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
 381     const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
 382
 383     label nFineCoeffs = upperAddr.size();
 384
 385 #   ifdef FULLDEBUG
 386     if (child_.size() != matrix_.lduAddr().size())
 387     {
 388         FatalErrorIn
 389         (
 390             "autoPtr<amgMatrix> pamgPolicy::restrictMatrix() const"
 391         )   << "Child array does not correspond to fine level. " << endl
 392             << " Child size: " << child_.size()
 393             << " number of equations: " << matrix_.lduAddr().size()
 394             << abort(FatalError);
 395     }
 396 #   endif
 397
 398
 399     // Does the matrix have solo equations
 400     bool soloEqns = nSolo_ > 0;
 401
 402     // Storage for block neighbours and coefficients
 403
 404     // Guess initial maximum number of neighbours in block
 405     label maxNnbrs = 10;
 406
 407     // Number of neighbours per block
 408     labelList blockNnbrs(nCoarseEqns_, 0);
 409
 410     // Setup initial packed storage for neighbours and coefficients
 411     labelList blockNbrsData(maxNnbrs*nCoarseEqns_);
 412
 413     scalarList blockCoeffsData(maxNnbrs*nCoarseEqns_, 0.0);
 414
 415     // Create face-restriction addressing
 416     // Note: value of coeffRestrictAddr for off-diagonal coefficients
 417     // touching solo cells will be invalid
 418     // HJ, 7/Apr/2015
 419     labelList coeffRestrictAddr(nFineCoeffs);
 420
 421     // Initial neighbour array (not in upper-triangle order)
 422     labelList initCoarseNeighb(nFineCoeffs);
 423
 424     // Counter for coarse coeffs
 425     label nCoarseCoeffs = 0;
 426
 427     // Loop through all fine coeffs
 428     forAll (upperAddr, fineCoeffI)
 429     {
 430         label rmUpperAddr = child_[upperAddr[fineCoeffI]];
 431         label rmLowerAddr = child_[lowerAddr[fineCoeffI]];
 432
 433         // If the coefficient touches block zero and solo equations are
 434         // present, skip it
 435         if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0))
 436         {
 437             continue;
 438         }
 439
 440         if (rmUpperAddr == rmLowerAddr)
 441         {
 442             // For each fine coeff inside of a coarse cluster keep the address
 443             // of the cluster corresponding to the coeff in the
 444             // coeffRestrictAddr as a negative index
 445             coeffRestrictAddr[fineCoeffI] = -(rmUpperAddr + 1);
 446         }
 447         else
 448         {
 449             // This coeff is a part of a coarse coeff
 450
 451             label cOwn = rmUpperAddr;
 452             label cNei = rmLowerAddr;
 453
 454             // Get coarse owner and neighbour
 455             if (rmUpperAddr > rmLowerAddr)
 456             {
 457                 cOwn = rmLowerAddr;
 458                 cNei = rmUpperAddr;
 459             }
 460
 461             // Check the neighbour to see if this coeff has already been found
 462             label* ccCoeffs = &blockNbrsData[maxNnbrs*cOwn];
 463
 464             bool nbrFound = false;
 465             label& ccnCoeffs = blockNnbrs[cOwn];
 466
 467             for (int i = 0; i < ccnCoeffs; i++)
 468             {
 469                 if (initCoarseNeighb[ccCoeffs[i]] == cNei)
 470                 {
 471                     nbrFound = true;
 472                     coeffRestrictAddr[fineCoeffI] = ccCoeffs[i];
 473                     break;
 474                 }
 475             }
 476
 477             if (!nbrFound)
 478             {
 479                 if (ccnCoeffs >= maxNnbrs)
 480                 {
 481                     label oldMaxNnbrs = maxNnbrs;
 482                     maxNnbrs *= 2;
 483
 484                     blockNbrsData.setSize(maxNnbrs*nCoarseEqns_);
 485
 486                     forAllReverse(blockNnbrs, i)
 487                     {
 488                         label* oldCcNbrs = &blockNbrsData[oldMaxNnbrs*i];
 489                         label* newCcNbrs = &blockNbrsData[maxNnbrs*i];
 490
 491                         for (int j=0; j<blockNnbrs[i]; j++)
 492                         {
 493                             newCcNbrs[j] = oldCcNbrs[j];
 494                         }
 495                     }
 496
 497                     ccCoeffs = &blockNbrsData[maxNnbrs*cOwn];
 498                 }
 499
 500                 ccCoeffs[ccnCoeffs] = nCoarseCoeffs;
 501                 initCoarseNeighb[nCoarseCoeffs] = cNei;
 502                 coeffRestrictAddr[fineCoeffI] = nCoarseCoeffs;
 503                 ccnCoeffs++;
 504
 505                 // New coarse coeff created
 506                 nCoarseCoeffs++;
 507             }
 508         }
 509     } // End for all fine coeffs
 510
 511
 512     // Renumber into upper-triangular order
 513
 514     // All coarse owner-neighbour storage
 515     labelList coarseOwner(nCoarseCoeffs);
 516     labelList coarseNeighbour(nCoarseCoeffs);
 517     labelList coarseCoeffMap(nCoarseCoeffs);
 518
 519     label coarseCoeffi = 0;
 520
 521     forAll (blockNnbrs, cci)
 522     {
 523         label* cCoeffs = &blockNbrsData[maxNnbrs*cci];
 524         label ccnCoeffs = blockNnbrs[cci];
 525
 526         for (int i = 0; i < ccnCoeffs; i++)
 527         {
 528             coarseOwner[coarseCoeffi] = cci;
 529             coarseNeighbour[coarseCoeffi] = initCoarseNeighb[cCoeffs[i]];
 530             coarseCoeffMap[cCoeffs[i]] = coarseCoeffi;
 531             coarseCoeffi++;
 532         }
 533     }
 534
 535     forAll (coeffRestrictAddr, fineCoeffI)
 536     {
 537         label rmUpperAddr = child_[upperAddr[fineCoeffI]];
 538         label rmLowerAddr = child_[lowerAddr[fineCoeffI]];
 539
 540         // If the coefficient touches block zero and solo equations are
 541         // present, skip it
 542         if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0))
 543         {
 544             continue;
 545         }
 546
 547         if (coeffRestrictAddr[fineCoeffI] >= 0)
 548         {
 549             coeffRestrictAddr[fineCoeffI] =
 550                 coarseCoeffMap[coeffRestrictAddr[fineCoeffI]];
 551         }
 552     }
 553
 554     // Clear the temporary storage for the coarse matrix data
 555     blockNnbrs.setSize(0);
 556     blockNbrsData.setSize(0);
 557     initCoarseNeighb.setSize(0);
 558     coarseCoeffMap.setSize(0);
 559
 560
 561     // Create coarse-level coupled interfaces
 562
 563     // Create coarse interfaces, addressing and coefficients
 564
 565     // Set the coarse interfaces and coefficients
 566     lduInterfacePtrsList* coarseInterfacesPtr =
 567         new lduInterfacePtrsList(interfaceFields.size());
 568     lduInterfacePtrsList& coarseInterfaces = *coarseInterfacesPtr;
 569
 570     // Set the coarse interfaceFields and coefficients
 571     lduInterfaceFieldPtrsList* coarseInterfaceFieldsPtr =
 572         new lduInterfaceFieldPtrsList(interfaceFields.size());
 573     lduInterfaceFieldPtrsList& coarseInterfaceFields =
 574         *coarseInterfaceFieldsPtr;
 575
 576     FieldField<Field, scalar>* coarseBouCoeffsPtr =
 577         new FieldField<Field, scalar>(interfaceFields.size());
 578     FieldField<Field, scalar>& coarseBouCoeffs = *coarseBouCoeffsPtr;
 579
 580     FieldField<Field, scalar>* coarseIntCoeffsPtr =
 581         new FieldField<Field, scalar>(interfaceFields.size());
 582     FieldField<Field, scalar>& coarseIntCoeffs = *coarseIntCoeffsPtr;
 583
 584     labelListList coarseInterfaceAddr(interfaceFields.size());
 585
 586     // Add the coarse level
 587
 588     // Set the coarse ldu addressing onto the list
 589     lduPrimitiveMesh* coarseAddrPtr =
 590         new lduPrimitiveMesh
 591         (
 592             nCoarseEqns_,
 593             coarseOwner,
 594             coarseNeighbour,
 595             true
 596         );
 597
 598     // Initialise transfer of restrict addressing on the interface
 599     forAll (interfaceFields, intI)
 600     {
 601         if (interfaceFields.set(intI))
 602         {
 603             interfaceFields[intI].coupledInterface().initInternalFieldTransfer
 604             (
 605                 Pstream::blocking,
 606                 child_
 607             );
 608         }
 609     }
 610
 611     // Store coefficients to avoid tangled communications
 612     // HJ, 1/Apr/2009
 613     FieldField<Field, label> fineInterfaceAddr(interfaceFields.size());
 614
 615     forAll (interfaceFields, intI)
 616     {
 617         if (interfaceFields.set(intI))
 618         {
 619             const lduInterface& fineInterface =
 620                 interfaceFields[intI].coupledInterface();
 621
 622             fineInterfaceAddr.set
 623             (
 624                 intI,
 625                 new labelField
 626                 (
 627                     fineInterface.internalFieldTransfer
 628                     (
 629                         Pstream::blocking,
 630                         child_
 631                     )
 632                 )
 633             );
 634         }
 635     }
 636
 637     // Create GAMG interfaces
 638     forAll (interfaceFields, intI)
 639     {
 640         if (interfaceFields.set(intI))
 641         {
 642             const lduInterface& fineInterface =
 643                 interfaceFields[intI].coupledInterface();
 644
 645             coarseInterfaces.set
 646             (
 647                 intI,
 648                 GAMGInterface::New
 649                 (
 650                     *coarseAddrPtr,
 651                     fineInterface,
 652                     fineInterface.interfaceInternalField(child_),
 653                     fineInterfaceAddr[intI]
 654                 ).ptr()
 655             );
 656         }
 657     }
 658
 659     forAll (interfaceFields, intI)
 660     {
 661         if (interfaceFields.set(intI))
 662         {
 663             const GAMGInterface& coarseInterface =
 664                 refCast<const GAMGInterface>(coarseInterfaces[intI]);
 665
 666             coarseInterfaceFields.set
 667             (
 668                 intI,
 669                 GAMGInterfaceField::New
 670                 (
 671                     coarseInterface,
 672                     interfaceFields[intI]
 673                 ).ptr()
 674             );
 675
 676             coarseBouCoeffs.set
 677             (
 678                 intI,
 679                 coarseInterface.agglomerateCoeffs(bouCoeffs[intI])
 680             );
 681
 682             coarseIntCoeffs.set
 683             (
 684                 intI,
 685                 coarseInterface.agglomerateCoeffs(intCoeffs[intI])
 686             );
 687
 688             coarseInterfaceAddr[intI] = coarseInterface.faceCells();
 689         }
 690     }
 691
 692     // Add interfaces
 693     coarseAddrPtr->addInterfaces
 694     (
 695         *coarseInterfacesPtr,
 696         coarseInterfaceAddr,
 697         matrix_.patchSchedule()
 698     );
 699
 700     // Matrix restriction done!
 701
 702     // Set the coarse level matrix
 703     lduMatrix* coarseMatrixPtr = new lduMatrix(*coarseAddrPtr);
 704     lduMatrix& coarseMatrix = *coarseMatrixPtr;
 705
 706     // Coarse matrix diagonal initialised by restricting the
 707     // finer mesh diagonal
 708     scalarField& coarseDiag = coarseMatrix.diag();
 709     restrictResidual(matrix_.diag(), coarseDiag);
 710
 711     // Check if matrix is assymetric and if so agglomerate both upper and lower
 712     // coefficients ...
 713     if (matrix_.hasLower())
 714     {
 715         // Get off-diagonal matrix coefficients
 716         const scalarField& fineUpper = matrix_.upper();
 717         const scalarField& fineLower = matrix_.lower();
 718
 719         // Coarse matrix upper coefficients
 720         scalarField& coarseUpper = coarseMatrix.upper();
 721         scalarField& coarseLower = coarseMatrix.lower();
 722
 723         forAll (coeffRestrictAddr, fineCoeffI)
 724         {
 725             label rmUpperAddr = child_[upperAddr[fineCoeffI]];
 726             label rmLowerAddr = child_[lowerAddr[fineCoeffI]];
 727
 728             // If the coefficient touches block zero and solo equations are
 729             // present, skip it
 730             if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0))
 731             {
 732                 continue;
 733             }
 734
 735             label cCoeff = coeffRestrictAddr[fineCoeffI];
 736
 737             if (cCoeff >= 0)
 738             {
 739                 coarseUpper[cCoeff] += fineUpper[fineCoeffI];
 740                 coarseLower[cCoeff] += fineLower[fineCoeffI];
 741             }
 742             else
 743             {
 744                 // Add the fine face coefficients into the diagonal.
 745                 coarseDiag[-1 - cCoeff] +=
 746                     fineUpper[fineCoeffI] + fineLower[fineCoeffI];
 747             }
 748         }
 749     }
 750     else // ... Otherwise it is symmetric so agglomerate just the upper
 751     {
 752         // Get off-diagonal matrix coefficients
 753         const scalarField& fineUpper = matrix_.upper();
 754
 755         // Coarse matrix upper coefficients
 756         scalarField& coarseUpper = coarseMatrix.upper();
 757
 758         forAll (coeffRestrictAddr, fineCoeffI)
 759         {
 760             label rmUpperAddr = child_[upperAddr[fineCoeffI]];
 761             label rmLowerAddr = child_[lowerAddr[fineCoeffI]];
 762
 763             // If the coefficient touches block zero and solo equations are
 764             // present, skip it
 765             if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0))
 766             {
 767                 continue;
 768             }
 769
 770             label cCoeff = coeffRestrictAddr[fineCoeffI];
 771
 772             if (cCoeff >= 0)
 773             {
 774                 coarseUpper[cCoeff] += fineUpper[fineCoeffI];
 775             }
 776             else
 777             {
 778                 // Add the fine face coefficient into the diagonal.
 779                 coarseDiag[-1 - cCoeff] += 2*fineUpper[fineCoeffI];
 780             }
 781         }
 782     }
 783
 784     // Create and return amgMatrix
 785     return autoPtr<amgMatrix>
 786     (
 787         new amgMatrix
 788         (
 789             coarseAddrPtr,
 790             coarseInterfacesPtr,
 791             coarseMatrixPtr,
 792             coarseBouCoeffsPtr,
 793             coarseIntCoeffsPtr,
 794             coarseInterfaceFieldsPtr
 795         )
 796     );
 797 }
 798
 799
 800 void Foam::pamgPolicy::restrictResidual
 801 (
 802     const scalarField& res,
 803     scalarField& coarseRes
 804 ) const
 805 {
 806     coarseRes = 0;
 807
 808     forAll (res, i)
 809     {
 810         coarseRes[child_[i]] += res[i];
 811     }
 812 }
 813
 814
 815 void Foam::pamgPolicy::prolongateCorrection
 816 (
 817     scalarField& x,
 818     const scalarField& coarseX
 819 ) const
 820 {
 821     forAll (x, i)
 822     {
 823         x[i] += coarseX[child_[i]];
 824     }
 825 }
 826
 827
 828 // ************************************************************************* //