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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "InflationInjection.H"
  27 #include "mathematicalConstants.H"
  28 #include "PackedBoolList.H"
  29 #include "cellSet.H"
  30 #include "ListListOps.H"
  31
  32 using namespace Foam::constant::mathematical;
  33
  34 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  35
  36 template<class CloudType>
  37 Foam::InflationInjection<CloudType>::InflationInjection
  38 (
  39     const dictionary& dict,
  40     CloudType& owner
  41 )
  42 :
  43     InjectionModel<CloudType>(dict, owner, typeName),
  44     generationSetName_(this->coeffDict().lookup("generationCellSet")),
  45     inflationSetName_(this->coeffDict().lookup("inflationCellSet")),
  46     generationCells_(),
  47     inflationCells_(),
  48     duration_(readScalar(this->coeffDict().lookup("duration"))),
  49     flowRateProfile_
  50     (
  51         DataEntry<scalar>::New
  52         (
  53             "flowRateProfile",
  54             this->coeffDict()
  55         )
  56     ),
  57     growthRate_
  58     (
  59         DataEntry<scalar>::New
  60         (
  61             "growthRate",
  62             this->coeffDict()
  63         )
  64     ),
  65     newParticles_(),
  66     volumeAccumulator_(0.0),
  67     fraction_(1.0),
  68     selfSeed_(this->coeffDict().lookupOrDefault("selfSeed", false)),
  69     dSeed_(SMALL),
  70     sizeDistribution_
  71     (
  72         distributionModels::distributionModel::New
  73         (
  74             this->coeffDict().subDict("sizeDistribution"),
  75             owner.rndGen()
  76         )
  77     )
  78 {
  79     if (selfSeed_)
  80     {
  81         dSeed_ = readScalar(this->coeffDict().lookup("dSeed"));
  82     }
  83
  84     cellSet generationCells(this->owner().mesh(), generationSetName_);
  85
  86     generationCells_ = generationCells.toc();
  87
  88     cellSet inflationCells(this->owner().mesh(), inflationSetName_);
  89
  90     // Union of cellSets
  91     inflationCells |= generationCells;
  92
  93     inflationCells_ = inflationCells.toc();
  94
  95     if (Pstream::parRun())
  96     {
  97         scalar generationVolume = 0.0;
  98
  99         forAll(generationCells_, gCI)
 100         {
 101             label cI = generationCells_[gCI];
 102
 103             generationVolume += this->owner().mesh().cellVolumes()[cI];
 104         }
 105
 106         scalar totalGenerationVolume = generationVolume;
 107
 108         reduce(totalGenerationVolume, sumOp<scalar>());
 109
 110         fraction_ = generationVolume/totalGenerationVolume;
 111     }
 112
 113     // Set total volume/mass to inject
 114     this->volumeTotal_ = fraction_*flowRateProfile_().integrate(0.0, duration_);
 115     this->massTotal_ *= fraction_;
 116 }
 117
 118
 119 template<class CloudType>
 120 Foam::InflationInjection<CloudType>::InflationInjection
 121 (
 122     const Foam::InflationInjection<CloudType>& im
 123 )
 124 :
 125     InjectionModel<CloudType>(im),
 126     generationSetName_(im.generationSetName_),
 127     inflationSetName_(im.inflationSetName_),
 128     generationCells_(im.generationCells_),
 129     inflationCells_(im.inflationCells_),
 130     duration_(im.duration_),
 131     flowRateProfile_(im.flowRateProfile_().clone().ptr()),
 132     growthRate_(im.growthRate_().clone().ptr()),
 133     newParticles_(im.newParticles_),
 134     volumeAccumulator_(im.volumeAccumulator_),
 135     fraction_(im.fraction_),
 136     selfSeed_(im.selfSeed_),
 137     dSeed_(im.dSeed_),
 138     sizeDistribution_(im.sizeDistribution_().clone().ptr())
 139 {}
 140
 141
 142 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 143
 144 template<class CloudType>
 145 Foam::InflationInjection<CloudType>::~InflationInjection()
 146 {}
 147
 148
 149 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 150
 151 template<class CloudType>
 152 Foam::scalar Foam::InflationInjection<CloudType>::timeEnd() const
 153 {
 154     return this->SOI_ + duration_;
 155 }
 156
 157
 158 template<class CloudType>
 159 Foam::label Foam::InflationInjection<CloudType>::parcelsToInject
 160 (
 161     const scalar time0,
 162     const scalar time1
 163 )
 164 {
 165     const polyMesh& mesh = this->owner().mesh();
 166
 167     List<DynamicList<typename CloudType::parcelType*> >& cellOccupancy =
 168         this->owner().cellOccupancy();
 169
 170     scalar gR = growthRate_().value(time1);
 171
 172     scalar dT = time1 - time0;
 173
 174     // Inflate existing particles
 175
 176     forAll(inflationCells_, iCI)
 177     {
 178         label cI = inflationCells_[iCI];
 179
 180         typename CloudType::parcelType* pPtr = NULL;
 181
 182         forAll(cellOccupancy[cI], cPI)
 183         {
 184             pPtr = cellOccupancy[cI][cPI];
 185
 186             scalar dTarget = pPtr->dTarget();
 187
 188             pPtr->d() = min(dTarget, pPtr->d() + gR*dT);
 189         }
 190     }
 191
 192     // Generate new particles
 193
 194     newParticles_.clear();
 195
 196     cachedRandom& rnd = this->owner().rndGen();
 197
 198     // Diameter factor, when splitting particles into 4, this is the
 199     // factor that modifies the diameter.
 200     scalar dFact = sqrt(2.0)/(sqrt(3.0) + sqrt(2.0));
 201
 202     if ((time0 >= 0.0) && (time0 < duration_))
 203     {
 204          volumeAccumulator_ +=
 205              fraction_*flowRateProfile_().integrate(time0, time1);
 206     }
 207
 208     labelHashSet cellCentresUsed;
 209
 210     // Loop escape counter
 211     label maxIterations = max
 212     (
 213         1,
 214         (10*volumeAccumulator_)
 215        /CloudType::parcelType::volume(sizeDistribution_().minValue())
 216     );
 217
 218     label iterationNo = 0;
 219
 220     // Info<< "Accumulated volume to inject: "
 221     //     << returnReduce(volumeAccumulator_, sumOp<scalar>()) << endl;
 222
 223     while (!generationCells_.empty() && volumeAccumulator_ > 0)
 224     {
 225         if (iterationNo > maxIterations)
 226         {
 227             WarningIn
 228             (
 229                 "Foam::label "
 230                 "Foam::InflationInjection<CloudType>::parcelsToInject"
 231                 "("
 232                     "const scalar time0, "
 233                     "const scalar time1"
 234                 ")"
 235             )
 236                 << "Maximum particle split iterations ("
 237                 << maxIterations << ") exceeded" << endl;
 238
 239             break;
 240         }
 241
 242         label cI =
 243             generationCells_[rnd.position(0, generationCells_.size() - 1)];
 244
 245         // Pick a particle at random from the cell - if there are
 246         // none, insert one at the cell centre.  Otherwise, split an
 247         // existing particle into four new ones.
 248
 249         if (cellOccupancy[cI].empty())
 250         {
 251             if (selfSeed_ && !cellCentresUsed.found(cI))
 252             {
 253                 scalar dNew = sizeDistribution_().sample();
 254
 255                 newParticles_.append
 256                 (
 257                     vectorPairScalarPair
 258                     (
 259                         Pair<vector>(mesh.cellCentres()[cI], vector::zero),
 260                         Pair<scalar>(dSeed_, dNew)
 261                     )
 262                 );
 263
 264                 volumeAccumulator_ -= CloudType::parcelType::volume(dNew);
 265
 266                 cellCentresUsed.insert(cI);
 267             }
 268         }
 269         else
 270         {
 271             label cPI = rnd.position(0, cellOccupancy[cI].size() - 1);
 272
 273             // This has to be a reference to the pointer so that it
 274             // can be set to NULL when the particle is deleted.
 275             typename CloudType::parcelType*& pPtr = cellOccupancy[cI][cPI];
 276
 277             if (pPtr != NULL)
 278             {
 279                 scalar pD = pPtr->d();
 280
 281                 // Select bigger particles by preference
 282                 if ((pD/pPtr->dTarget()) < rnd.sample01<scalar>())
 283                 {
 284                     continue;
 285                 }
 286
 287                 const point& pP = pPtr->position();
 288                 const vector& pU = pPtr->U();
 289
 290                 // Generate a tetrahedron of new positions with the
 291                 // four new spheres fitting inside the old one, where
 292                 // a is the diameter of the new spheres, and is
 293                 // related to the diameter of the enclosing sphere, A,
 294                 // by a = sqrt(2)*A/(sqrt(3) + sqrt(2));
 295
 296                 // Positions around the origin, which is the
 297                 // tetrahedron centroid (centre of old sphere).
 298
 299                 // x = a/sqrt(3)
 300                 // r = a/(2*sqrt(6))
 301                 // R = sqrt(3)*a/(2*sqrt(2))
 302                 // d = a/(2*sqrt(3))
 303
 304                 // p0(x, 0, -r)
 305                 // p1(-d, a/2, -r)
 306                 // p2(-d, -a/2, -r)
 307                 // p3(0, 0, R)
 308
 309                 scalar a = pD*dFact;
 310
 311                 scalar x = a/sqrt(3.0);
 312                 scalar r = a/(2.0*sqrt(6.0));
 313                 scalar R = sqrt(3.0)*a/(2.0*sqrt(2.0));
 314                 scalar d = a/(2.0*sqrt(3.0));
 315
 316                 scalar dNew = sizeDistribution_().sample();
 317                 scalar volNew = CloudType::parcelType::volume(dNew);
 318
 319                 newParticles_.append
 320                 (
 321                     vectorPairScalarPair
 322                     (
 323                         Pair<vector>(vector(x, 0, -r)  + pP, pU),
 324                         Pair<scalar>(a, dNew)
 325                     )
 326                 );
 327                 volumeAccumulator_ -= volNew;
 328
 329                 dNew = sizeDistribution_().sample();
 330                 newParticles_.append
 331                 (
 332                     vectorPairScalarPair
 333                     (
 334                         Pair<vector>(vector(-d, a/2, -r) + pP, pU),
 335                         Pair<scalar>(a, dNew)
 336                     )
 337                 );
 338                 volumeAccumulator_ -= volNew;
 339
 340                 dNew = sizeDistribution_().sample();
 341                 newParticles_.append
 342                 (
 343                     vectorPairScalarPair
 344                     (
 345                         Pair<vector>(vector(-d, -a/2, -r) + pP, pU),
 346                         Pair<scalar>(a, dNew)
 347                     )
 348                 );
 349                 volumeAccumulator_ -= volNew;
 350
 351                 dNew = sizeDistribution_().sample();
 352                 newParticles_.append
 353                 (
 354                     vectorPairScalarPair
 355                     (
 356                         Pair<vector>(vector(0, 0, R) + pP, pU),
 357                         Pair<scalar>(a, dNew)
 358                     )
 359                 );
 360                 volumeAccumulator_ -= volNew;
 361
 362                 // Account for the lost volume of the particle which
 363                 // is to be deleted
 364                 volumeAccumulator_ += CloudType::parcelType::volume
 365                 (
 366                     pPtr->dTarget()
 367                 );
 368
 369                 this->owner().deleteParticle(*pPtr);
 370
 371                 pPtr = NULL;
 372             }
 373         }
 374
 375         iterationNo++;
 376     }
 377
 378     if (Pstream::parRun())
 379     {
 380         List<List<vectorPairScalarPair> > gatheredNewParticles
 381         (
 382             Pstream::nProcs()
 383         );
 384
 385         gatheredNewParticles[Pstream::myProcNo()] = newParticles_;
 386
 387         // Gather data onto master
 388         Pstream::gatherList(gatheredNewParticles);
 389
 390         // Combine
 391         List<vectorPairScalarPair> combinedNewParticles
 392         (
 393             ListListOps::combine<List<vectorPairScalarPair> >
 394             (
 395                 gatheredNewParticles,
 396                 accessOp<List<vectorPairScalarPair> >()
 397             )
 398         );
 399
 400         if (Pstream::master())
 401         {
 402             newParticles_ = combinedNewParticles;
 403         }
 404
 405         Pstream::scatter(newParticles_);
 406     }
 407
 408     return newParticles_.size();
 409 }
 410
 411
 412 template<class CloudType>
 413 Foam::scalar Foam::InflationInjection<CloudType>::volumeToInject
 414 (
 415     const scalar time0,
 416     const scalar time1
 417 )
 418 {
 419     if ((time0 >= 0.0) && (time0 < duration_))
 420     {
 421         return fraction_*flowRateProfile_().integrate(time0, time1);
 422     }
 423     else
 424     {
 425         return 0.0;
 426     }
 427 }
 428
 429
 430 template<class CloudType>
 431 void Foam::InflationInjection<CloudType>::setPositionAndCell
 432 (
 433     const label parcelI,
 434     const label,
 435     const scalar,
 436     vector& position,
 437     label& cellOwner,
 438     label& tetFaceI,
 439     label& tetPtI
 440 )
 441 {
 442     position = newParticles_[parcelI].first().first();
 443
 444     this->findCellAtPosition
 445     (
 446         cellOwner,
 447         tetFaceI,
 448         tetPtI,
 449         position,
 450         false
 451     );
 452 }
 453
 454
 455 template<class CloudType>
 456 void Foam::InflationInjection<CloudType>::setProperties
 457 (
 458     const label parcelI,
 459     const label,
 460     const scalar,
 461     typename CloudType::parcelType& parcel
 462 )
 463 {
 464     parcel.U() = newParticles_[parcelI].first().second();
 465
 466     parcel.d() = newParticles_[parcelI].second().first();
 467
 468     parcel.dTarget() = newParticles_[parcelI].second().second();
 469 }
 470
 471
 472 template<class CloudType>
 473 bool Foam::InflationInjection<CloudType>::fullyDescribed() const
 474 {
 475     return false;
 476 }
 477
 478
 479 template<class CloudType>
 480 bool Foam::InflationInjection<CloudType>::validInjection(const label)
 481 {
 482     return true;
 483 }
 484
 485
 486 // ************************************************************************* //