src/postProcessing/functionObjects/check/mixingPlaneCheck/mixingPlaneCheckFunctionObject.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   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 Author
  25     Martin Beaudoin, Hydro-Quebec, 2009.  All rights reserved
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "mixingPlaneCheckFunctionObject.H"
  30 #include "addToRunTimeSelectionTable.H"
  31 #include "mixingPlaneFvsPatchFields.H"
  32 #include "surfaceFields.H"
  33 #include "mixingPlanePolyPatch.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39     defineTypeNameAndDebug(mixingPlaneCheckFunctionObject, 0);
  40
  41     addToRunTimeSelectionTable
  42     (
  43         functionObject,
  44         mixingPlaneCheckFunctionObject,
  45         dictionary
  46     );
  47 }
  48
  49
  50 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  51
  52 Foam::mixingPlaneCheckFunctionObject::mixingPlaneCheckFunctionObject
  53 (
  54     const word& name,
  55     const Time& t,
  56     const dictionary& dict
  57 )
  58 :
  59     functionObject(name),
  60     time_(t),
  61     regionName_(polyMesh::defaultRegion),
  62     phiName_(dict.lookup("phi"))
  63 {
  64     if (dict.found("region"))
  65     {
  66         dict.lookup("region") >> regionName_;
  67     }
  68
  69     Info << "Creating mixingPlane check functionObject" << endl;
  70 }
  71
  72
  73 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  74
  75 bool Foam::mixingPlaneCheckFunctionObject::start()
  76 {
  77     return true;
  78 }
  79
  80
  81 bool Foam::mixingPlaneCheckFunctionObject::execute()
  82 {
  83     const polyMesh& mesh =
  84         time_.lookupObject<polyMesh>(regionName_);
  85
  86     boolList visited(mesh.boundaryMesh().size(), false);
  87
  88     forAll (mesh.boundaryMesh(), patchI)
  89     {
  90         if (
  91             isA<mixingPlanePolyPatch>(mesh.boundaryMesh()[patchI]) &&
  92             mesh.boundaryMesh()[patchI].size()
  93         )
  94         {
  95             if (!visited[patchI])
  96             {
  97                 visited[patchI] = true;
  98
  99                 const mixingPlanePolyPatch& mixingMaster =
 100                     refCast<const mixingPlanePolyPatch>
 101                     (
 102                         mesh.boundaryMesh()[patchI]
 103                     );
 104
 105                 const label shadowPatchI = mixingMaster.shadowIndex();
 106
 107                 visited[shadowPatchI] = true;
 108
 109                 const mixingPlanePolyPatch& mixingShadow =
 110                     mixingMaster.shadow();
 111
 112                 // Get access to the mixing plane patch
 113                 const standAlonePatch& mixingPlanePatch =
 114                     mixingMaster.patchToPatch().mixingPlanePatch();
 115
 116                 // Calculate areas of the mixing patch
 117                 scalarField mixingPlanePatchAreas(mixingPlanePatch.size());
 118                 const vectorField& mixingPlanePatchPoints =
 119                     mixingPlanePatch.points();
 120
 121                 forAll (mixingPlanePatch, faceI)
 122                 {
 123                     mixingPlanePatchAreas[faceI] =
 124                         mixingPlanePatch[faceI].mag(mixingPlanePatchPoints);
 125                 }
 126
 127                 const scalarField masterAreas = mag(mixingMaster.faceAreas());
 128                 const scalarField shadowAreas = mag(mixingShadow.faceAreas());
 129
 130                 // Until the mixingPlane is fully parallelized, we stick with
 131                 // the serial version of sum. The interface is residing on a
 132                 // single processor when running in parallel
 133                 scalar sumMasterAreas = sum(masterAreas);
 134                 scalar sumShadowAreas = sum(shadowAreas);
 135 //                 scalar sumMixingAreas = sum(mixingPlanePatchAreas);
 136
 137 #if 0   // Remove this for now
 138                 Info<< "Mixing plane functionObject: area check " << nl
 139                     << "     "  << "Master " << mixingMaster.name()
 140                     << " = " << sumMasterAreas << nl
 141                     << "     "  << "Shadow  " << mixingMaster.shadow().name()
 142                     << " = " << sumShadowAreas << nl
 143                     << "     "  << "Mixing = " << sumMixingAreas << nl
 144                     << "     "  << "mixingPlanePatchAreas = " << mixingPlanePatchAreas
 145                     //<< "     "  << "mixingPlanePatchPoints = " << mixingPlanePatchPoints
 146                     << endl;
 147 #endif
 148
 149                 // Calculate master to strip sum
 150                 scalarField masterToStripsAreas(mixingPlanePatch.size(), 0);
 151
 152                 const labelListList& mppAddr =
 153                     mixingMaster.patchToPatch().masterProfileToPatchAddr();
 154
 155                 const scalarListList& mppWeights =
 156                     mixingMaster.patchToPatch().masterProfileToPatchWeights();
 157
 158                 forAll (masterAreas, masterI)
 159                 {
 160                     const labelList& curMppAddr = mppAddr[masterI];
 161                     const scalarList& curMppWeights = mppWeights[masterI];
 162
 163                     forAll (curMppAddr, i)
 164                     {
 165                         masterToStripsAreas[curMppAddr[i]] +=
 166                             curMppWeights[i]*masterAreas[masterI];
 167                     }
 168                 }
 169
 170 #if 0
 171                 // Does not work when in cylindrical coordinates
 172                 Info<< "Master scaling = "
 173                     << mixingPlanePatchAreas/masterToStripsAreas << endl;
 174 #endif
 175
 176                 // Calculate shadow to strip sum
 177                 scalarField shadowToStripsAreas(mixingPlanePatch.size(), 0);
 178
 179                 const labelListList& sppAddr =
 180                     mixingMaster.patchToPatch().slaveProfileToPatchAddr();
 181
 182                 const scalarListList& sppWeights =
 183                     mixingMaster.patchToPatch().slaveProfileToPatchWeights();
 184
 185                 forAll (shadowAreas, shadowI)
 186                 {
 187                     const labelList& curSppAddr = sppAddr[shadowI];
 188                     const scalarList& curSppWeights = sppWeights[shadowI];
 189
 190                     forAll (curSppAddr, i)
 191                     {
 192                         shadowToStripsAreas[curSppAddr[i]] +=
 193                             curSppWeights[i]*shadowAreas[shadowI];
 194                     }
 195                 }
 196
 197 #if 0
 198                 // Does not work when in cylindrical coordinates
 199                 Info<< "Shadow scaling = "
 200                     << mixingPlanePatchAreas/shadowToStripsAreas << endl;
 201 #endif
 202
 203                // Old way of computing phi balance
 204
 205                 if( mesh.foundObject<surfaceScalarField>(phiName_) )
 206                 {
 207                     const surfaceScalarField& phi =
 208                         mesh.lookupObject<surfaceScalarField>(phiName_);
 209
 210                     // Calculate local and shadow flux
 211                     scalar masterPatchScaleFactor_ = 1.0;
 212                     scalar shadowPatchScaleFactor_ = sumMasterAreas/sumShadowAreas;
 213                     scalar localFlux    = masterPatchScaleFactor_ * sum(phi.boundaryField()[patchI]);
 214                     scalar localFluxMag = mag(localFlux);
 215
 216                     scalar shadowFlux    = shadowPatchScaleFactor_ * sum(phi.boundaryField()[shadowPatchI]);
 217 //                     scalar shadowFluxMag = mag(shadowFlux);
 218
 219                     Info<< "Mixing plane pair "
 220                         << "(" << mixingMaster.name() << ", "
 221                         << mixingMaster.shadow().name() << ") : "
 222                         << localFlux << " " << shadowFlux
 223                         << " Diff = " << localFlux + shadowFlux << " or "
 224                         << mag(localFlux + shadowFlux)/(localFluxMag + SMALL)*100
 225                         << " %" << endl;
 226                 }
 227             }
 228         }
 229     }
 230
 231
 232 //         {
 233 //             word patchName = phi.boundaryField()[patchI].patch().name();
 234
 235 //             if (patchName == masterPatchName_ && !visited[patchI])
 236 //             {
 237 //                 visited[patchI] = true;
 238
 239 //                 // Calculate local and shadow flux
 240 //                 scalar localFlux    = masterPatchScaleFactor_ * sum(phi.boundaryField()[patchI]);
 241 //                 //scalar localFluxMag = masterPatchScaleFactor_ * sumMag(phi.boundaryField()[patchI]);
 242 //                 scalar localFluxMag = mag(localFlux);
 243
 244 //                 const mixingPlanePolyPatch& mixingPlanePatch =
 245 //                     refCast<const mixingPlanePolyPatch>
 246 //                     (
 247 //                         phi.boundaryField()[patchI].patch().patch()
 248 //                     );
 249
 250 //                 const label shadowPatchI = mixingPlanePatch.shadowIndex();
 251
 252 //                 visited[shadowPatchI] = true;
 253
 254 //                 scalar shadowFlux    = shadowPatchScaleFactor_ * sum(phi.boundaryField()[shadowPatchI]);
 255 //                 //scalar shadowFluxMag = shadowPatchScaleFactor_ * sumMag(phi.boundaryField()[shadowPatchI]);
 256 //                 scalar shadowFluxMag = mag(shadowFlux);
 257
 258 //                 Info<< "mixingPlane pair " << name_ << " (" << mixingPlanePatch.name() << ", " << mixingPlanePatch.shadow().name() << ") : "
 259 //                     << " flux: " << localFlux << " " << shadowFlux
 260 //                     << " : mag: " <<  localFluxMag << " " << shadowFluxMag
 261 //                     << " Diff = " << localFlux + shadowFlux << " or "
 262 //                     << mag(localFlux + shadowFlux)/(localFluxMag + SMALL)*100
 263 //                     << " %" << endl;
 264 //             }
 265 //         }
 266 //     }
 267
 268     return true;
 269 }
 270
 271
 272 bool Foam::mixingPlaneCheckFunctionObject::read(const dictionary& dict)
 273 {
 274     return false;
 275 }
 276
 277
 278 // ************************************************************************* //