All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
steadyParticleTracks.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-2021 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 Application
25  steadyParticleTracks
26 
27 Description
28  Generates a VTK file of particle tracks for cases that were computed using
29  a steady-state cloud
30  NOTE: case must be re-constructed (if running in parallel) before use
31 
32 \*---------------------------------------------------------------------------*/
33 
34 #include "argList.H"
35 #include "Cloud.H"
36 #include "IOdictionary.H"
37 #include "fvMesh.H"
38 #include "Time.H"
39 #include "timeSelector.H"
40 #include "OFstream.H"
41 #include "passiveParticleCloud.H"
42 #include "systemDict.H"
43 
44 #include "SortableList.H"
45 #include "IOobjectList.H"
46 #include "PtrList.H"
47 #include "Field.H"
49 
50 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
51 
52 using namespace Foam;
53 
54 namespace Foam
55 {
56 
57 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
58 
59 label validateFields
60 (
61  const List<word>& userFields,
62  const IOobjectList& cloudObjs
63 )
64 {
65  List<bool> ok(userFields.size(), false);
66 
67  forAll(userFields, i)
68  {
69  ok[i] = ok[i] || fieldOk<label>(cloudObjs, userFields[i]);
70  ok[i] = ok[i] || fieldOk<scalar>(cloudObjs, userFields[i]);
71  ok[i] = ok[i] || fieldOk<vector>(cloudObjs, userFields[i]);
72  ok[i] = ok[i] || fieldOk<sphericalTensor>(cloudObjs, userFields[i]);
73  ok[i] = ok[i] || fieldOk<symmTensor>(cloudObjs, userFields[i]);
74  ok[i] = ok[i] || fieldOk<tensor>(cloudObjs, userFields[i]);
75  }
76 
77  label nOk = 0;
78  forAll(ok, i)
79  {
80  if (ok[i])
81  {
82  nOk++;
83  }
84  else
85  {
86  Info << "\n*** Warning: user specified field '" << userFields[i]
87  << "' unavailable" << endl;
88  }
89  }
90 
91  return nOk;
92 }
93 
94 
95 template<>
96 void writeVTK(OFstream& os, const label& value)
97 {
98  os << value;
99 }
100 
101 
102 template<>
103 void writeVTK(OFstream& os, const scalar& value)
104 {
105  os << value;
106 }
107 
108 }
109 
110 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
111 
112 int main(int argc, char *argv[])
113 {
116  #include "addRegionOption.H"
117  #include "addDictOption.H"
118 
119  #include "setRootCase.H"
120 
121  #include "createTime.H"
123  #include "createNamedMesh.H"
124  #include "createFields.H"
125 
126  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
127 
128  const fileName vtkPath(runTime.rootPath()/runTime.globalCaseName()/"VTK");
129  mkDir(vtkPath);
130 
131  typedef HashTable<label, labelPair, labelPair::Hash<>> trackTableType;
132 
133  forAll(timeDirs, timeI)
134  {
135  runTime.setTime(timeDirs[timeI], timeI);
136  Info<< "Time = " << runTime.timeName() << endl;
137 
138  const fileName vtkTimePath(vtkPath/runTime.timeName());
139  mkDir(vtkTimePath);
140 
141  Info<< " Reading particle positions" << endl;
142 
143  PtrList<passiveParticle> particles(0);
144 
145  // Transfer particles to (more convenient) list
146  {
148  Info<< "\n Read " << returnReduce(ppc.size(), sumOp<label>())
149  << " particles" << endl;
150 
151  particles.setSize(ppc.size());
152 
153  label i = 0;
154  forAllIter(passiveParticleCloud, ppc, iter)
155  {
156  particles.set(i++, ppc.remove(&iter()));
157  }
158 
159  // myCloud should now be empty
160  }
161 
162  List<label> particleToTrack(particles.size());
163  label nTracks = 0;
164 
165  {
166  trackTableType trackTable;
167  forAll(particles, i)
168  {
169  const label origProc = particles[i].origProc();
170  const label origId = particles[i].origId();
171 
172  const trackTableType::const_iterator& iter =
173  trackTable.find(labelPair(origProc, origId));
174 
175  if (iter == trackTable.end())
176  {
177  particleToTrack[i] = nTracks;
178  trackTable.insert(labelPair(origProc, origId), nTracks);
179  nTracks++;
180  }
181  else
182  {
183  particleToTrack[i] = iter();
184  }
185  }
186  }
187 
188 
189  if (nTracks == 0)
190  {
191  Info<< "\n No track data" << endl;
192  }
193  else
194  {
195  Info<< "\n Generating " << nTracks << " tracks" << endl;
196 
197  // Determine length of each track
198  labelList trackLengths(nTracks, 0);
199  forAll(particleToTrack, i)
200  {
201  const label trackI = particleToTrack[i];
202  trackLengths[trackI]++;
203  }
204 
205  // Particle "age" property used to sort the tracks
206  List<SortableList<scalar>> agePerTrack(nTracks);
207  List<List<label>> particleMap(nTracks);
208 
209  forAll(trackLengths, i)
210  {
211  const label length = trackLengths[i];
212  agePerTrack[i].setSize(length);
213  particleMap[i].setSize(length);
214  }
215 
216  // Store the particle age per track
217  IOobjectList cloudObjs
218  (
219  mesh,
220  runTime.timeName(),
222  );
223 
224  // TODO: gather age across all procs
225  {
226  tmp<scalarField> tage =
227  readParticleField<scalar>("age", cloudObjs);
228 
229  const scalarField& age = tage();
230 
231  List<label> trackSamples(nTracks, 0);
232 
233  forAll(particleToTrack, i)
234  {
235  const label trackI = particleToTrack[i];
236  const label sampleI = trackSamples[trackI];
237  agePerTrack[trackI][sampleI] = age[i];
238  particleMap[trackI][sampleI] = i;
239  trackSamples[trackI]++;
240  }
241  tage.clear();
242  }
243 
244 
245  if (Pstream::master())
246  {
247  OFstream os(vtkTimePath/"particleTracks.vtk");
248 
249  Info<< "\n Writing particle tracks to " << os.name() << endl;
250 
251  label nPoints = sum(trackLengths);
252 
253  os << "# vtk DataFile Version 2.0" << nl
254  << "particleTracks" << nl
255  << "ASCII" << nl
256  << "DATASET POLYDATA" << nl
257  << "POINTS " << nPoints << " float" << nl;
258 
259  Info<< "\n Writing points" << endl;
260 
261  {
262  forAll(agePerTrack, i)
263  {
264  agePerTrack[i].sort();
265 
266  const labelList& ids = agePerTrack[i].indices();
267  labelList& particleIds = particleMap[i];
268 
269  {
270  // Update addressing
271  List<label> sortedIds(ids);
272  forAll(sortedIds, j)
273  {
274  sortedIds[j] = particleIds[ids[j]];
275  }
276  particleIds = sortedIds;
277  }
278 
279  forAll(ids, j)
280  {
281  const label localId = particleIds[j];
282  const vector& pos = particles[localId].position();
283  os << pos.x() << ' ' << pos.y() << ' ' << pos.z()
284  << nl;
285  }
286  }
287  }
288 
289 
290  // Write track (line) connectivity to file
291 
292  Info<< "\n Writing track lines" << endl;
293  os << "\nLINES " << nTracks << ' ' << nPoints + nTracks << nl;
294 
295  // Write ids of track points to file
296  {
297  label globalPtI = 0;
298  forAll(particleMap, i)
299  {
300  os << particleMap[i].size() << nl;
301 
302  forAll(particleMap[i], j)
303  {
304  os << ' ' << globalPtI++;
305 
306  if (((j + 1) % 10 == 0) && (j != 0))
307  {
308  os << nl;
309  }
310  }
311 
312  os << nl;
313  }
314  }
315 
316 
317  const label nFields = validateFields(userFields, cloudObjs);
318 
319  os << "POINT_DATA " << nPoints << nl
320  << "FIELD attributes " << nFields << nl;
321 
322  Info<< "\n Processing fields" << nl << endl;
323 
324  processFields<label>(os, particleMap, userFields, cloudObjs);
325  processFields<scalar>(os, particleMap, userFields, cloudObjs);
326  processFields<vector>(os, particleMap, userFields, cloudObjs);
327  processFields<sphericalTensor>
328  (os, particleMap, userFields, cloudObjs);
329  processFields<symmTensor>
330  (os, particleMap, userFields, cloudObjs);
331  processFields<tensor>(os, particleMap, userFields, cloudObjs);
332  }
333  }
334  Info<< endl;
335  }
336 
337  Info<< "End\n" << endl;
338 
339  return 0;
340 }
341 
342 
343 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
A class for handling file names.
Definition: fileName.H:79
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:237
List of IOobjects with searching and retrieving facilities.
Definition: IOobjectList.H:50
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
Output to file stream.
Definition: OFstream.H:82
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
static void noParallel()
Remove the parallel options.
Definition: argList.C:175
const Cmpt & z() const
Definition: VectorI.H:87
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:636
A Cloud of passive particles.
const Cmpt & y() const
Definition: VectorI.H:81
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
dynamicFvMesh & mesh
List< word > userFields(propsDict.lookup("fields"))
const fileName & rootPath() const
Return root path.
Definition: TimePaths.H:95
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
virtual void setTime(const Time &)
Reset the time and time-index to those of the given time.
Definition: Time.C:899
An STL-conforming hash table.
Definition: HashTable.H:61
const Cmpt & x() const
Definition: VectorI.H:75
const fileName & globalCaseName() const
Return global case name.
Definition: TimePaths.H:101
static const char nl
Definition: Ostream.H:260
static const word prefix
The prefix to local: lagrangian.
Definition: cloud.H:62
static instantList select0(Time &runTime, const argList &args)
Return the set of times selected based on the argList options.
Definition: timeSelector.C:252
bool mkDir(const fileName &, mode_t=0777)
Make a directory and return an error if it could not be created.
Definition: POSIX.C:290
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
messageStream Info
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const word cloudName(propsDict.lookup("cloudName"))
void writeVTK(OFstream &os, const Type &value)
A class for managing temporary objects.
Definition: PtrList.H:53
Foam::argList args(argc, argv)
static void addOptions(const bool constant=true, const bool withZero=false)
Add the options handled by timeSelector to argList::validOptions.
Definition: timeSelector.C:114
Namespace for OpenFOAM.