Cloud.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 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 "polyTopoChangeMap.H"
31 #include "Time.H"
32 #include "OFstream.H"
33 #include "wallPolyPatch.H"
34 #include "cyclicAMIPolyPatch.H"
36 
37 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
38 
39 template<class ParticleType>
41 {
42  const polyBoundaryMesh& pbm = polyMesh_.boundaryMesh();
43  bool ok = true;
44  forAll(pbm, patchi)
45  {
46  if (isA<cyclicAMIPolyPatch>(pbm[patchi]))
47  {
48  const cyclicAMIPolyPatch& cami =
49  refCast<const cyclicAMIPolyPatch>(pbm[patchi]);
50 
51  ok = ok && cami.singlePatchProc() != -1;
52  }
53  }
54 
55  if (!ok)
56  {
58  << "Particle tracking across AMI patches is only currently "
59  << "supported for cases where the AMI patches reside on a "
60  << "single processor" << abort(FatalError);
61  }
62 }
63 
64 
65 template<class ParticleType>
67 (
68  const polyMesh& pMesh
69 )
70 {
71  const polyBoundaryMesh& pbm = pMesh.boundaryMesh();
72 
73  labelList result(pbm.size(), -1);
74 
75  if (Pstream::parRun())
76  {
77  forAll(pbm, patchi)
78  {
79  if (isA<processorPolyPatch>(pbm[patchi]))
80  {
81  const processorPolyPatch& ppp =
82  refCast<const processorPolyPatch>(pbm[patchi]);
83 
84  result[patchi] = ppp.neighbProcNo();
85  }
86  }
87  }
88 
89  return result;
90 }
91 
92 
93 template<class ParticleType>
95 (
96  const polyMesh& pMesh
97 )
98 {
99  const polyBoundaryMesh& pbm = pMesh.boundaryMesh();
100 
101  labelList result(pbm.size(), -1);
102 
103  if (Pstream::parRun())
104  {
105  PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
106 
107  forAll(pbm, patchi)
108  {
109  if (isA<processorPolyPatch>(pbm[patchi]))
110  {
111  const processorPolyPatch& ppp =
112  refCast<const processorPolyPatch>(pbm[patchi]);
113 
114  UOPstream(ppp.neighbProcNo(), pBufs)()
115  << ppp.index();
116  }
117  }
118 
119  pBufs.finishedSends();
120 
121  forAll(pbm, patchi)
122  {
123  if (isA<processorPolyPatch>(pbm[patchi]))
124  {
125  const processorPolyPatch& ppp =
126  refCast<const processorPolyPatch>(pbm[patchi]);
127 
128  UIPstream(ppp.neighbProcNo(), pBufs)()
129  >> result[patchi];
130  }
131  }
132  }
133 
134  return result;
135 }
136 
137 
138 template<class ParticleType>
140 (
141  const polyMesh& pMesh
142 )
143 {
144  const polyBoundaryMesh& pbm = pMesh.boundaryMesh();
145 
146  labelListList result(pbm.size(), labelList());
147 
148  forAll(pbm, patchi)
149  {
150  if (isA<nonConformalCyclicPolyPatch>(pbm[patchi]))
151  {
152  const nonConformalCyclicPolyPatch& nccPp =
153  refCast<const nonConformalCyclicPolyPatch>(pbm[patchi]);
154 
155  result[nccPp.origPatchID()].append(patchi);
156  }
157  }
158 
159  return result;
160 }
161 
162 
163 template<class ParticleType>
165 {
166  const polyBoundaryMesh& pbm = polyMesh_.boundaryMesh();
167 
168  forAll(patchNonConformalCyclicPatches_, patchi)
169  {
170  forAll(patchNonConformalCyclicPatches_[patchi], i)
171  {
172  const label nccPatchi =
173  patchNonConformalCyclicPatches_[patchi][i];
174 
175  const nonConformalCyclicPolyPatch& nccPp =
176  refCast<const nonConformalCyclicPolyPatch>(pbm[nccPatchi]);
177 
178  if (nccPp.owner()) nccPp.rays();
179  }
180  }
181 }
182 
183 
184 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
185 
186 template<class ParticleType>
188 (
189  const polyMesh& pMesh,
190  const word& cloudName,
191  const IDLList<ParticleType>& particles
192 )
193 :
194  cloud(pMesh, cloudName),
196  polyMesh_(pMesh),
197  patchNbrProc_(patchNbrProc(pMesh)),
198  patchNbrProcPatch_(patchNbrProcPatch(pMesh)),
199  patchNonConformalCyclicPatches_(patchNonConformalCyclicPatches(pMesh)),
200  globalPositionsPtr_()
201 {
202  checkPatches();
203 
204  // Ask for the tetBasePtIs and oldCellCentres to trigger all processors to
205  // build them, otherwise, if some processors have no particles then there
206  // is a comms mismatch.
207  polyMesh_.tetBasePtIs();
208  polyMesh_.oldCellCentres();
209 
210  if (particles.size())
211  {
213  }
214 }
215 
216 
217 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
218 
219 template<class ParticleType>
221 {
222  this->append(pPtr);
223 }
224 
225 
226 template<class ParticleType>
228 {
229  delete(this->remove(&p));
230 }
231 
232 
233 template<class ParticleType>
235 {
236  forAllIter(typename Cloud<ParticleType>, *this, pIter)
237  {
238  ParticleType& p = pIter();
239 
240  if (p.cell() == -1)
241  {
243  << "deleting lost particle at position " << p.position()
244  << endl;
245 
246  deleteParticle(p);
247  }
248  }
249 }
250 
251 
252 template<class ParticleType>
254 {
255  // Reset particle count and particles only
256  // - not changing the cloud object registry or reference to the polyMesh
257  ParticleType::particleCount_ = 0;
259 }
260 
261 
262 template<class ParticleType>
263 template<class TrackCloudType>
265 (
266  TrackCloudType& cloud,
267  typename ParticleType::trackingData& td,
268  const scalar trackTime
269 )
270 {
271  // Clear the global positions as these are about to change
272  globalPositionsPtr_.clear();
273 
274  // Ensure rays are available for non conformal transfers
275  storeRays();
276 
277  // Initialise the stepFraction moved for the particles
278  forAllIter(typename Cloud<ParticleType>, *this, pIter)
279  {
280  pIter().reset(0);
281  }
282 
283  // Create transfer buffers
284  PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
285 
286  // Create lists of particles and patch indices to transfer
287  List<IDLList<ParticleType>> sendParticles(Pstream::nProcs());
288  List<DynamicList<label>> sendPatchIndices(Pstream::nProcs());
289 
290  // While there are particles to transfer
291  while (true)
292  {
293  // Clear the transfer lists
294  forAll(sendParticles, proci)
295  {
296  sendParticles[proci].clear();
297  sendPatchIndices[proci].clear();
298  }
299 
300  // Loop over all particles
301  forAllIter(typename Cloud<ParticleType>, *this, pIter)
302  {
303  ParticleType& p = pIter();
304 
305  // Move the particle
306  bool keepParticle = p.move(cloud, td, trackTime);
307 
308  // If the particle is to be kept
309  if (keepParticle)
310  {
311  if (td.sendToProc != -1)
312  {
313  #ifdef FULLDEBUG
314  if (!Pstream::parRun() || !p.onBoundaryFace())
315  {
317  << "Switch processor flag is true when no parallel "
318  << "transfer is possible. This is a bug."
319  << exit(FatalError);
320  }
321  #endif
322 
323  p.prepareForParallelTransfer(td);
324 
325  sendParticles[td.sendToProc].append(this->remove(&p));
326 
327  sendPatchIndices[td.sendToProc].append(td.sendToPatch);
328  }
329  }
330  else
331  {
332  deleteParticle(p);
333  }
334  }
335 
336  // If running in serial then everything has been moved, so finish
337  if (!Pstream::parRun())
338  {
339  break;
340  }
341 
342  // Clear transfer buffers
343  pBufs.clear();
344 
345  // Stream into send buffers
346  forAll(sendParticles, proci)
347  {
348  if (sendParticles[proci].size())
349  {
350  UOPstream particleStream(proci, pBufs);
351 
352  particleStream
353  << sendPatchIndices[proci]
354  << sendParticles[proci];
355  }
356  }
357 
358  // Start sending. Sets number of bytes transferred.
359  labelList receiveSizes(Pstream::nProcs());
360  pBufs.finishedSends(receiveSizes);
361 
362  // Determine if any particles were transferred. If not, then finish.
363  bool transferred = false;
364  forAll(receiveSizes, proci)
365  {
366  if (receiveSizes[proci])
367  {
368  transferred = true;
369  break;
370  }
371  }
372  reduce(transferred, orOp<bool>());
373  if (!transferred)
374  {
375  break;
376  }
377 
378  // Retrieve from receive buffers and add into the cloud
379  forAll(receiveSizes, proci)
380  {
381  if (receiveSizes[proci])
382  {
383  UIPstream particleStream(proci, pBufs);
384 
385  const labelList receivePatchIndices(particleStream);
386 
387  IDLList<ParticleType> newParticles
388  (
389  particleStream,
390  typename ParticleType::iNew(polyMesh_)
391  );
392 
393  label i = 0;
394 
395  forAllIter(typename Cloud<ParticleType>, newParticles, iter)
396  {
397  const label patchi = receivePatchIndices[i ++];
398 
399  ParticleType& p = iter();
400 
401  td.sendToPatch = patchi;
402 
403  p.correctAfterParallelTransfer(td);
404 
405  addParticle(newParticles.remove(&p));
406  }
407  }
408  }
409  }
410 }
411 
412 
413 template<class ParticleType>
415 {
416  if (!globalPositionsPtr_.valid())
417  {
419  << "Global positions are not available. "
420  << "Cloud::storeGlobalPositions has not been called."
421  << exit(FatalError);
422  }
423 
424  // Ask for the tetBasePtIs to trigger all processors to build
425  // them, otherwise, if some processors have no particles then
426  // there is a comms mismatch.
427  polyMesh_.tetBasePtIs();
428  polyMesh_.oldCellCentres();
429 
430  const vectorField& positions = globalPositionsPtr_();
431 
432  label i = 0;
433  forAllIter(typename Cloud<ParticleType>, *this, iter)
434  {
435  iter().autoMap(positions[i], mapper);
436  ++ i;
437  }
438 }
439 
440 
441 template<class ParticleType>
443 {
444  OFstream pObj
445  (
446  this->db().time().path()/this->name() + "_positions.obj"
447  );
448 
449  forAllConstIter(typename Cloud<ParticleType>, *this, pIter)
450  {
451  const ParticleType& p = pIter();
452  pObj<< "v " << p.position().x() << " " << p.position().y() << " "
453  << p.position().z() << nl;
454  }
455 
456  pObj.flush();
457 }
458 
459 
460 template<class ParticleType>
462 {
463  // Store the global positions for later use by autoMap. It would be
464  // preferable not to need this. If the polyTopoChangeMap object passed to
465  // autoMap had a copy of the old mesh then the global positions could be
466  // recovered within autoMap, and this pre-processing would not be necessary.
467 
468  globalPositionsPtr_.reset(new vectorField(this->size()));
469 
470  vectorField& positions = globalPositionsPtr_();
471 
472  label i = 0;
473  forAllConstIter(typename Cloud<ParticleType>, *this, iter)
474  {
475  positions[i] = iter().position();
476  ++ i;
477  }
478 }
479 
480 
481 // * * * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * //
482 
483 #include "CloudIO.C"
484 
485 // ************************************************************************* //
Template class for intrusive linked lists.
Definition: ILList.H:50
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
void cloudReset(const Cloud< ParticleType > &c)
Reset the particles.
Definition: Cloud.C:253
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
void finishedSends(const bool block=true)
Mark all sends as having been done. This will start receives.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
void deleteParticle(ParticleType &)
Remove particle from cloud and delete.
Definition: Cloud.C:227
const labelList & patchNbrProcPatch() const
Map from patch index to the corresponding patch index on the.
Definition: Cloud.H:178
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
const labelList & patchNbrProc() const
Map from patch index to the neighbouring processor index.
Definition: Cloud.H:171
Output to file stream.
Definition: OFstream.H:82
Cloud(const polyMesh &mesh, const word &cloudName, const IDLList< ParticleType > &particles)
Construct from mesh and a list of particles.
Definition: Cloud.C:188
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
void autoMap(const polyTopoChangeMap &)
Remap the cells of particles corresponding to the.
Definition: Cloud.C:414
volVectorField vectorField(fieldObject, mesh)
const dimensionedScalar c
Speed of light in a vacuum.
Combination-Reduction operation for a parallel run. The information from all nodes is collected on th...
Input inter-processor communications stream operating on external buffer.
Definition: UIPstream.H:53
void addParticle(ParticleType *pPtr)
Transfer particle to cloud.
Definition: Cloud.C:220
const labelListList & patchNonConformalCyclicPatches() const
Return map from patch index to connected non-conformal cyclics.
Definition: Cloud.H:184
void storeGlobalPositions() const
Call this before a topology change. Stores the particles global.
Definition: Cloud.C:461
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:125
A class for handling words, derived from string.
Definition: word.H:59
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
virtual void flush()
Flush stream.
Definition: OSstream.C:207
List< label > labelList
A List of labels.
Definition: labelList.H:56
T * remove(T *p)
Remove and return element.
Definition: UILList.H:142
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Base cloud calls templated on particle type.
Definition: Cloud.H:52
Output inter-processor communications stream operating on external buffer.
Definition: UOPstream.H:54
static const char nl
Definition: Ostream.H:260
Buffers for inter-processor communications streams (UOPstream, UIPstream).
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
void deleteLostParticles()
Remove lost particles from cloud and delete.
Definition: Cloud.C:234
void clear()
Clear storage and reset.
void move(TrackCloudType &cloud, typename ParticleType::trackingData &td, const scalar trackTime)
Move the particles.
Definition: Cloud.C:265
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
volScalarField & p
void writePositions() const
Write positions to <cloudName>_positions.obj file.
Definition: Cloud.C:442