applications/test/nearWallDist-wave/Test-WallDist2.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 Description
  25     Wave propagation of nearwall distance through grid. Every iteration
  26     information goes through one layer of cells.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "fvCFD.H"
  31 #include "wallFvPatch.H"
  32 #include "FaceCellWave.H"
  33 #include "wallPoint.H"
  34
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37 // Main program:
  38
  39 int main(int argc, char *argv[])
  40 {
  41 #   include "setRootCase.H"
  42 #   include "createTime.H"
  43 #   include "createMesh.H"
  44
  45     Info<< "Mesh read in = "
  46         << runTime.cpuTimeIncrement()
  47         << " s\n" << endl << endl;
  48
  49     Info<< "Creating field wDistNC\n" << endl;
  50     volScalarField wallDistUncorrected
  51     (
  52         IOobject
  53         (
  54             "wDistNC",
  55             runTime.timeName(),
  56             mesh,
  57             IOobject::NO_READ,
  58             IOobject::AUTO_WRITE
  59         ),
  60         mesh,
  61         dimensionedScalar
  62         (
  63             "wallDist",
  64             dimensionSet(0, 1, 0, 0, 0),
  65             0.0
  66         )
  67     );
  68
  69     //
  70     // Set initial changed faces: set wallPoint for wall faces to wall centre
  71     //
  72
  73     // Count walls
  74     label nWalls = 0;
  75     forAll(mesh.boundary(), patchI)
  76     {
  77         const fvPatch& patch = mesh.boundary()[patchI];
  78
  79         if (isA<wallFvPatch>(patch))
  80         {
  81             nWalls += patch.size();
  82         }
  83     }
  84
  85     List<wallPoint> faceDist(nWalls);
  86     labelList changedFaces(nWalls);
  87
  88     label nChangedFaces = 0;
  89     forAll(mesh.boundary(), patchI)
  90     {
  91         const fvPatch& patch = mesh.boundary()[patchI];
  92
  93         if (isA<wallFvPatch>(patch))
  94         {
  95             forAll(patch.Cf(), patchFaceI)
  96             {
  97                 const polyPatch& polyPatch = mesh.boundaryMesh()[patchI];
  98
  99                 label meshFaceI = polyPatch.start() + patchFaceI;
 100
 101                 changedFaces[nChangedFaces] = meshFaceI;
 102
 103                 faceDist[nChangedFaces] =
 104                     wallPoint(patch.Cf()[patchFaceI], 0.0);
 105
 106                 nChangedFaces++;
 107             }
 108         }
 109     }
 110
 111     List<wallPoint> allFaceInfo(mesh.nFaces());
 112     List<wallPoint> allCellInfo(mesh.nCells());
 113
 114     FaceCellWave<wallPoint> wallDistCalc
 115     (
 116         mesh,
 117         changedFaces,
 118         faceDist,
 119         allFaceInfo,
 120         allCellInfo,
 121         0             // max iterations
 122     );
 123
 124     Info<< "\nStarting time loop\n" << endl;
 125
 126     while (runTime.loop())
 127     {
 128         Info<< "Time = " << runTime.timeName() << endl;
 129
 130
 131         label nCells = wallDistCalc.faceToCell();
 132
 133         Info<< "    Total changed cells   : " << nCells << endl;
 134
 135         if (nCells == 0)
 136         {
 137             break;
 138         }
 139
 140
 141         label nFaces = wallDistCalc.cellToFace();
 142
 143         Info<< "    Total changed faces   : " << nFaces << endl;
 144
 145         if (nFaces == 0)
 146         {
 147             break;
 148         }
 149
 150
 151         //
 152         // Copy face and cell values into field
 153         //
 154
 155         label nIllegal = 0;
 156
 157         // Copy cell values
 158         forAll(allCellInfo, cellI)
 159         {
 160             scalar dist = allCellInfo[cellI].distSqr();
 161             if (allCellInfo[cellI].valid(wallDistCalc.data()))
 162             {
 163                 wallDistUncorrected[cellI] = Foam::sqrt(dist);
 164             }
 165             else
 166             {
 167                 wallDistUncorrected[cellI] = -1;
 168                 nIllegal++;
 169             }
 170         }
 171
 172         // Copy boundary values
 173         forAll(wallDistUncorrected.boundaryField(), patchI)
 174         {
 175             fvPatchScalarField& patchField =
 176                 wallDistUncorrected.boundaryField()[patchI];
 177
 178             forAll(patchField, patchFaceI)
 179             {
 180                 const label meshFaceI = patchField.patch().start() + patchFaceI;
 181
 182                 scalar dist = allFaceInfo[meshFaceI].distSqr();
 183                 if (allFaceInfo[meshFaceI].valid(wallDistCalc.data()))
 184                 {
 185                     patchField[patchFaceI] = Foam::sqrt(dist);
 186                 }
 187                 else
 188                 {
 189                     patchField[patchFaceI] = dist;
 190                     nIllegal++;
 191                 }
 192             }
 193         }
 194
 195         Info<< "nIllegal:" << nIllegal << endl;
 196
 197
 198         //
 199         // Write it
 200         //
 201
 202         wallDistUncorrected.write();
 203
 204         Info<< "ExecutionTime = "
 205             << runTime.elapsedCpuTime()
 206             << " s\n" << endl << endl;
 207     }
 208
 209     Info<< "End\n" << endl;
 210
 211     return 0;
 212 }
 213
 214
 215 // ************************************************************************* //