src/postProcessing/functionObjects/field/streamLine/streamLineParticle.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 "streamLineParticle.H"
  27 #include "vectorFieldIOField.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33 //    defineParticleTypeNameAndDebug(streamLineParticle, 0);
  34 }
  35
  36
  37 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  38
  39 Foam::scalar Foam::streamLineParticle::calcSubCycleDeltaT
  40 (
  41     trackingData& td,
  42     const scalar dt,
  43     const vector& U
  44 ) const
  45 {
  46     streamLineParticle testParticle(*this);
  47
  48     bool oldKeepParticle = td.keepParticle;
  49     bool oldSwitchProcessor = td.switchProcessor;
  50     scalar fraction = testParticle.trackToFace(position()+dt*U, td);
  51     td.keepParticle = oldKeepParticle;
  52     td.switchProcessor = oldSwitchProcessor;
  53     // Adapt the dt to subdivide the trajectory into substeps.
  54     return dt*fraction/td.nSubCycle_;
  55 }
  56
  57
  58 Foam::vector Foam::streamLineParticle::interpolateFields
  59 (
  60     const trackingData& td,
  61     const point& position,
  62     const label cellI
  63 )
  64 {
  65     if (cellI == -1)
  66     {
  67         FatalErrorIn("streamLineParticle::interpolateFields(..)")
  68             << "Cell:" << cellI << abort(FatalError);
  69     }
  70
  71     sampledScalars_.setSize(td.vsInterp_.size());
  72     forAll(td.vsInterp_, scalarI)
  73     {
  74         sampledScalars_[scalarI].append
  75         (
  76             td.vsInterp_[scalarI].interpolate
  77             (
  78                 position,
  79                 cellI
  80             )
  81         );
  82     }
  83
  84     sampledVectors_.setSize(td.vvInterp_.size());
  85     forAll(td.vvInterp_, vectorI)
  86     {
  87         sampledVectors_[vectorI].append
  88         (
  89             td.vvInterp_[vectorI].interpolate
  90             (
  91                 position,
  92                 cellI
  93             )
  94         );
  95     }
  96
  97     const DynamicList<vector>& U = sampledVectors_[td.UIndex_];
  98
  99     return U.last();
 100 }
 101
 102
 103 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 104
 105 Foam::streamLineParticle::streamLineParticle
 106 (
 107     const polyMesh& mesh,
 108     const vector& position,
 109     const label cellI,
 110     const label lifeTime
 111 )
 112 :
 113     particle(mesh, position, cellI),
 114     lifeTime_(lifeTime)
 115 {}
 116
 117
 118 Foam::streamLineParticle::streamLineParticle
 119 (
 120     const polyMesh& mesh,
 121     Istream& is,
 122     bool readFields
 123 )
 124 :
 125     particle(mesh, is, readFields)
 126 {
 127     if (readFields)
 128     {
 129         //if (is.format() == IOstream::ASCII)
 130         List<scalarList> sampledScalars;
 131         List<vectorList> sampledVectors;
 132
 133         is  >> lifeTime_ >> sampledPositions_ >> sampledScalars
 134             >> sampledVectors;
 135
 136         sampledScalars_.setSize(sampledScalars.size());
 137         forAll(sampledScalars, i)
 138         {
 139             sampledScalars_[i].transfer(sampledScalars[i]);
 140         }
 141         sampledVectors_.setSize(sampledVectors.size());
 142         forAll(sampledVectors, i)
 143         {
 144             sampledVectors_[i].transfer(sampledVectors[i]);
 145         }
 146     }
 147
 148     // Check state of Istream
 149     is.check
 150     (
 151         "streamLineParticle::streamLineParticle"
 152         "(const Cloud<streamLineParticle>&, Istream&, bool)"
 153     );
 154 }
 155
 156
 157 Foam::streamLineParticle::streamLineParticle
 158 (
 159     const streamLineParticle& p
 160 )
 161 :
 162     particle(p),
 163     lifeTime_(p.lifeTime_),
 164     sampledPositions_(p.sampledPositions_),
 165     sampledScalars_(p.sampledScalars_)
 166 {}
 167
 168
 169 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 170
 171 bool Foam::streamLineParticle::move
 172 (
 173     trackingData& td,
 174     const scalar trackTime
 175 )
 176 {
 177     streamLineParticle& p = static_cast<streamLineParticle&>(*this);
 178
 179     td.switchProcessor = false;
 180     td.keepParticle = true;
 181
 182     scalar tEnd = (1.0 - stepFraction())*trackTime;
 183     scalar maxDt = mesh_.bounds().mag();
 184
 185     while
 186     (
 187         td.keepParticle
 188     && !td.switchProcessor
 189     && lifeTime_ > 0
 190     )
 191     {
 192         // set the lagrangian time-step
 193         scalar dt = maxDt;
 194
 195         // Cross cell in steps:
 196         // - at subiter 0 calculate dt to cross cell in nSubCycle steps
 197         // - at the last subiter do all of the remaining track
 198         for (label subIter = 0; subIter < td.nSubCycle_; subIter++)
 199         {
 200             --lifeTime_;
 201
 202             // Store current position and sampled velocity.
 203             sampledPositions_.append(position());
 204             vector U = interpolateFields(td, position(), cell());
 205
 206             if (!td.trackForward_)
 207             {
 208                 U = -U;
 209             }
 210
 211             scalar magU = mag(U);
 212
 213             if (magU < SMALL)
 214             {
 215                 // Stagnant particle. Might as well stop
 216                 lifeTime_ = 0;
 217                 break;
 218             }
 219
 220             U /= magU;
 221
 222             if (subIter == 0 && td.nSubCycle_ > 1)
 223             {
 224                 // Adapt dt to cross cell in a few steps
 225                 dt = calcSubCycleDeltaT(td, dt, U);
 226             }
 227             else if (subIter == td.nSubCycle_ - 1)
 228             {
 229                 // Do full step on last subcycle
 230                 dt = maxDt;
 231             }
 232
 233
 234             scalar fraction = trackToFace(position() + dt*U, td);
 235             dt *= fraction;
 236
 237             tEnd -= dt;
 238             stepFraction() = 1.0 - tEnd/trackTime;
 239
 240             if (tEnd <= ROOTVSMALL)
 241             {
 242                 // Force removal
 243                 lifeTime_ = 0;
 244             }
 245
 246             if
 247             (
 248                 face() != -1
 249             || !td.keepParticle
 250             ||  td.switchProcessor
 251             ||  lifeTime_ == 0
 252             )
 253             {
 254                 break;
 255             }
 256         }
 257     }
 258
 259
 260     if (!td.keepParticle || lifeTime_ == 0)
 261     {
 262         if (lifeTime_ == 0)
 263         {
 264             if (debug)
 265             {
 266                 Pout<< "streamLineParticle : Removing stagnant particle:"
 267                     << p << " sampled positions:" << sampledPositions_.size()
 268                     << endl;
 269             }
 270             td.keepParticle = false;
 271         }
 272         else
 273         {
 274             // Normal exit. Store last position and fields
 275             sampledPositions_.append(position());
 276             interpolateFields(td, position(), cell());
 277
 278             if (debug)
 279             {
 280                 Pout<< "streamLineParticle : Removing particle:"
 281                     << p << " sampled positions:" << sampledPositions_.size()
 282                     << endl;
 283             }
 284         }
 285
 286         // Transfer particle data into trackingData.
 287         //td.allPositions_.append(sampledPositions_);
 288         td.allPositions_.append(vectorList());
 289         vectorList& top = td.allPositions_.last();
 290         top.transfer(sampledPositions_);
 291
 292         forAll(sampledScalars_, i)
 293         {
 294             //td.allScalars_[i].append(sampledScalars_[i]);
 295             td.allScalars_[i].append(scalarList());
 296             scalarList& top = td.allScalars_[i].last();
 297             top.transfer(sampledScalars_[i]);
 298         }
 299         forAll(sampledVectors_, i)
 300         {
 301             //td.allVectors_[i].append(sampledVectors_[i]);
 302             td.allVectors_[i].append(vectorList());
 303             vectorList& top = td.allVectors_[i].last();
 304             top.transfer(sampledVectors_[i]);
 305         }
 306     }
 307
 308     return td.keepParticle;
 309 }
 310
 311
 312 bool Foam::streamLineParticle::hitPatch
 313 (
 314     const polyPatch&,
 315     trackingData& td,
 316     const label patchI,
 317     const scalar trackFraction,
 318     const tetIndices& tetIs
 319 )
 320 {
 321     // Disable generic patch interaction
 322     return false;
 323 }
 324
 325
 326 void Foam::streamLineParticle::hitWedgePatch
 327 (
 328     const wedgePolyPatch& pp,
 329     trackingData& td
 330 )
 331 {
 332     // Remove particle
 333     td.keepParticle = false;
 334 }
 335
 336
 337 void Foam::streamLineParticle::hitSymmetryPatch
 338 (
 339     const symmetryPolyPatch& pp,
 340     trackingData& td
 341 )
 342 {
 343     // Remove particle
 344     td.keepParticle = false;
 345 }
 346
 347
 348 void Foam::streamLineParticle::hitCyclicPatch
 349 (
 350     const cyclicPolyPatch& pp,
 351     trackingData& td
 352 )
 353 {
 354     // Remove particle
 355     td.keepParticle = false;
 356 }
 357
 358
 359 void Foam::streamLineParticle::hitProcessorPatch
 360 (
 361     const processorPolyPatch&,
 362     trackingData& td
 363 )
 364 {
 365     // Switch particle
 366     td.switchProcessor = true;
 367 }
 368
 369
 370 void Foam::streamLineParticle::hitWallPatch
 371 (
 372     const wallPolyPatch& wpp,
 373     trackingData& td,
 374     const tetIndices&
 375 )
 376 {
 377     // Remove particle
 378     td.keepParticle = false;
 379 }
 380
 381
 382 void Foam::streamLineParticle::hitPatch
 383 (
 384     const polyPatch& wpp,
 385     trackingData& td
 386 )
 387 {
 388     // Remove particle
 389     td.keepParticle = false;
 390 }
 391
 392
 393 void Foam::streamLineParticle::readFields(Cloud<streamLineParticle>& c)
 394 {
 395     if (!c.size())
 396     {
 397         return;
 398     }
 399
 400     particle::readFields(c);
 401
 402     IOField<label> lifeTime
 403     (
 404         c.fieldIOobject("lifeTime", IOobject::MUST_READ)
 405     );
 406     c.checkFieldIOobject(c, lifeTime);
 407
 408     vectorFieldIOField sampledPositions
 409     (
 410         c.fieldIOobject("sampledPositions", IOobject::MUST_READ)
 411     );
 412     c.checkFieldIOobject(c, sampledPositions);
 413
 414 //    vectorFieldIOField sampleVelocity
 415 //    (
 416 //        c.fieldIOobject("sampleVelocity", IOobject::MUST_READ)
 417 //    );
 418 //    c.checkFieldIOobject(c, sampleVelocity);
 419
 420     label i = 0;
 421     forAllIter(Cloud<streamLineParticle>, c, iter)
 422     {
 423         iter().lifeTime_ = lifeTime[i];
 424         iter().sampledPositions_.transfer(sampledPositions[i]);
 425 //        iter().sampleVelocity_.transfer(sampleVelocity[i]);
 426         i++;
 427     }
 428 }
 429
 430
 431 void Foam::streamLineParticle::writeFields(const Cloud<streamLineParticle>& c)
 432 {
 433     particle::writeFields(c);
 434
 435     label np =  c.size();
 436
 437     IOField<label> lifeTime
 438     (
 439         c.fieldIOobject("lifeTime", IOobject::NO_READ),
 440         np
 441     );
 442     vectorFieldIOField sampledPositions
 443     (
 444         c.fieldIOobject("sampledPositions", IOobject::NO_READ),
 445         np
 446     );
 447 //    vectorFieldIOField sampleVelocity
 448 //    (
 449 //        c.fieldIOobject("sampleVelocity", IOobject::NO_READ),
 450 //        np
 451 //    );
 452
 453     label i = 0;
 454     forAllConstIter(Cloud<streamLineParticle>, c, iter)
 455     {
 456         lifeTime[i] = iter().lifeTime_;
 457         sampledPositions[i] = iter().sampledPositions_;
 458 //        sampleVelocity[i] = iter().sampleVelocity_;
 459         i++;
 460     }
 461
 462     lifeTime.write();
 463     sampledPositions.write();
 464 //    sampleVelocity.write();
 465 }
 466
 467
 468 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 469
 470 Foam::Ostream& Foam::operator<<(Ostream& os, const streamLineParticle& p)
 471 {
 472     os  << static_cast<const particle&>(p)
 473         << token::SPACE << p.lifeTime_
 474         << token::SPACE << p.sampledPositions_
 475         << token::SPACE << p.sampledScalars_
 476         << token::SPACE << p.sampledVectors_;
 477
 478     // Check state of Ostream
 479     os.check("Ostream& operator<<(Ostream&, const streamLineParticle&)");
 480
 481     return os;
 482 }
 483
 484
 485 // ************************************************************************* //