particleTracks.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2021-2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 particleTracks
29
30Group
31 grpPostProcessingUtilities
32
33Description
34 Generate particle tracks for cases that were computed using a
35 tracked-parcel-type cloud.
36
37\*---------------------------------------------------------------------------*/
38
39#include "argList.H"
40#include "Cloud.H"
41#include "IOdictionary.H"
42#include "fvMesh.H"
43#include "Time.H"
44#include "timeSelector.H"
45#include "coordSetWriter.H"
48
49// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
50
51using namespace Foam;
52
53template<class Type>
54void writeTrackFields
55(
57 HashTable<List<DynamicList<Type>>>& fieldTable
58)
59{
60 for (const word& fieldName : fieldTable.sortedToc())
61 {
62 // Steal and reshape from List<DynamicList> to List<Field>
63 auto& input = fieldTable[fieldName];
64
65 List<Field<Type>> fields(input.size());
66 forAll(input, tracki)
67 {
68 fields[tracki].transfer(input[tracki]);
69 }
70
71 writer.write(fieldName, fields);
72 }
73}
74
75
76// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
77
78int main(int argc, char *argv[])
79{
80 argList::addNote
81 (
82 "Generate a file of particle tracks for cases that were"
83 " computed using a tracked-parcel-type cloud"
84 );
85 timeSelector::addOptions();
86 #include "addRegionOption.H"
87
88 // Less frequently used - reduce some clutter
89 argList::setAdvanced("decomposeParDict");
90 argList::setAdvanced("noFunctionObjects");
91
92 argList::addOption
93 (
94 "dict",
95 "file",
96 "Alternative particleTracksProperties dictionary"
97 );
98 argList::addOption
99 (
100 "stride",
101 "int",
102 "Override the sample-frequency"
103 );
104 argList::addOption
105 (
106 "fields",
107 "wordRes",
108 "Specify single or multiple fields to write "
109 "(default: all or 'fields' from dictionary)\n"
110 "Eg, 'T' or '( \"U.*\" )'"
111 );
112 argList::addOption
113 (
114 "format",
115 "name",
116 "The writer format "
117 "(default: vtk or 'setFormat' from dictionary)"
118 );
119 argList::addVerboseOption("Additional verbosity");
120
121 #include "setRootCase.H"
122 #include "createTime.H"
123 instantList timeDirs = timeSelector::select0(runTime, args);
124 #include "createNamedMesh.H"
125
126 // ------------------------------------------------------------------------
127 // Control properties
128
129 #include "createControls.H"
130
132 args.readIfPresent("format", setFormat);
133
135 sampleFrequency = max(1, sampleFrequency); // sanity
136
137 // Setup the writer
138 auto writerPtr =
139 coordSetWriter::New
140 (
141 setFormat,
142 formatOptions.optionalSubDict(setFormat, keyType::LITERAL)
143 );
144
145 writerPtr().useTracks(true);
146
147 if (args.verbose())
148 {
149 writerPtr().verbose(true);
150 }
151
152 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
153
154 particleTracksSampler trackSampler;
155
156 Info<< "Scanning times to determine track data for cloud " << cloudName
157 << nl << endl;
158
159 {
160 labelList maxIds(Pstream::nProcs(), -1);
161 forAll(timeDirs, timei)
162 {
163 runTime.setTime(timeDirs[timei], timei);
164 Info<< "Time = " << runTime.timeName() << endl;
165
167
168 Info<< " Read " << returnReduce(myCloud.size(), sumOp<label>())
169 << " particles" << endl;
170
171 for (const passiveParticle& p : myCloud)
172 {
173 const label origId = p.origId();
174 const label origProc = p.origProc();
175
176 // Handle case where processing particle data generated in
177 // parallel using more cores than for this application.
178 if (origProc >= maxIds.size())
179 {
180 maxIds.resize(origProc+1, -1);
181 }
182
183 maxIds[origProc] = max(maxIds[origProc], origId);
184 }
185 }
186
187 const label maxNProcs = returnReduce(maxIds.size(), maxOp<label>());
188 maxIds.resize(maxNProcs, -1);
189
190 Pstream::listCombineAllGather(maxIds, maxEqOp<label>());
191
192 // From ids to count
193 const labelList numIds = maxIds + 1;
194
195 // Set parcel addressing
196 trackSampler.reset(numIds);
197
198 Info<< nl
199 << "Detected particles originating from "
200 << maxNProcs << " processors." << nl
201 << "Particle statistics:" << endl;
202
203 if (Pstream::master())
204 {
205 const globalIndex& parcelAddr = trackSampler.parcelAddr();
206
207 for (const label proci : parcelAddr.allProcs())
208 {
209 Info<< " Found " << parcelAddr.localSize(proci)
210 << " particles originating"
211 << " from processor " << proci << nl;
212 }
213 }
214 }
215
217
218
219 // Number of tracks to generate
220 const label nTracks = trackSampler.nTracks();
221
222
223 // Storage for all particle tracks
224 List<DynamicList<point>> allTrackPos(nTracks);
225 List<DynamicList<scalar>> allTrackTimes(nTracks);
226
227 // Track field values by name/type
228 HashTable<List<DynamicList<label>>> labelFields; // <- mostly unused
231 HashTable<List<DynamicList<sphericalTensor>>> sphericalTensorFields;
234
235 Info<< nl
236 << "Generating " << nTracks
237 << " particle tracks for cloud " << cloudName << nl << endl;
238
239 forAll(timeDirs, timei)
240 {
241 runTime.setTime(timeDirs[timei], timei);
242 Info<< "Time = " << runTime.timeName() << " (processing)" << endl;
243
244 // Read particles. Will be size 0 if no particles.
246
247 pointField localPositions(myCloud.size());
248 scalarField localTimes(myCloud.size(), runTime.value());
249
250 // Gather track data from all processors that have positions
251 trackSampler.resetCloud(myCloud.size());
252 {
253 labelField& origIds = trackSampler.origParcelIds_;
254 labelField& origProcs = trackSampler.origProcIds_;
255
256 label np = 0;
257 for (const passiveParticle& p : myCloud)
258 {
259 origIds[np] = p.origId();
260 origProcs[np] = p.origProc();
261 localPositions[np] = p.position();
262 ++np;
263 }
264
265 trackSampler.gatherInplace(origIds);
266 trackSampler.gatherInplace(origProcs);
267 }
268
269
270 // Read cloud fields (from disk) into object registry
272 (
274 (
275 "cloudFields",
276 runTime.timeName(),
277 runTime
278 )
279 );
280
281 myCloud.readFromFiles(obr, acceptFields, excludeFields);
282
283 // Create track positions and track time fields
284 // (not registered as IOFields)
285
286 trackSampler.createTrackField(localPositions, allTrackPos);
287 trackSampler.createTrackField(localTimes, allTrackTimes);
288
289 // Create the track fields
291 trackSampler.setTrackFields(obr, scalarFields);
292 trackSampler.setTrackFields(obr, vectorFields);
293 trackSampler.setTrackFields(obr, sphericalTensorFields);
294 trackSampler.setTrackFields(obr, symmTensorFields);
295 trackSampler.setTrackFields(obr, tensorFields);
296 }
297
298 const label nFields =
299 (
300 labelFields.size()
301 + scalarFields.size()
302 + vectorFields.size()
303 + sphericalTensorFields.size()
304 + symmTensorFields.size()
305 + tensorFields.size()
306 );
307
308 Info<< nl
309 << "Extracted " << nFields << " cloud fields" << nl;
310
311 if (nFields)
312 {
313 #undef doLocalCode
314 #define doLocalCode(FieldContent) \
315 if (!FieldContent.empty()) \
316 { \
317 Info<< " "; \
318 for (const word& fieldName : FieldContent.sortedToc()) \
319 { \
320 Info<< ' ' << fieldName; \
321 } \
322 Info<< nl; \
323 }
324
325 doLocalCode(labelFields);
326 doLocalCode(scalarFields);
327 doLocalCode(vectorFields);
328 doLocalCode(sphericalTensorFields);
329 doLocalCode(symmTensorFields);
330 doLocalCode(tensorFields);
331 #undef doLocalCode
332 }
333
334 Info<< nl
335 << "Writing particle tracks (" << setFormat << " format)" << nl;
336
337 if (Pstream::master())
338 {
339 PtrList<coordSet> tracks(allTrackPos.size());
340 List<scalarField> times(allTrackPos.size());
341 forAll(tracks, tracki)
342 {
343 tracks.set
344 (
345 tracki,
346 new coordSet("track" + Foam::name(tracki), "xyz")
347 );
348 tracks[tracki].transfer(allTrackPos[tracki]);
349 times[tracki].transfer(allTrackTimes[tracki]);
350
351 if (!tracki) tracks[0].rename("tracks");
352 }
353
366
367
368 // Write tracks with fields
369 if (nFields)
370 {
371 auto& writer = *writerPtr;
372
373 const fileName outputPath
374 (
375 functionObject::outputPrefix/cloud::prefix/cloudName
376 / "particleTracks" / "tracks"
377 );
378
379 Info<< " "
380 << runTime.relativePath(outputPath) << endl;
381
382 writer.open(tracks, outputPath);
383 writer.setTrackTimes(times);
384 writer.nFields(nFields);
385
387 writeTrackFields(writer, scalarFields);
388 writeTrackFields(writer, vectorFields);
389 writeTrackFields(writer, sphericalTensorFields);
390 writeTrackFields(writer, symmTensorFields);
391 writeTrackFields(writer, tensorFields);
392 }
393 else
394 {
395 Info<< "Warning: no fields, did not write" << endl;
396 }
397 }
398
399 Info<< nl << "End\n" << endl;
400
401 return 0;
402}
403
404
405// ************************************************************************* //
Y[inertIndex] max(0.0)
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
A HashTable similar to std::unordered_map.
Definition: HashTable.H:123
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
int verbose() const noexcept
Return the verbose flag.
Definition: argListI.H:128
bool readListIfPresent(const word &optName, List< T > &list) const
Definition: argListI.H:394
bool readIfPresent(const word &optName, T &val) const
Read a value from the named option if present.
Definition: argListI.H:323
Base class for writing coordSet(s) and tracks with fields.
Holds list of sampling positions.
Definition: coordSet.H:56
A class for handling file names.
Definition: fileName.H:76
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
label localSize() const
My local size.
Definition: globalIndexI.H:207
labelRange allProcs() const noexcept
Range of process indices for all addressed offsets (processes)
Definition: globalIndexI.H:151
Registry of regIOobjects.
Helper class when generating particle tracks. The interface is fairly rudimentary.
labelField origProcIds_
The originating processor ids.
void resetCloud(const label localCloudSize)
const globalIndex & parcelAddr() const noexcept
The original parcel addressing.
void reset(const labelUList &origParcelCounts)
Define the orig parcel mappings.
label setSampleRate(const label sampleFreq, const label maxPositions, const label maxTracks=-1)
Set the sampling stride, upper limits.
label nTracks() const noexcept
Number of tracks to generate.
void createTrackField(const UList< Type > &values, List< DynamicList< Type > > &trackValues) const
void gatherInplace(List< Type > &fld) const
label setTrackFields(const objectRegistry &obr, HashTable< List< DynamicList< Type > > > &fieldTable) const
labelField origParcelIds_
The originating parcel ids.
A Cloud of passive particles.
Copy of base particle.
A wordRe is a Foam::word, but can contain a regular expression for matching words or strings.
Definition: wordRe.H:83
A class for handling words, derived from Foam::string.
Definition: word.H:68
volScalarField & p
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Foam::argList args(argc, argv)
multivariateSurfaceInterpolationScheme< scalar >::fieldTable fields
Definition: createFields.H:97
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
label maxTracks(propsDict.getOrDefault< label >("maxTracks", -1))
label maxPositions(propsDict.get< label >("maxPositions"))
wordRes excludeFields
label sampleFrequency(propsDict.get< label >("sampleFrequency"))
word setFormat(propsDict.getOrDefault< word >("setFormat", "vtk"))
wordRes acceptFields
const dictionary formatOptions(propsDict.subOrEmptyDict("formatOptions", keyType::LITERAL))
const word cloudName(propsDict.get< word >("cloud"))
#define doLocalCode(GeoField)