src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
   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 \*---------------------------------------------------------------------------*/
  25
  26 #include "CellZoneInjection.H"
  27 #include "mathematicalConstants.H"
  28 #include "polyMeshTetDecomposition.H"
  29 #include "globalIndex.H"
  30 #include "Pstream.H"
  31
  32 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  33
  34 template<class CloudType>
  35 void Foam::CellZoneInjection<CloudType>::setPositions
  36 (
  37     const labelList& cellZoneCells
  38 )
  39 {
  40     const fvMesh& mesh = this->owner().mesh();
  41     const scalarField& V = mesh.V();
  42     const label nCells = cellZoneCells.size();
  43     cachedRandom& rnd = this->owner().rndGen();
  44
  45     DynamicList<vector> positions(nCells);          // initial size only
  46     DynamicList<label> injectorCells(nCells);       // initial size only
  47     DynamicList<label> injectorTetFaces(nCells);    // initial size only
  48     DynamicList<label> injectorTetPts(nCells);      // initial size only
  49
  50     scalar newParticlesTotal = 0.0;
  51     label addParticlesTotal = 0;
  52
  53     forAll(cellZoneCells, i)
  54     {
  55         const label cellI = cellZoneCells[i];
  56
  57         // Calc number of particles to add
  58         const scalar newParticles = V[cellI]*numberDensity_;
  59         newParticlesTotal += newParticles;
  60         label addParticles = floor(newParticles);
  61         addParticlesTotal += addParticles;
  62
  63         const scalar diff = newParticlesTotal - addParticlesTotal;
  64         if (diff > 1)
  65         {
  66             label corr = floor(diff);
  67             addParticles += corr;
  68             addParticlesTotal += corr;
  69         }
  70
  71         // Construct cell tet indices
  72         const List<tetIndices> cellTetIs =
  73             polyMeshTetDecomposition::cellTetIndices(mesh, cellI);
  74
  75         // Construct cell tet volume fractions
  76         scalarList cTetVFrac(cellTetIs.size(), 0.0);
  77         for (label tetI = 1; tetI < cellTetIs.size() - 1; tetI++)
  78         {
  79             cTetVFrac[tetI] =
  80                 (cTetVFrac[tetI-1] + cellTetIs[tetI].tet(mesh).mag())/V[cellI];
  81         }
  82         cTetVFrac.last() = 1.0;
  83
  84         // Set new particle position and cellId
  85         for (label pI = 0; pI < addParticles; pI++)
  86         {
  87             const scalar volFrac = rnd.sample01<scalar>();
  88             label tetI = 0;
  89             forAll(cTetVFrac, vfI)
  90             {
  91                 if (cTetVFrac[vfI] > volFrac)
  92                 {
  93                     tetI = vfI;
  94                     break;
  95                 }
  96             }
  97             positions.append(cellTetIs[tetI].tet(mesh).randomPoint(rnd));
  98
  99             injectorCells.append(cellI);
 100             injectorTetFaces.append(cellTetIs[tetI].face());
 101             injectorTetPts.append(cellTetIs[tetI].faceBasePt());
 102         }
 103     }
 104
 105     // Parallel operation manipulations
 106     globalIndex globalPositions(positions.size());
 107     List<vector> allPositions(globalPositions.size(), point::max);
 108     List<label> allInjectorCells(globalPositions.size(), -1);
 109     List<label> allInjectorTetFaces(globalPositions.size(), -1);
 110     List<label> allInjectorTetPts(globalPositions.size(), -1);
 111
 112     // Gather all positions on to all processors
 113     SubList<vector>
 114     (
 115         allPositions,
 116         globalPositions.localSize(Pstream::myProcNo()),
 117         globalPositions.offset(Pstream::myProcNo())
 118     ).assign(positions);
 119
 120     Pstream::listCombineGather(allPositions, minEqOp<point>());
 121     Pstream::listCombineScatter(allPositions);
 122
 123     // Gather local cell tet and tet-point Ids, but leave non-local ids set -1
 124     SubList<label>
 125     (
 126         allInjectorCells,
 127         globalPositions.localSize(Pstream::myProcNo()),
 128         globalPositions.offset(Pstream::myProcNo())
 129     ).assign(injectorCells);
 130     SubList<label>
 131     (
 132         allInjectorTetFaces,
 133         globalPositions.localSize(Pstream::myProcNo()),
 134         globalPositions.offset(Pstream::myProcNo())
 135     ).assign(injectorTetFaces);
 136     SubList<label>
 137     (
 138         allInjectorTetPts,
 139         globalPositions.localSize(Pstream::myProcNo()),
 140         globalPositions.offset(Pstream::myProcNo())
 141     ).assign(injectorTetPts);
 142
 143     // Transfer data
 144     positions_.transfer(allPositions);
 145     injectorCells_.transfer(allInjectorCells);
 146     injectorTetFaces_.transfer(allInjectorTetFaces);
 147     injectorTetPts_.transfer(allInjectorTetPts);
 148
 149
 150     if (debug)
 151     {
 152         OFstream points("points.obj");
 153         forAll(positions_, i)
 154         {
 155             meshTools::writeOBJ(points, positions_[i]);
 156         }
 157     }
 158 }
 159
 160
 161 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 162
 163 template<class CloudType>
 164 Foam::CellZoneInjection<CloudType>::CellZoneInjection
 165 (
 166     const dictionary& dict,
 167     CloudType& owner
 168 )
 169 :
 170     InjectionModel<CloudType>(dict, owner, typeName),
 171     cellZoneName_(this->coeffDict().lookup("cellZone")),
 172     numberDensity_(readScalar(this->coeffDict().lookup("numberDensity"))),
 173     positions_(),
 174     injectorCells_(),
 175     injectorTetFaces_(),
 176     injectorTetPts_(),
 177     diameters_(),
 178     U0_(this->coeffDict().lookup("U0")),
 179     sizeDistribution_
 180     (
 181         distributionModels::distributionModel::New
 182         (
 183             this->coeffDict().subDict("sizeDistribution"), owner.rndGen()
 184         )
 185     )
 186 {
 187     const fvMesh& mesh = owner.mesh();
 188     const label zoneI = mesh.cellZones().findZoneID(cellZoneName_);
 189
 190     if (zoneI < 0)
 191     {
 192         FatalErrorIn
 193         (
 194             "Foam::CellZoneInjection<CloudType>::CellZoneInjection"
 195             "("
 196                 "const dictionary&, "
 197                 "CloudType&"
 198             ")"
 199         )   << "Unknown cell zone name: " << cellZoneName_
 200             << ". Valid cell zones are: " << mesh.cellZones().names()
 201             << nl << exit(FatalError);
 202     }
 203
 204     const labelList& cellZoneCells = mesh.cellZones()[zoneI];
 205     const label nCells = cellZoneCells.size();
 206     const scalar nCellsTotal = returnReduce(nCells, sumOp<label>());
 207     const scalar VCells = sum(scalarField(mesh.V(), cellZoneCells));
 208     const scalar VCellsTotal = returnReduce(VCells, sumOp<scalar>());
 209     Info<< "    cell zone size      = " << nCellsTotal << endl;
 210     Info<< "    cell zone volume    = " << VCellsTotal << endl;
 211
 212     if ((nCells == 0) || (VCellsTotal*numberDensity_ < 1))
 213     {
 214         WarningIn
 215         (
 216             "Foam::CellZoneInjection<CloudType>::CellZoneInjection"
 217             "("
 218                 "const dictionary&, "
 219                 "CloudType&"
 220             ")"
 221         )   << "Number of particles to be added to cellZone " << cellZoneName_
 222             << " is zero" << endl;
 223     }
 224     else
 225     {
 226         setPositions(cellZoneCells);
 227
 228         Info<< "    number density      = " << numberDensity_ << nl
 229             << "    number of particles = " << positions_.size() << endl;
 230
 231         // Construct parcel diameters
 232         diameters_.setSize(positions_.size());
 233         forAll(diameters_, i)
 234         {
 235             diameters_[i] = sizeDistribution_->sample();
 236         }
 237     }
 238
 239     // Determine volume of particles to inject
 240     this->volumeTotal_ = sum(pow3(diameters_))*constant::mathematical::pi/6.0;
 241 }
 242
 243
 244 template<class CloudType>
 245 Foam::CellZoneInjection<CloudType>::CellZoneInjection
 246 (
 247     const CellZoneInjection<CloudType>& im
 248 )
 249 :
 250     InjectionModel<CloudType>(im),
 251     cellZoneName_(im.cellZoneName_),
 252     numberDensity_(im.numberDensity_),
 253     positions_(im.positions_),
 254     injectorCells_(im.injectorCells_),
 255     injectorTetFaces_(im.injectorTetFaces_),
 256     injectorTetPts_(im.injectorTetPts_),
 257     diameters_(im.diameters_),
 258     U0_(im.U0_),
 259     sizeDistribution_(im.sizeDistribution_().clone().ptr())
 260 {}
 261
 262
 263 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 264
 265 template<class CloudType>
 266 Foam::CellZoneInjection<CloudType>::~CellZoneInjection()
 267 {}
 268
 269
 270 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 271
 272 template<class CloudType>
 273 Foam::scalar Foam::CellZoneInjection<CloudType>::timeEnd() const
 274 {
 275     // Not used
 276     return this->SOI_;
 277 }
 278
 279
 280 template<class CloudType>
 281 Foam::label Foam::CellZoneInjection<CloudType>::parcelsToInject
 282 (
 283     const scalar time0,
 284     const scalar time1
 285 )
 286 {
 287     if ((0.0 >= time0) && (0.0 < time1))
 288     {
 289         return positions_.size();
 290     }
 291     else
 292     {
 293         return 0;
 294     }
 295 }
 296
 297
 298 template<class CloudType>
 299 Foam::scalar Foam::CellZoneInjection<CloudType>::volumeToInject
 300 (
 301     const scalar time0,
 302     const scalar time1
 303 )
 304 {
 305     // All parcels introduced at SOI
 306     if ((0.0 >= time0) && (0.0 < time1))
 307     {
 308         return this->volumeTotal_;
 309     }
 310     else
 311     {
 312         return 0.0;
 313     }
 314 }
 315
 316
 317 template<class CloudType>
 318 void Foam::CellZoneInjection<CloudType>::setPositionAndCell
 319 (
 320     const label parcelI,
 321     const label nParcels,
 322     const scalar time,
 323     vector& position,
 324     label& cellOwner,
 325     label& tetFaceI,
 326     label& tetPtI
 327 )
 328 {
 329     position = positions_[parcelI];
 330     cellOwner = injectorCells_[parcelI];
 331     tetFaceI = injectorTetFaces_[parcelI];
 332     tetPtI = injectorTetPts_[parcelI];
 333 }
 334
 335
 336 template<class CloudType>
 337 void Foam::CellZoneInjection<CloudType>::setProperties
 338 (
 339     const label parcelI,
 340     const label,
 341     const scalar,
 342     typename CloudType::parcelType& parcel
 343 )
 344 {
 345     // set particle velocity
 346     parcel.U() = U0_;
 347
 348     // set particle diameter
 349     parcel.d() = diameters_[parcelI];
 350 }
 351
 352
 353 template<class CloudType>
 354 bool Foam::CellZoneInjection<CloudType>::fullyDescribed() const
 355 {
 356     return false;
 357 }
 358
 359
 360 template<class CloudType>
 361 bool Foam::CellZoneInjection<CloudType>::validInjection(const label)
 362 {
 363     return true;
 364 }
 365
 366
 367 // ************************************************************************* //