applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "domainDecomposition.H"
  27 #include "decompositionMethod.H"
  28 #include "cpuTime.H"
  29 #include "cellSet.H"
  30 #include "regionSplit.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 void Foam::domainDecomposition::distributeCells()
  35 {
  36     Info<< "\nCalculating distribution of cells" << endl;
  37
  38     cpuTime decompositionTime;
  39
  40
  41     // See if any faces need to have owner and neighbour on same processor
  42     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  43
  44     labelHashSet sameProcFaces;
  45
  46     if (decompositionDict_.found("preservePatches"))
  47     {
  48         wordList pNames(decompositionDict_.lookup("preservePatches"));
  49
  50         Info<< "Keeping owner of faces in patches " << pNames
  51             << " on same processor. This only makes sense for cyclics." << endl;
  52
  53         const polyBoundaryMesh& patches = boundaryMesh();
  54
  55         forAll(pNames, i)
  56         {
  57             const label patchI = patches.findPatchID(pNames[i]);
  58
  59             if (patchI == -1)
  60             {
  61                 FatalErrorIn("domainDecomposition::distributeCells()")
  62                     << "Unknown preservePatch " << pNames[i]
  63                     << endl << "Valid patches are " << patches.names()
  64                     << exit(FatalError);
  65             }
  66
  67             const polyPatch& pp = patches[patchI];
  68
  69             forAll(pp, i)
  70             {
  71                 sameProcFaces.insert(pp.start() + i);
  72             }
  73         }
  74     }
  75     if (decompositionDict_.found("preserveFaceZones"))
  76     {
  77         wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
  78
  79         Info<< "Keeping owner and neighbour of faces in zones " << zNames
  80             << " on same processor" << endl;
  81
  82         const faceZoneMesh& fZones = faceZones();
  83
  84         forAll(zNames, i)
  85         {
  86             label zoneI = fZones.findZoneID(zNames[i]);
  87
  88             if (zoneI == -1)
  89             {
  90                 FatalErrorIn("domainDecomposition::distributeCells()")
  91                     << "Unknown preserveFaceZone " << zNames[i]
  92                     << endl << "Valid faceZones are " << fZones.names()
  93                     << exit(FatalError);
  94             }
  95
  96             const faceZone& fz = fZones[zoneI];
  97
  98             forAll(fz, i)
  99             {
 100                 sameProcFaces.insert(fz[i]);
 101             }
 102         }
 103     }
 104
 105
 106     // Construct decomposition method and either do decomposition on
 107     // cell centres or on agglomeration
 108
 109
 110     autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
 111     (
 112         decompositionDict_
 113     );
 114
 115     if (sameProcFaces.empty())
 116     {
 117         if (decompositionDict_.found("weightField"))
 118         {
 119             word weightName = decompositionDict_.lookup("weightField");
 120
 121             volScalarField weights
 122             (
 123                 IOobject
 124                 (
 125                     weightName,
 126                     time().timeName(),
 127                     *this,
 128                     IOobject::MUST_READ,
 129                     IOobject::NO_WRITE
 130                 ),
 131                 *this
 132             );
 133
 134             cellToProc_ = decomposePtr().decompose
 135             (
 136                 *this,
 137                 cellCentres(),
 138                 weights.internalField()
 139             );
 140         }
 141         else
 142         {
 143             cellToProc_ = decomposePtr().decompose(*this, cellCentres());
 144         }
 145
 146     }
 147     else
 148     {
 149         Info<< "Selected " << sameProcFaces.size()
 150             << " faces whose owner and neighbour cell should be kept on the"
 151             << " same processor" << endl;
 152
 153         // Faces where owner and neighbour are not 'connected' (= all except
 154         // sameProcFaces)
 155         boolList blockedFace(nFaces(), true);
 156
 157         forAllConstIter(labelHashSet, sameProcFaces, iter)
 158         {
 159             blockedFace[iter.key()] = false;
 160         }
 161
 162         // Connect coupled boundary faces
 163         const polyBoundaryMesh& patches =  boundaryMesh();
 164
 165         forAll(patches, patchI)
 166         {
 167             const polyPatch& pp = patches[patchI];
 168
 169             if (pp.coupled())
 170             {
 171                 forAll(pp, i)
 172                 {
 173                     blockedFace[pp.start()+i] = false;
 174                 }
 175             }
 176         }
 177
 178         // Determine global regions, separated by blockedFaces
 179         regionSplit globalRegion(*this, blockedFace);
 180
 181
 182         // Determine region cell centres
 183         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 184
 185         // This just takes the first cell in the region. Otherwise the problem
 186         // is with cyclics - if we'd average the region centre might be
 187         // somewhere in the middle of the domain which might not be anywhere
 188         // near any of the cells.
 189
 190         pointField regionCentres(globalRegion.nRegions(), point::max);
 191
 192         forAll(globalRegion, cellI)
 193         {
 194             label regionI = globalRegion[cellI];
 195
 196             if (regionCentres[regionI] == point::max)
 197             {
 198                 regionCentres[regionI] = cellCentres()[cellI];
 199             }
 200         }
 201
 202         // Do decomposition on agglomeration
 203         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 204         if (decompositionDict_.found("weightField"))
 205         {
 206             scalarField regionWeights(globalRegion.nRegions(), 0);
 207
 208             word weightName = decompositionDict_.lookup("weightField");
 209
 210             volScalarField weights
 211             (
 212                 IOobject
 213                 (
 214                     weightName,
 215                     time().timeName(),
 216                     *this,
 217                     IOobject::MUST_READ,
 218                     IOobject::NO_WRITE
 219                 ),
 220                 *this
 221             );
 222
 223             forAll(globalRegion, cellI)
 224             {
 225                 label regionI = globalRegion[cellI];
 226
 227                 regionWeights[regionI] += weights.internalField()[cellI];
 228             }
 229
 230             cellToProc_ = decomposePtr().decompose
 231             (
 232                 *this,
 233                 globalRegion,
 234                 regionCentres,
 235                 regionWeights
 236             );
 237         }
 238         else
 239         {
 240             cellToProc_ = decomposePtr().decompose
 241             (
 242                 *this,
 243                 globalRegion,
 244                 regionCentres
 245             );
 246         }
 247     }
 248
 249     Info<< "\nFinished decomposition in "
 250         << decompositionTime.elapsedCpuTime()
 251         << " s" << endl;
 252 }
 253
 254
 255 // ************************************************************************* //