src/engine/engineTopoChangerMesh/twoStrokeEngine/checkMotionFluxes.H

   1 {
   2     const fvsPatchScalarField& pistonMeshPhi =
   3         phi().boundaryField()[boundaryMesh().findPatchID("piston")];
   4
   5     Info << "max piston mesh phi 2 = " << max(mag(pistonMeshPhi)) << endl;
   6
   7     const labelList& cutFaces =
   8         faceZones()[faceZones().findZoneID("cutFaceZone")];
   9
  10     forAll(cutFaces, facei)
  11     {
  12         correctedMeshPhi[cutFaces[facei]] = 0.0;
  13     }
  14
  15 //     setPhi(correctedMeshPhi);
  16
  17     Info<< "max mesh phi, corrected cut faces = "
  18         << max(mag(phi().internalField())) << nl
  19         << "max mesh phi, corrected cut faces = "
  20         << min(mag(phi().internalField())) << endl;
  21
  22     const polyPatch& wallPatch =
  23         boundaryMesh()[boundaryMesh().findPatchID("wall")];
  24
  25     label wallSize = wallPatch.size();
  26
  27     fvsPatchScalarField& wallMeshPhi =
  28         correctedMeshPhi.boundaryField()
  29         [
  30             boundaryMesh().findPatchID("wall")
  31         ];
  32
  33     forAll(wallMeshPhi, k)
  34     {
  35         if(mag(wallMeshPhi[k]) > SMALL)
  36         {
  37             wallMeshPhi[k] = 0;
  38         }
  39     }
  40
  41     fvsPatchScalarField& linerMeshPhi =
  42         correctedMeshPhi.boundaryField()[boundaryMesh().findPatchID("liner")];
  43
  44     forAll(linerMeshPhi, k)
  45     {
  46         if(mag(linerMeshPhi[k]) > SMALL)
  47         {
  48             linerMeshPhi[k] = 0;
  49         }
  50     }
  51
  52 //        setPhi(correctedMeshPhi);
  53
  54     Info<< "max mesh phi after correction = " << max(mag(phi())) << nl
  55         << "min mesh phi after correction = " << min(mag(phi())) << endl;
  56
  57 //    if(1 == 0)
  58     {
  59         scalarField dVdt = (1.0 - V0()/V())/engTime().deltaT();
  60
  61         scalarField sumMeshPhi(V().size(), 0.0);
  62         scalarField sumMeshPhiCorrected(V().size(), 0.0);
  63
  64         forAll(owner(), facei)
  65         {
  66             sumMeshPhi[owner()[facei]] += phi()[facei];
  67             sumMeshPhi[neighbour()[facei]] -= phi()[facei];
  68             sumMeshPhiCorrected[owner()[facei]] += correctedMeshPhi[facei];
  69             sumMeshPhiCorrected[neighbour()[facei]] -= correctedMeshPhi[facei];
  70         }
  71
  72         forAll(boundary(), patchi)
  73         {
  74             const unallocLabelList& pFaceCells =
  75                 boundary()[patchi].faceCells();
  76
  77             const fvsPatchField<scalar>& pssf = phi().boundaryField()[patchi];
  78             const fvsPatchField<scalar>& pssfCorrected =
  79                 correctedMeshPhi.boundaryField()[patchi];
  80
  81             forAll(boundary()[patchi], facei)
  82             {
  83                 sumMeshPhi[pFaceCells[facei]] += pssf[facei];
  84                 sumMeshPhiCorrected[pFaceCells[facei]] += pssfCorrected[facei];
  85             }
  86         }
  87
  88         sumMeshPhi /= V();
  89
  90         sumMeshPhiCorrected /= V();
  91
  92         scalarField sumCheck = dVdt - sumMeshPhi;
  93         scalarField sumCheckCorrected = dVdt - sumMeshPhiCorrected;
  94
  95         label nOutCheck = 0;
  96         label nOutCheckCorrected = 0;
  97
  98         forAll(sumCheck, celli)
  99         {
 100             if(mag(sumCheck[celli]) > 1)
 101             {
 102                 nOutCheck++;
 103             }
 104         }
 105
 106         forAll(sumCheckCorrected, celli)
 107         {
 108             if (mag(sumCheckCorrected[celli]) > 1)
 109             {
 110                 nOutCheckCorrected++;
 111             }
 112         }
 113
 114         Info<< "found " << nOutCheck
 115             << " cells with inconsistent motion fluxes real" << endl;
 116         Info<< "found " << nOutCheckCorrected
 117             << " cells with inconsistent motion fluxes CORRECTED" << endl;
 118     }
 119
 120     /*
 121         newPoints = allPoints();
 122
 123 //        movePoints(newPoints);
 124
 125 //        pointField mappedOldPointsNew(newPoints.size());
 126         pointField mappedOldPointsNew(points().size());
 127
 128         mappedOldPointsNew.map(oldPointsNew, topoChangeMap4->pointMap());
 129
 130         // Solve the correct mesh motion to make sure motion fluxes
 131         // are solved for and not mapped
 132         movePoints(mappedOldPointsNew);
 133
 134 //        resetMotion();
 135
 136         Info << "max mesh phi after reset motion = " << max(mag(phi())) << endl;
 137         Info << "max mesh phi after reset motion = " << min(mag(phi())) << endl;
 138
 139         setV0();
 140         movePoints(newPoints);
 141 //        pointField newPointsNew = allPoints();
 142 //        movePoints(newPointsNew);
 143
 144         Info << "max mesh phi after corr = " << max(mag(phi())) << endl;
 145         Info << "max mesh phi after corr = " << min(mag(phi())) << endl;
 146
 147         */
 148
 149         // Changing topology by hand
 150 //        autoPtr<mapPolyMesh> topoChangeMap3 = topoChanger_.changeMesh();
 151
 152     Info << "Sliding interfaces coupled: " << attached() << endl;
 153     Info << "max mesh phi after correction = " << max(mag(phi())) << endl;
 154 }