src/engine/engineTopoChangerMesh/twoStrokeEngine/twoStrokeEngineCheckMotionFluxes.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   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 \*---------------------------------------------------------------------------*/
  25
  26
  27 #include "twoStrokeEngine.H"
  28 #include "slidingInterface.H"
  29 #include "layerAdditionRemoval.H"
  30 #include "surfaceFields.H"
  31 #include "regionSplit.H"
  32 #include "componentMixedTetPolyPatchVectorField.H"
  33 #include "mapPolyMesh.H"
  34
  35 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  36 void Foam::twoStrokeEngine::checkMotionFluxes()
  37 {
  38 /*
  39
  40     //checking the motion fluxes for the cutFaceZone
  41
  42     Info << "Checking the motion fluxes in the cut-face zone" << endl;
  43
  44     const labelList& cutFaceZoneAddressing =
  45         faceZones()[faceZones().findZoneID("cutFaceZone")];
  46
  47     boolList calculatedMeshPhi(V().size(), false);
  48     scalarField sumMeshPhi(V().size(), 0.0);
  49
  50     forAll(cutFaceZoneAddressing, i)
  51     {
  52         label facei = cutFaceZoneAddressing[i];
  53
  54         calculatedMeshPhi[owner()[facei]] = true;
  55         calculatedMeshPhi[neighbour()[facei]] = true;
  56     }
  57
  58     Info << "Checking the motion fluxes in the liner zone" << endl;
  59
  60     const labelList& linerAddressing =
  61         faceZones()[faceZones().findZoneID(scavInCylPatchName_ + "Zone")];
  62
  63     forAll(linerAddressing, i)
  64     {
  65         label facei = linerAddressing[i];
  66
  67         calculatedMeshPhi[owner()[facei]] = true;
  68         calculatedMeshPhi[neighbour()[facei]] = true;
  69     }
  70
  71     const polyPatch& wallPatch =
  72         boundaryMesh()[boundaryMesh().findPatchID("wall")];
  73
  74     labelList wallLabels(wallPatch.size());
  75
  76     forAll (wallLabels, i)
  77     {
  78         wallLabels[i] = wallPatch.start() + i;
  79     }
  80
  81     forAll(wallLabels, i)
  82     {
  83         label facei = wallLabels[i];
  84
  85         calculatedMeshPhi[owner()[facei]] = true;
  86         calculatedMeshPhi[neighbour()[facei]] = true;
  87     }
  88
  89     forAll(owner(), facei)
  90     {
  91         sumMeshPhi[owner()[facei]] += phi()[facei];
  92         sumMeshPhi[neighbour()[facei]] -= phi()[facei];
  93     }
  94
  95     forAll(boundary(), patchi)
  96     {
  97         const unallocLabelList& pFaceCells =
  98             boundary()[patchi].faceCells();
  99
 100         const fvsPatchField<scalar>& pssf = phi().boundaryField()[patchi];
 101
 102         forAll(boundary()[patchi], facei)
 103         {
 104             sumMeshPhi[pFaceCells[facei]] += pssf[facei];
 105         }
 106     }
 107
 108     sumMeshPhi /= V();
 109
 110     DynamicList<label> checkMeshPhi;
 111     label checkMeshPhiSize = 0;
 112
 113     forAll(calculatedMeshPhi, i)
 114     {
 115         if(calculatedMeshPhi[i] == true)
 116         {
 117             checkMeshPhi.append(i);
 118             checkMeshPhiSize++;
 119         }
 120     }
 121
 122     Info << "checkMeshPhiSize = " << checkMeshPhiSize << endl;
 123
 124     checkMeshPhi.setSize(checkMeshPhiSize);
 125
 126     scalarField dVdt(checkMeshPhiSize, 0.0);
 127
 128     forAll(dVdt, i)
 129     {
 130         label celli = checkMeshPhi[i];
 131
 132         dVdt[i] = (1.0 - V0()[celli]/V()[celli])/engTime().deltaT().value();
 133     }
 134
 135     Info << "checking mesh flux in the cutFaces" << endl;
 136
 137     label nOutCheck = 0;
 138
 139     forAll(checkMeshPhi, i)
 140     {
 141         scalar sumCheck = dVdt[i] - sumMeshPhi[checkMeshPhi[i]];
 142         if(mag(sumCheck) > 1)
 143         {
 144             Info<< "LOCAL sumCheck[" << checkMeshPhi[i] <<  "] = "
 145                 << sumCheck << endl;
 146             nOutCheck++;
 147         }
 148     }
 149
 150     Info<< "found " << nOutCheck << " cells with inconsistent motion fluxes"
 151         << endl;
 152     Info << "end " <<  endl;
 153
 154 */
 155     {
 156
 157         scalarField dVdt = (1.0 - V0()/V())/engTime().deltaT();
 158
 159         scalarField sumMeshPhi(V().size(), 0.0);
 160
 161         forAll(owner(), facei)
 162         {
 163             sumMeshPhi[owner()[facei]] += phi()[facei];
 164             sumMeshPhi[neighbour()[facei]] -= phi()[facei];
 165         }
 166
 167         forAll(boundary(), patchi)
 168         {
 169             const unallocLabelList& pFaceCells =
 170                 boundary()[patchi].faceCells();
 171
 172             const fvsPatchField<scalar>& pssf = phi().boundaryField()[patchi];
 173
 174             forAll(boundary()[patchi], facei)
 175             {
 176                 sumMeshPhi[pFaceCells[facei]] += pssf[facei];
 177             }
 178         }
 179
 180         sumMeshPhi /= V();
 181
 182         scalarField sumCheck = dVdt - sumMeshPhi;
 183
 184         label nOutCheck = 0;
 185
 186         forAll(sumCheck, celli)
 187         {
 188             if(mag(sumCheck[celli]) > 1)
 189             {
 190                 Info<< "sumCheck GLOBAL[" << celli << "] = "
 191                     << sumCheck[celli] << endl;
 192                 nOutCheck++;
 193             }
 194         }
 195
 196         Info<< "found " << nOutCheck
 197             << " cells with inconsistent motion fluxes GLOBAL" << endl;
 198
 199     }
 200 }