src/engine/engineTopoChangerMesh/twoStrokeEngine/twoStrokeEngineCheckMotionFluxes.C

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