Cloud.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-2014 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 "Cloud.H"
27 #include "processorPolyPatch.H"
28 #include "globalMeshData.H"
30 #include "mapPolyMesh.H"
31 #include "Time.H"
32 #include "OFstream.H"
33 #include "wallPolyPatch.H"
34 #include "cyclicAMIPolyPatch.H"
35 
36 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
37 
38 template<class ParticleType>
40 {
41  const polyBoundaryMesh& pbm = polyMesh_.boundaryMesh();
42  bool ok = true;
43  forAll(pbm, patchI)
44  {
45  if (isA<cyclicAMIPolyPatch>(pbm[patchI]))
46  {
47  const cyclicAMIPolyPatch& cami =
48  refCast<const cyclicAMIPolyPatch>(pbm[patchI]);
49 
50  if (cami.owner())
51  {
52  ok = ok && (cami.AMI().singlePatchProc() != -1);
53  }
54  }
55  }
56 
57  if (!ok)
58  {
59  FatalErrorIn("void Foam::Cloud<ParticleType>::initCloud(const bool)")
60  << "Particle tracking across AMI patches is only currently "
61  << "supported for cases where the AMI patches reside on a "
62  << "single processor" << abort(FatalError);
63  }
64 }
65 
66 
67 template<class ParticleType>
69 {
70  cellWallFacesPtr_.reset(new PackedBoolList(pMesh().nCells(), false));
71 
72  PackedBoolList& cellWallFaces = cellWallFacesPtr_();
73 
74  const polyBoundaryMesh& patches = polyMesh_.boundaryMesh();
75 
76  forAll(patches, patchI)
77  {
78  if (isA<wallPolyPatch>(patches[patchI]))
79  {
80  const polyPatch& patch = patches[patchI];
81 
82  const labelList& pFaceCells = patch.faceCells();
83 
84  forAll(pFaceCells, pFCI)
85  {
86  cellWallFaces[pFaceCells[pFCI]] = true;
87  }
88  }
89  }
90 }
91 
92 
93 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
94 
95 template<class ParticleType>
97 (
98  const polyMesh& pMesh,
99  const IDLList<ParticleType>& particles
100 )
101 :
102  cloud(pMesh),
104  polyMesh_(pMesh),
105  labels_(),
106  nTrackingRescues_(),
107  cellWallFacesPtr_()
108 {
109  checkPatches();
110 
111  // Ask for the tetBasePtIs to trigger all processors to build
112  // them, otherwise, if some processors have no particles then
113  // there is a comms mismatch.
114  polyMesh_.tetBasePtIs();
115 
117 }
118 
119 
120 template<class ParticleType>
122 (
123  const polyMesh& pMesh,
124  const word& cloudName,
125  const IDLList<ParticleType>& particles
126 )
127 :
128  cloud(pMesh, cloudName),
130  polyMesh_(pMesh),
131  labels_(),
132  nTrackingRescues_(),
133  cellWallFacesPtr_()
134 {
135  checkPatches();
136 
137  // Ask for the tetBasePtIs to trigger all processors to build
138  // them, otherwise, if some processors have no particles then
139  // there is a comms mismatch.
140  polyMesh_.tetBasePtIs();
141 
143 }
144 
145 
146 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
147 
148 template<class ParticleType>
150 const
151 {
152  if (!cellWallFacesPtr_.valid())
153  {
154  calcCellWallFaces();
155  }
156 
157  return cellWallFacesPtr_();
158 }
159 
160 
161 template<class ParticleType>
163 {
164  this->append(pPtr);
165 }
166 
167 
168 template<class ParticleType>
170 {
171  delete(this->remove(&p));
172 }
173 
174 
175 template<class ParticleType>
177 {
178  // Reset particle cound and particles only
179  // - not changing the cloud object registry or reference to the polyMesh
180  ParticleType::particleCount_ = 0;
182 }
183 
184 
185 template<class ParticleType>
186 template<class TrackData>
187 void Foam::Cloud<ParticleType>::move(TrackData& td, const scalar trackTime)
188 {
189  const polyBoundaryMesh& pbm = pMesh().boundaryMesh();
190  const globalMeshData& pData = polyMesh_.globalData();
191 
192  // Which patches are processor patches
193  const labelList& procPatches = pData.processorPatches();
194 
195  // Indexing of patches into the procPatches list
196  const labelList& procPatchIndices = pData.processorPatchIndices();
197 
198  // Indexing of equivalent patch on neighbour processor into the
199  // procPatches list on the neighbour
200  const labelList& procPatchNeighbours = pData.processorPatchNeighbours();
201 
202  // Which processors this processor is connected to
203  const labelList& neighbourProcs = pData[Pstream::myProcNo()];
204 
205  // Indexing from the processor number into the neighbourProcs list
206  labelList neighbourProcIndices(Pstream::nProcs(), -1);
207 
208  forAll(neighbourProcs, i)
209  {
210  neighbourProcIndices[neighbourProcs[i]] = i;
211  }
212 
213  // Initialise the stepFraction moved for the particles
214  forAllIter(typename Cloud<ParticleType>, *this, pIter)
215  {
216  pIter().stepFraction() = 0;
217  }
218 
219  // Reset nTrackingRescues
220  nTrackingRescues_ = 0;
221 
222 
223  // List of lists of particles to be transfered for all of the
224  // neighbour processors
225  List<IDLList<ParticleType> > particleTransferLists
226  (
227  neighbourProcs.size()
228  );
229 
230  // List of destination processorPatches indices for all of the
231  // neighbour processors
232  List<DynamicList<label> > patchIndexTransferLists
233  (
234  neighbourProcs.size()
235  );
236 
237  // Allocate transfer buffers
238  PstreamBuffers pBufs(Pstream::nonBlocking);
239 
240 
241  // While there are particles to transfer
242  while (true)
243  {
244  particleTransferLists = IDLList<ParticleType>();
245  forAll(patchIndexTransferLists, i)
246  {
247  patchIndexTransferLists[i].clear();
248  }
249 
250  // Loop over all particles
251  forAllIter(typename Cloud<ParticleType>, *this, pIter)
252  {
253  ParticleType& p = pIter();
254 
255  // Move the particle
256  bool keepParticle = p.move(td, trackTime);
257 
258  // If the particle is to be kept
259  // (i.e. it hasn't passed through an inlet or outlet)
260  if (keepParticle)
261  {
262  // If we are running in parallel and the particle is on a
263  // boundary face
264  if (Pstream::parRun() && p.face() >= pMesh().nInternalFaces())
265  {
266  label patchI = pbm.whichPatch(p.face());
267 
268  // ... and the face is on a processor patch
269  // prepare it for transfer
270  if (procPatchIndices[patchI] != -1)
271  {
272  label n = neighbourProcIndices
273  [
274  refCast<const processorPolyPatch>
275  (
276  pbm[patchI]
277  ).neighbProcNo()
278  ];
279 
280  p.prepareForParallelTransfer(patchI, td);
281 
282  particleTransferLists[n].append(this->remove(&p));
283 
284  patchIndexTransferLists[n].append
285  (
286  procPatchNeighbours[patchI]
287  );
288  }
289  }
290  }
291  else
292  {
293  deleteParticle(p);
294  }
295  }
296 
297  if (!Pstream::parRun())
298  {
299  break;
300  }
301 
302 
303  // Clear transfer buffers
304  pBufs.clear();
305 
306  // Stream into send buffers
307  forAll(particleTransferLists, i)
308  {
309  if (particleTransferLists[i].size())
310  {
311  UOPstream particleStream
312  (
313  neighbourProcs[i],
314  pBufs
315  );
316 
317  particleStream
318  << patchIndexTransferLists[i]
319  << particleTransferLists[i];
320  }
321  }
322 
323 
324  // Start sending. Sets number of bytes transferred
325  labelListList allNTrans(Pstream::nProcs());
326  pBufs.finishedSends(allNTrans);
327 
328 
329  bool transfered = false;
330 
331  forAll(allNTrans, i)
332  {
333  forAll(allNTrans[i], j)
334  {
335  if (allNTrans[i][j])
336  {
337  transfered = true;
338  break;
339  }
340  }
341  }
342 
343  if (!transfered)
344  {
345  break;
346  }
347 
348  // Retrieve from receive buffers
349  forAll(neighbourProcs, i)
350  {
351  label neighbProci = neighbourProcs[i];
352 
353  label nRec = allNTrans[neighbProci][Pstream::myProcNo()];
354 
355  if (nRec)
356  {
357  UIPstream particleStream(neighbProci, pBufs);
358 
359  labelList receivePatchIndex(particleStream);
360 
361  IDLList<ParticleType> newParticles
362  (
363  particleStream,
364  typename ParticleType::iNew(polyMesh_)
365  );
366 
367  label pI = 0;
368 
369  forAllIter(typename Cloud<ParticleType>, newParticles, newpIter)
370  {
371  ParticleType& newp = newpIter();
372 
373  label patchI = procPatches[receivePatchIndex[pI++]];
374 
375  newp.correctAfterParallelTransfer(patchI, td);
376 
377  addParticle(newParticles.remove(&newp));
378  }
379  }
380  }
381  }
382 
383  if (cloud::debug)
384  {
385  reduce(nTrackingRescues_, sumOp<label>());
386 
387  if (nTrackingRescues_ > 0)
388  {
389  Info<< nTrackingRescues_ << " tracking rescue corrections" << endl;
390  }
391  }
392 }
393 
394 
395 template<class ParticleType>
396 template<class TrackData>
398 (
399  TrackData& td,
400  const mapPolyMesh& mapper
401 )
402 {
403  if (cloud::debug)
404  {
405  Info<< "Cloud<ParticleType>::autoMap(TrackData&, const mapPolyMesh&) "
406  << "for lagrangian cloud " << cloud::name() << endl;
407  }
408 
409  const labelList& reverseCellMap = mapper.reverseCellMap();
410  const labelList& reverseFaceMap = mapper.reverseFaceMap();
411 
412  // Reset stored data that relies on the mesh
413 // polyMesh_.clearCellTree();
414  cellWallFacesPtr_.clear();
415 
416  // Ask for the tetBasePtIs to trigger all processors to build
417  // them, otherwise, if some processors have no particles then
418  // there is a comms mismatch.
419  polyMesh_.tetBasePtIs();
420 
421 
422  forAllIter(typename Cloud<ParticleType>, *this, pIter)
423  {
424  ParticleType& p = pIter();
425 
426  if (reverseCellMap[p.cell()] >= 0)
427  {
428  p.cell() = reverseCellMap[p.cell()];
429 
430  if (p.face() >= 0 && reverseFaceMap[p.face()] >= 0)
431  {
432  p.face() = reverseFaceMap[p.face()];
433  }
434  else
435  {
436  p.face() = -1;
437  }
438 
439  p.initCellFacePt();
440  }
441  else
442  {
443  label trackStartCell = mapper.mergedCell(p.cell());
444 
445  if (trackStartCell < 0)
446  {
447  trackStartCell = 0;
448  p.cell() = 0;
449  }
450  else
451  {
452  p.cell() = trackStartCell;
453  }
454 
455  vector pos = p.position();
456 
457  const_cast<vector&>(p.position()) =
458  polyMesh_.cellCentres()[trackStartCell];
459 
460  p.stepFraction() = 0;
461 
462  p.initCellFacePt();
463 
464  p.track(pos, td);
465  }
466  }
467 }
468 
469 
470 template<class ParticleType>
472 {
473  OFstream pObj
474  (
475  this->db().time().path()/this->name() + "_positions.obj"
476  );
477 
478  forAllConstIter(typename Cloud<ParticleType>, *this, pIter)
479  {
480  const ParticleType& p = pIter();
481  pObj<< "v " << p.position().x() << " " << p.position().y() << " "
482  << p.position().z() << nl;
483  }
484 
485  pObj.flush();
486 }
487 
488 
489 // * * * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * //
490 
491 #include "CloudIO.C"
492 
493 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Output to file stream.
Definition: OFstream.H:81
Buffers for inter-processor communications streams (UOPstream, UIPstream).
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:497
const labelList & processorPatchIndices() const
Return list of indices into processorPatches_ for each patch.
void autoMap(TrackData &td, const mapPolyMesh &)
Remap the cells of particles corresponding to the.
Definition: Cloud.C:398
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Output inter-processor communications stream operating on external buffer.
Definition: UOPstream.H:54
Input inter-processor communications stream operating on external buffer.
Definition: UIPstream.H:53
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
#define forAllIter(Container, container, iter)
Definition: UList.H:440
A bit-packed bool list.
A class for handling words, derived from string.
Definition: word.H:59
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
messageStream Info
T * remove(T *p)
Remove and return element.
Definition: UILList.H:139
void move(TrackData &td, const scalar trackTime)
Move the particles.
Definition: Cloud.C:187
const labelList & processorPatchNeighbours() const
Return processorPatchIndices of the neighbours.
const PackedBoolList & cellHasWallFaces() const
Whether each cell has any wall faces (demand driven data)
Definition: Cloud.C:149
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
label n
void clear()
Clear the list, i.e. set size to zero.
Definition: List.C:379
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual void flush()
Flush stream.
Definition: OSstream.C:246
Foam::polyBoundaryMesh.
volScalarField & p
Definition: createFields.H:51
#define forAll(list, i)
Definition: UList.H:421
label mergedCell(const label oldCellI) const
If cell is removed return cell (on new mesh) it merged into.
Definition: mapPolyMesh.H:534
void deleteParticle(ParticleType &)
Remove particle from cloud and delete.
Definition: Cloud.C:169
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
A cloud is a collection of lagrangian particles.
Definition: cloud.H:51
errorManip< error > abort(error &err)
Definition: errorManip.H:131
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
Various mesh related information for a parallel run. Upon construction, constructs all info using par...
label nInternalFaces() const
void cloudReset(const Cloud< ParticleType > &c)
Reset the particles.
Definition: Cloud.C:176
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
const dimensionedScalar c
Speed of light in a vacuum.
const labelList & reverseCellMap() const
Reverse cell map.
Definition: mapPolyMesh.H:528
Base cloud calls templated on particle type.
Definition: Cloud.H:52
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:97
dimensionedScalar pos(const dimensionedScalar &ds)
Cloud(const polyMesh &mesh, const IDLList< ParticleType > &particles)
Construct from mesh and a list of particles.
Definition: Cloud.C:97
void addParticle(ParticleType *pPtr)
Transfer particle to cloud.
Definition: Cloud.C:162
const labelList & processorPatches() const
Return list of processor patch labels.
Combination-Reduction operation for a parallel run. The information from all nodes is collected on th...
void writePositions() const
Write positions to <cloudName>_positions.obj file.
Definition: Cloud.C:471