src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.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 "KinematicCloud.H"
  27 #include "IntegrationScheme.H"
  28 #include "interpolation.H"
  29 #include "subCycleTime.H"
  30
  31 #include "DispersionModel.H"
  32 #include "InjectionModel.H"
  33 #include "PatchInteractionModel.H"
  34 #include "SurfaceFilmModel.H"
  35
  36 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
  37
  38 template<class CloudType>
  39 void Foam::KinematicCloud<CloudType>::setModels()
  40 {
  41     dispersionModel_.reset
  42     (
  43         DispersionModel<KinematicCloud<CloudType> >::New
  44         (
  45             subModelProperties_,
  46             *this
  47         ).ptr()
  48     );
  49
  50     injectionModel_.reset
  51     (
  52         InjectionModel<KinematicCloud<CloudType> >::New
  53         (
  54             subModelProperties_,
  55             *this
  56         ).ptr()
  57     );
  58
  59     patchInteractionModel_.reset
  60     (
  61         PatchInteractionModel<KinematicCloud<CloudType> >::New
  62         (
  63             subModelProperties_,
  64             *this
  65         ).ptr()
  66     );
  67
  68     surfaceFilmModel_.reset
  69     (
  70         SurfaceFilmModel<KinematicCloud<CloudType> >::New
  71         (
  72             subModelProperties_,
  73             *this,
  74             g_
  75         ).ptr()
  76     );
  77
  78     UIntegrator_.reset
  79     (
  80         vectorIntegrationScheme::New
  81         (
  82             "U",
  83             solution_.integrationSchemes()
  84         ).ptr()
  85     );
  86 }
  87
  88
  89 template<class CloudType>
  90 template<class TrackData>
  91 void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
  92 {
  93     if (solution_.steadyState())
  94     {
  95         td.cloud().storeState();
  96
  97         td.cloud().preEvolve();
  98
  99         evolveCloud(td);
 100
 101         if (solution_.coupled())
 102         {
 103             td.cloud().relaxSources(td.cloud().cloudCopy());
 104         }
 105     }
 106     else
 107     {
 108         td.cloud().preEvolve();
 109
 110         evolveCloud(td);
 111
 112         if (solution_.coupled())
 113         {
 114             td.cloud().scaleSources();
 115         }
 116     }
 117
 118     td.cloud().info();
 119
 120     td.cloud().postEvolve();
 121
 122     if (solution_.steadyState())
 123     {
 124         td.cloud().restoreState();
 125     }
 126 }
 127
 128
 129 template<class CloudType>
 130 void Foam::KinematicCloud<CloudType>::buildCellOccupancy()
 131 {
 132     if (cellOccupancyPtr_.empty())
 133     {
 134         cellOccupancyPtr_.reset
 135         (
 136             new List<DynamicList<parcelType*> >(mesh_.nCells())
 137         );
 138     }
 139     else if (cellOccupancyPtr_().size() != mesh_.nCells())
 140     {
 141         // If the size of the mesh has changed, reset the
 142         // cellOccupancy size
 143
 144         cellOccupancyPtr_().setSize(mesh_.nCells());
 145     }
 146
 147     List<DynamicList<parcelType*> >& cellOccupancy = cellOccupancyPtr_();
 148
 149     forAll(cellOccupancy, cO)
 150     {
 151         cellOccupancy[cO].clear();
 152     }
 153
 154     forAllIter(typename KinematicCloud<CloudType>, *this, iter)
 155     {
 156         cellOccupancy[iter().cell()].append(&iter());
 157     }
 158 }
 159
 160
 161 template<class CloudType>
 162 void Foam::KinematicCloud<CloudType>::updateCellOccupancy()
 163 {
 164     // Only build the cellOccupancy if the pointer is set, i.e. it has
 165     // been requested before.
 166
 167     if (cellOccupancyPtr_.valid())
 168     {
 169         buildCellOccupancy();
 170     }
 171 }
 172
 173
 174 template<class CloudType>
 175 template<class TrackData>
 176 void Foam::KinematicCloud<CloudType>::evolveCloud(TrackData& td)
 177 {
 178     if (solution_.coupled())
 179     {
 180         td.cloud().resetSourceTerms();
 181     }
 182
 183     if (solution_.transient())
 184     {
 185         label preInjectionSize = this->size();
 186
 187         this->surfaceFilm().inject(td);
 188
 189         // Update the cellOccupancy if the size of the cloud has changed
 190         // during the injection.
 191         if (preInjectionSize != this->size())
 192         {
 193             updateCellOccupancy();
 194
 195             preInjectionSize = this->size();
 196         }
 197         this->injection().inject(td);
 198
 199
 200         // Assume that motion will update the cellOccupancy as necessary
 201         // before it is required.
 202         td.cloud().motion(td);
 203     }
 204     else
 205     {
 206 //        this->surfaceFilm().injectSteadyState(td);
 207
 208         this->injection().injectSteadyState(td, solution_.trackTime());
 209
 210         td.part() = TrackData::tpLinearTrack;
 211         CloudType::move(td,  solution_.trackTime());
 212     }
 213 }
 214
 215
 216 template<class CloudType>
 217 void Foam::KinematicCloud<CloudType>::postEvolve()
 218 {
 219     Info<< endl;
 220
 221     if (debug)
 222     {
 223         this->writePositions();
 224     }
 225
 226     this->dispersion().cacheFields(false);
 227
 228     forces_.cacheFields(false);
 229     functions_.postEvolve();
 230
 231     solution_.nextIter();
 232 }
 233
 234
 235 template<class CloudType>
 236 void Foam::KinematicCloud<CloudType>::cloudReset(KinematicCloud<CloudType>& c)
 237 {
 238     CloudType::cloudReset(c);
 239
 240     rndGen_ = c.rndGen_;
 241
 242     forces_.transfer(c.forces_);
 243
 244     functions_.transfer(c.functions_);
 245
 246     dispersionModel_.reset(c.dispersionModel_.ptr());
 247     injectionModel_.reset(c.injectionModel_.ptr());
 248     patchInteractionModel_.reset(c.patchInteractionModel_.ptr());
 249     surfaceFilmModel_.reset(c.surfaceFilmModel_.ptr());
 250
 251     UIntegrator_.reset(c.UIntegrator_.ptr());
 252 }
 253
 254
 255 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 256
 257 template<class CloudType>
 258 Foam::KinematicCloud<CloudType>::KinematicCloud
 259 (
 260     const word& cloudName,
 261     const volScalarField& rho,
 262     const volVectorField& U,
 263     const volScalarField& mu,
 264     const dimensionedVector& g,
 265     bool readFields
 266 )
 267 :
 268     CloudType(rho.mesh(), cloudName, false),
 269     kinematicCloud(),
 270     cloudCopyPtr_(NULL),
 271     mesh_(rho.mesh()),
 272     particleProperties_
 273     (
 274         IOobject
 275         (
 276             cloudName + "Properties",
 277             rho.mesh().time().constant(),
 278             rho.mesh(),
 279             IOobject::MUST_READ_IF_MODIFIED,
 280             IOobject::NO_WRITE
 281         )
 282     ),
 283     solution_(mesh_, particleProperties_.subDict("solution")),
 284     constProps_(particleProperties_, solution_.active()),
 285     subModelProperties_
 286     (
 287         particleProperties_.subOrEmptyDict("subModels", solution_.active())
 288     ),
 289     rndGen_
 290     (
 291         label(0),
 292         solution_.steadyState() ?
 293         particleProperties_.lookupOrDefault<label>("randomSampleSize", 100000)
 294       : -1
 295     ),
 296     cellOccupancyPtr_(),
 297     rho_(rho),
 298     U_(U),
 299     mu_(mu),
 300     g_(g),
 301     pAmbient_(0.0),
 302     forces_
 303     (
 304         *this,
 305         mesh_,
 306         subModelProperties_.subOrEmptyDict
 307         (
 308             "particleForces",
 309             solution_.active()
 310         ),
 311         solution_.active()
 312     ),
 313     functions_
 314     (
 315         *this,
 316         particleProperties_.subOrEmptyDict("cloudFunctions"),
 317         solution_.active()
 318     ),
 319     dispersionModel_(NULL),
 320     injectionModel_(NULL),
 321     patchInteractionModel_(NULL),
 322     surfaceFilmModel_(NULL),
 323     UIntegrator_(NULL),
 324     UTrans_
 325     (
 326         new DimensionedField<vector, volMesh>
 327         (
 328             IOobject
 329             (
 330                 this->name() + "UTrans",
 331                 this->db().time().timeName(),
 332                 this->db(),
 333                 IOobject::READ_IF_PRESENT,
 334                 IOobject::AUTO_WRITE
 335             ),
 336             mesh_,
 337             dimensionedVector("zero", dimMass*dimVelocity, vector::zero)
 338         )
 339     ),
 340     UCoeff_
 341     (
 342         new DimensionedField<scalar, volMesh>
 343         (
 344             IOobject
 345             (
 346                 this->name() + "UCoeff",
 347                 this->db().time().timeName(),
 348                 this->db(),
 349                 IOobject::READ_IF_PRESENT,
 350                 IOobject::AUTO_WRITE
 351             ),
 352             mesh_,
 353             dimensionedScalar("zero",  dimMass, 0.0)
 354         )
 355     )
 356 {
 357     if (solution_.active())
 358     {
 359         setModels();
 360
 361         if (readFields)
 362         {
 363             parcelType::readFields(*this);
 364         }
 365     }
 366
 367     if (solution_.resetSourcesOnStartup())
 368     {
 369         resetSourceTerms();
 370     }
 371 }
 372
 373
 374 template<class CloudType>
 375 Foam::KinematicCloud<CloudType>::KinematicCloud
 376 (
 377     KinematicCloud<CloudType>& c,
 378     const word& name
 379 )
 380 :
 381     CloudType(c.mesh_, name, c),
 382     kinematicCloud(),
 383     cloudCopyPtr_(NULL),
 384     mesh_(c.mesh_),
 385     particleProperties_(c.particleProperties_),
 386     solution_(c.solution_),
 387     constProps_(c.constProps_),
 388     subModelProperties_(c.subModelProperties_),
 389     rndGen_(c.rndGen_, true),
 390     cellOccupancyPtr_(NULL),
 391     rho_(c.rho_),
 392     U_(c.U_),
 393     mu_(c.mu_),
 394     g_(c.g_),
 395     pAmbient_(c.pAmbient_),
 396     forces_(c.forces_),
 397     functions_(c.functions_),
 398     dispersionModel_(c.dispersionModel_->clone()),
 399     injectionModel_(c.injectionModel_->clone()),
 400     patchInteractionModel_(c.patchInteractionModel_->clone()),
 401     surfaceFilmModel_(c.surfaceFilmModel_->clone()),
 402     UIntegrator_(c.UIntegrator_->clone()),
 403     UTrans_
 404     (
 405         new DimensionedField<vector, volMesh>
 406         (
 407             IOobject
 408             (
 409                 this->name() + "UTrans",
 410                 this->db().time().timeName(),
 411                 this->db(),
 412                 IOobject::NO_READ,
 413                 IOobject::NO_WRITE,
 414                 false
 415             ),
 416             c.UTrans_()
 417         )
 418     ),
 419     UCoeff_
 420     (
 421         new DimensionedField<scalar, volMesh>
 422         (
 423             IOobject
 424             (
 425                 name + "UCoeff",
 426                 this->db().time().timeName(),
 427                 this->db(),
 428                 IOobject::NO_READ,
 429                 IOobject::NO_WRITE,
 430                 false
 431             ),
 432             c.UCoeff_()
 433         )
 434     )
 435 {}
 436
 437
 438 template<class CloudType>
 439 Foam::KinematicCloud<CloudType>::KinematicCloud
 440 (
 441     const fvMesh& mesh,
 442     const word& name,
 443     const KinematicCloud<CloudType>& c
 444 )
 445 :
 446     CloudType(mesh, name, IDLList<parcelType>()),
 447     kinematicCloud(),
 448     cloudCopyPtr_(NULL),
 449     mesh_(mesh),
 450     particleProperties_
 451     (
 452         IOobject
 453         (
 454             name + "Properties",
 455             mesh.time().constant(),
 456             mesh,
 457             IOobject::NO_READ,
 458             IOobject::NO_WRITE,
 459             false
 460         )
 461     ),
 462     solution_(mesh),
 463     constProps_(),
 464     subModelProperties_(dictionary::null),
 465     rndGen_(0, 0),
 466     cellOccupancyPtr_(NULL),
 467     rho_(c.rho_),
 468     U_(c.U_),
 469     mu_(c.mu_),
 470     g_(c.g_),
 471     pAmbient_(c.pAmbient_),
 472     forces_(*this, mesh),
 473     functions_(*this),
 474     dispersionModel_(NULL),
 475     injectionModel_(NULL),
 476     patchInteractionModel_(NULL),
 477     surfaceFilmModel_(NULL),
 478     UIntegrator_(NULL),
 479     UTrans_(NULL),
 480     UCoeff_(NULL)
 481 {}
 482
 483
 484 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 485
 486 template<class CloudType>
 487 Foam::KinematicCloud<CloudType>::~KinematicCloud()
 488 {}
 489
 490
 491 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 492
 493 template<class CloudType>
 494 bool Foam::KinematicCloud<CloudType>::hasWallImpactDistance() const
 495 {
 496     return true;
 497 }
 498
 499
 500 template<class CloudType>
 501 void Foam::KinematicCloud<CloudType>::setParcelThermoProperties
 502 (
 503     parcelType& parcel,
 504     const scalar lagrangianDt
 505 )
 506 {
 507     parcel.rho() = constProps_.rho0();
 508 }
 509
 510
 511 template<class CloudType>
 512 void Foam::KinematicCloud<CloudType>::checkParcelProperties
 513 (
 514     parcelType& parcel,
 515     const scalar lagrangianDt,
 516     const bool fullyDescribed
 517 )
 518 {
 519     const scalar carrierDt = mesh_.time().deltaTValue();
 520     parcel.stepFraction() = (carrierDt - lagrangianDt)/carrierDt;
 521     parcel.typeId() = constProps_.parcelTypeId();
 522 }
 523
 524
 525 template<class CloudType>
 526 void Foam::KinematicCloud<CloudType>::storeState()
 527 {
 528     cloudCopyPtr_.reset
 529     (
 530         static_cast<KinematicCloud<CloudType>*>
 531         (
 532             clone(this->name() + "Copy").ptr()
 533         )
 534     );
 535 }
 536
 537
 538 template<class CloudType>
 539 void Foam::KinematicCloud<CloudType>::restoreState()
 540 {
 541     cloudReset(cloudCopyPtr_());
 542     cloudCopyPtr_.clear();
 543 }
 544
 545
 546 template<class CloudType>
 547 void Foam::KinematicCloud<CloudType>::resetSourceTerms()
 548 {
 549     UTrans().field() = vector::zero;
 550     UCoeff().field() = 0.0;
 551 }
 552
 553
 554 template<class CloudType>
 555 template<class Type>
 556 void Foam::KinematicCloud<CloudType>::relax
 557 (
 558     DimensionedField<Type, volMesh>& field,
 559     const DimensionedField<Type, volMesh>& field0,
 560     const word& name
 561 ) const
 562 {
 563     const scalar coeff = solution_.relaxCoeff(name);
 564     field = field0 + coeff*(field - field0);
 565 }
 566
 567
 568 template<class CloudType>
 569 template<class Type>
 570 void Foam::KinematicCloud<CloudType>::scale
 571 (
 572     DimensionedField<Type, volMesh>& field,
 573     const word& name
 574 ) const
 575 {
 576     const scalar coeff = solution_.relaxCoeff(name);
 577     field *= coeff;
 578 }
 579
 580
 581 template<class CloudType>
 582 void Foam::KinematicCloud<CloudType>::relaxSources
 583 (
 584     const KinematicCloud<CloudType>& cloudOldTime
 585 )
 586 {
 587     this->relax(UTrans_(), cloudOldTime.UTrans(), "U");
 588     this->relax(UCoeff_(), cloudOldTime.UCoeff(), "U");
 589 }
 590
 591
 592 template<class CloudType>
 593 void Foam::KinematicCloud<CloudType>::scaleSources()
 594 {
 595     this->scale(UTrans_(), "U");
 596     this->scale(UCoeff_(), "U");
 597 }
 598
 599
 600 template<class CloudType>
 601 void Foam::KinematicCloud<CloudType>::preEvolve()
 602 {
 603     Info<< "\nSolving cloud " << this->name() << endl;
 604
 605     this->dispersion().cacheFields(true);
 606     forces_.cacheFields(true);
 607     updateCellOccupancy();
 608
 609     pAmbient_ = constProps_.dict().template
 610         lookupOrDefault<scalar>("pAmbient", pAmbient_);
 611
 612     functions_.preEvolve();
 613 }
 614
 615
 616 template<class CloudType>
 617 void Foam::KinematicCloud<CloudType>::evolve()
 618 {
 619     if (solution_.canEvolve())
 620     {
 621         typename parcelType::template
 622             TrackingData<KinematicCloud<CloudType> > td(*this);
 623
 624         solve(td);
 625     }
 626 }
 627
 628
 629 template<class CloudType>
 630 template<class TrackData>
 631 void  Foam::KinematicCloud<CloudType>::motion(TrackData& td)
 632 {
 633     td.part() = TrackData::tpLinearTrack;
 634     CloudType::move(td,  solution_.trackTime());
 635
 636     updateCellOccupancy();
 637 }
 638
 639
 640 template<class CloudType>
 641 void Foam::KinematicCloud<CloudType>::info() const
 642 {
 643     vector linearMomentum = linearMomentumOfSystem();
 644     reduce(linearMomentum, sumOp<vector>());
 645
 646     scalar linearKineticEnergy = linearKineticEnergyOfSystem();
 647     reduce(linearKineticEnergy, sumOp<scalar>());
 648
 649     scalar rotationalKineticEnergy = rotationalKineticEnergyOfSystem();
 650     reduce(rotationalKineticEnergy, sumOp<scalar>());
 651
 652     Info<< "Cloud: " << this->name() << nl
 653         << "    Current number of parcels       = "
 654         << returnReduce(this->size(), sumOp<label>()) << nl
 655         << "    Current mass in system          = "
 656         << returnReduce(massInSystem(), sumOp<scalar>()) << nl
 657         << "    Linear momentum                 = "
 658         << linearMomentum << nl
 659         << "   |Linear momentum|                = "
 660         << mag(linearMomentum) << nl
 661         << "    Linear kinetic energy           = "
 662         << linearKineticEnergy << nl
 663         << "    Rotational kinetic energy       = "
 664         << rotationalKineticEnergy << nl;
 665
 666     this->injection().info(Info);
 667     this->surfaceFilm().info(Info);
 668     this->patchInteraction().info(Info);
 669 }
 670
 671
 672 // ************************************************************************* //