streamLineParticle.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2013 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  particle 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  const label faceI
64 )
65 {
66  if (cellI == -1)
67  {
68  FatalErrorIn("streamLineParticle::interpolateFields(..)")
69  << "Cell:" << cellI << abort(FatalError);
70  }
71 
72  sampledScalars_.setSize(td.vsInterp_.size());
73  forAll(td.vsInterp_, scalarI)
74  {
75  sampledScalars_[scalarI].append
76  (
77  td.vsInterp_[scalarI].interpolate
78  (
79  position,
80  cellI,
81  faceI
82  )
83  );
84  }
85 
86  sampledVectors_.setSize(td.vvInterp_.size());
87  forAll(td.vvInterp_, vectorI)
88  {
89  sampledVectors_[vectorI].append
90  (
91  td.vvInterp_[vectorI].interpolate
92  (
93  position,
94  cellI,
95  faceI
96  )
97  );
98  }
99 
100  const DynamicList<vector>& U = sampledVectors_[td.UIndex_];
101 
102  return U.last();
103 }
104 
105 
106 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
107 
109 (
110  const polyMesh& mesh,
111  const vector& position,
112  const label cellI,
113  const label lifeTime
114 )
115 :
116  particle(mesh, position, cellI),
117  lifeTime_(lifeTime)
118 {}
119 
120 
122 (
123  const polyMesh& mesh,
124  Istream& is,
125  bool readFields
126 )
127 :
128  particle(mesh, is, readFields)
129 {
130  if (readFields)
131  {
132  //if (is.format() == IOstream::ASCII)
133  List<scalarList> sampledScalars;
134  List<vectorList> sampledVectors;
135 
136  is >> lifeTime_ >> sampledPositions_ >> sampledScalars
137  >> sampledVectors;
138 
139  sampledScalars_.setSize(sampledScalars.size());
140  forAll(sampledScalars, i)
141  {
142  sampledScalars_[i].transfer(sampledScalars[i]);
143  }
144  sampledVectors_.setSize(sampledVectors.size());
145  forAll(sampledVectors, i)
146  {
147  sampledVectors_[i].transfer(sampledVectors[i]);
148  }
149  }
150 
151  // Check state of Istream
152  is.check
153  (
154  "streamLineParticle::streamLineParticle"
155  "(const Cloud<streamLineParticle>&, Istream&, bool)"
156  );
157 }
158 
159 
161 (
162  const streamLineParticle& p
163 )
164 :
165  particle(p),
166  lifeTime_(p.lifeTime_),
167  sampledPositions_(p.sampledPositions_),
168  sampledScalars_(p.sampledScalars_)
169 {}
170 
171 
172 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
173 
175 (
176  trackingData& td,
177  const scalar trackTime
178 )
179 {
180  streamLineParticle& p = static_cast<streamLineParticle&>(*this);
181 
182  td.switchProcessor = false;
183  td.keepParticle = true;
184 
185  scalar tEnd = (1.0 - stepFraction())*trackTime;
186  scalar maxDt = mesh_.bounds().mag();
187 
188  while
189  (
190  td.keepParticle
191  && !td.switchProcessor
192  && lifeTime_ > 0
193  )
194  {
195  // set the lagrangian time-step
196  scalar dt = maxDt;
197 
198  // Cross cell in steps:
199  // - at subiter 0 calculate dt to cross cell in nSubCycle steps
200  // - at the last subiter do all of the remaining track
201  for (label subIter = 0; subIter < td.nSubCycle_; subIter++)
202  {
203  --lifeTime_;
204 
205  // Store current position and sampled velocity.
206  sampledPositions_.append(position());
207  vector U = interpolateFields(td, position(), cell(), face());
208 
209  if (!td.trackForward_)
210  {
211  U = -U;
212  }
213 
214  scalar magU = mag(U);
215 
216  if (magU < SMALL)
217  {
218  // Stagnant particle. Might as well stop
219  lifeTime_ = 0;
220  break;
221  }
222 
223  U /= magU;
224 
225  if (td.trackLength_ < GREAT)
226  {
227  dt = td.trackLength_;
228  //Pout<< " subiteration " << subIter
229  // << " : fixed length: updated dt:" << dt << endl;
230  }
231  else if (subIter == 0 && td.nSubCycle_ > 1)
232  {
233  // Adapt dt to cross cell in a few steps
234  dt = calcSubCycleDeltaT(td, dt, U);
235  }
236  else if (subIter == td.nSubCycle_ - 1)
237  {
238  // Do full step on last subcycle
239  dt = maxDt;
240  }
241 
242 
243  scalar fraction = trackToFace(position() + dt*U, td);
244  dt *= fraction;
245 
246  tEnd -= dt;
247  stepFraction() = 1.0 - tEnd/trackTime;
248 
249  if (tEnd <= ROOTVSMALL)
250  {
251  // Force removal
252  lifeTime_ = 0;
253  }
254 
255  if
256  (
257  face() != -1
258  || !td.keepParticle
259  || td.switchProcessor
260  || lifeTime_ == 0
261  )
262  {
263  break;
264  }
265  }
266  }
267 
268 
269  if (!td.keepParticle || lifeTime_ == 0)
270  {
271  if (lifeTime_ == 0)
272  {
273  if (debug)
274  {
275  Pout<< "streamLineParticle : Removing stagnant particle:"
276  << p.position()
277  << " sampled positions:" << sampledPositions_.size()
278  << endl;
279  }
280  td.keepParticle = false;
281  }
282  else
283  {
284  // Normal exit. Store last position and fields
285  sampledPositions_.append(position());
286  interpolateFields(td, position(), cell(), face());
287 
288  if (debug)
289  {
290  Pout<< "streamLineParticle : Removing particle:"
291  << p.position()
292  << " sampled positions:" << sampledPositions_.size()
293  << endl;
294  }
295  }
296 
297  // Transfer particle data into trackingData.
298  //td.allPositions_.append(sampledPositions_);
299  td.allPositions_.append(vectorList());
300  vectorList& top = td.allPositions_.last();
301  top.transfer(sampledPositions_);
302 
303  forAll(sampledScalars_, i)
304  {
305  //td.allScalars_[i].append(sampledScalars_[i]);
306  td.allScalars_[i].append(scalarList());
307  scalarList& top = td.allScalars_[i].last();
308  top.transfer(sampledScalars_[i]);
309  }
310  forAll(sampledVectors_, i)
311  {
312  //td.allVectors_[i].append(sampledVectors_[i]);
313  td.allVectors_[i].append(vectorList());
314  vectorList& top = td.allVectors_[i].last();
315  top.transfer(sampledVectors_[i]);
316  }
317  }
318 
319  return td.keepParticle;
320 }
321 
322 
324 (
325  const polyPatch&,
326  trackingData& td,
327  const label patchI,
328  const scalar trackFraction,
329  const tetIndices& tetIs
330 )
331 {
332  // Disable generic patch interaction
333  return false;
334 }
335 
336 
338 (
339  const wedgePolyPatch& pp,
340  trackingData& td
341 )
342 {
343  // Remove particle
344  td.keepParticle = false;
345 }
346 
347 
349 (
350  const symmetryPlanePolyPatch& pp,
351  trackingData& td
352 )
353 {
354  // Remove particle
355  td.keepParticle = false;
356 }
357 
358 
360 (
361  const symmetryPolyPatch& pp,
362  trackingData& td
363 )
364 {
365  // Remove particle
366  td.keepParticle = false;
367 }
368 
369 
371 (
372  const cyclicPolyPatch& pp,
373  trackingData& td
374 )
375 {
376  // Remove particle
377  td.keepParticle = false;
378 }
379 
380 
382 (
383  const processorPolyPatch&,
384  trackingData& td
385 )
386 {
387  // Switch particle
388  td.switchProcessor = true;
389 }
390 
391 
393 (
394  const wallPolyPatch& wpp,
395  trackingData& td,
396  const tetIndices&
397 )
398 {
399  // Remove particle
400  td.keepParticle = false;
401 }
402 
403 
405 (
406  const polyPatch& wpp,
407  trackingData& td
408 )
409 {
410  // Remove particle
411  td.keepParticle = false;
412 }
413 
414 
416 {
417  if (!c.size())
418  {
419  return;
420  }
421 
423 
424  IOField<label> lifeTime
425  (
426  c.fieldIOobject("lifeTime", IOobject::MUST_READ)
427  );
428  c.checkFieldIOobject(c, lifeTime);
429 
430  vectorFieldIOField sampledPositions
431  (
432  c.fieldIOobject("sampledPositions", IOobject::MUST_READ)
433  );
434  c.checkFieldIOobject(c, sampledPositions);
435 
436 // vectorFieldIOField sampleVelocity
437 // (
438 // c.fieldIOobject("sampleVelocity", IOobject::MUST_READ)
439 // );
440 // c.checkFieldIOobject(c, sampleVelocity);
441 
442  label i = 0;
444  {
445  iter().lifeTime_ = lifeTime[i];
446  iter().sampledPositions_.transfer(sampledPositions[i]);
447 // iter().sampleVelocity_.transfer(sampleVelocity[i]);
448  i++;
449  }
450 }
451 
452 
454 {
456 
457  label np = c.size();
458 
459  IOField<label> lifeTime
460  (
461  c.fieldIOobject("lifeTime", IOobject::NO_READ),
462  np
463  );
464  vectorFieldIOField sampledPositions
465  (
466  c.fieldIOobject("sampledPositions", IOobject::NO_READ),
467  np
468  );
469 // vectorFieldIOField sampleVelocity
470 // (
471 // c.fieldIOobject("sampleVelocity", IOobject::NO_READ),
472 // np
473 // );
474 
475  label i = 0;
477  {
478  lifeTime[i] = iter().lifeTime_;
479  sampledPositions[i] = iter().sampledPositions_;
480 // sampleVelocity[i] = iter().sampleVelocity_;
481  i++;
482  }
483 
484  lifeTime.write();
485  sampledPositions.write();
486 // sampleVelocity.write();
487 }
488 
489 
490 // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
491 
493 {
494  os << static_cast<const particle&>(p)
495  << token::SPACE << p.lifeTime_
496  << token::SPACE << p.sampledPositions_
497  << token::SPACE << p.sampledScalars_
498  << token::SPACE << p.sampledVectors_;
499 
500  // Check state of Ostream
501  os.check("Ostream& operator<<(Ostream&, const streamLineParticle&)");
502 
503  return os;
504 }
505 
506 
507 // ************************************************************************* //
label size() const
Return the number of elements in the VectorSpace = nCmpt.
Definition: VectorSpaceI.H:67
bool hitPatch(const polyPatch &, trackingData &td, const label patchI, const scalar trackFraction, const tetIndices &tetIs)
Overridable function to handle the particle hitting a patch.
vector point
Point is a vector.
Definition: point.H:41
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:390
dimensioned< scalar > mag(const dimensioned< Type > &)
Wedge front and back plane patch.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
Symmetry patch for non-planar or multi-plane patches.
void hitWedgePatch(const wedgePolyPatch &, trackingData &td)
Overridable function to handle the particle hitting a wedge.
bool keepParticle
Flag to indicate whether to keep particle (false = delete)
Definition: particle.H:112
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
#define forAllIter(Container, container, iter)
Definition: UList.H:440
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
static void writeFields(const Cloud< streamLineParticle > &)
Write.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Particle class that samples fields as it passes through. Used in streamline calculation.
bool move(trackingData &, const scalar trackTime)
Track all particles to their end point.
Ostream & operator<<(Ostream &, const edgeMesh &)
Definition: edgeMeshIO.C:133
Namespace for OpenFOAM.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
virtual bool check(const char *operation) const
Check IOstream status for given operation.
Definition: IOstream.C:92
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
void hitSymmetryPatch(const symmetryPolyPatch &, trackingData &td)
Overridable function to handle the particle hitting a.
virtual bool write() const
Write using setting from DB.
static void readFields(CloudType &c)
Read the fields associated with the owner cloud.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Neighbour processor patch.
static void writeFields(const CloudType &c)
Write the fields associated with the owner cloud.
Class used to pass tracking data to the trackToFace function.
volScalarField & p
Definition: createFields.H:51
const vector & position() const
Return current particle position.
Definition: particleI.H:603
List< DynamicList< vectorList > > & allVectors_
#define forAll(list, i)
Definition: UList.H:421
Base particle class.
Definition: particle.H:78
void hitCyclicPatch(const cyclicPolyPatch &, trackingData &td)
Overridable function to handle the particle hitting a cyclic.
static void readFields(Cloud< streamLineParticle > &)
Read.
Cyclic plane patch.
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
errorManip< error > abort(error &err)
Definition: errorManip.H:131
DynamicList< vectorList > & allPositions_
void hitWallPatch(const wallPolyPatch &, trackingData &td, const tetIndices &)
Overridable function to handle the particle hitting a wallPatch.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
void checkFieldIOobject(const Cloud< ParticleType > &c, const IOField< DataType > &data) const
Check lagrangian data field.
Definition: CloudIO.C:215
List< DynamicList< scalarList > > & allScalars_
streamLineParticle(const polyMesh &c, const vector &position, const label cellI, const label lifeTime)
Construct from components.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
List< vector > vectorList
A List of vectors.
Definition: vectorList.H:50
Foam::wallPolyPatch.
Definition: wallPolyPatch.H:48
label size() const
Definition: Cloud.H:175
bool switchProcessor
Flag to switch processor.
Definition: particle.H:109
error FatalError
void hitSymmetryPlanePatch(const symmetryPlanePolyPatch &, trackingData &td)
Overridable function to handle the particle hitting a.
void hitProcessorPatch(const processorPolyPatch &, trackingData &td)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
const dimensionedScalar c
Speed of light in a vacuum.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:73
U
Definition: pEqn.H:82
IOobject fieldIOobject(const word &fieldName, const IOobject::readOption r) const
Helper to construct IOobject for field and current time.
Definition: CloudIO.C:195
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
A primitive field of type <T> with automated input and output.
Definition: IOField.H:50