applications/utilities: Reorganized

This commit is contained in:
Henry Weller
2016-06-13 17:02:43 +01:00
parent fdf69b3c25
commit 029825b83c
25 changed files with 1 additions and 3 deletions

View File

@ -0,0 +1,3 @@
particleTracks.C
EXE = $(FOAM_APPBIN)/particleTracks

View File

@ -0,0 +1,12 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lfileFormats \
-lgenericPatchFields \
-llagrangian

View File

@ -0,0 +1,23 @@
IOdictionary propsDict
(
IOobject
(
"particleTrackProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED
)
);
const word cloudName(propsDict.lookup("cloudName"));
label sampleFrequency(readLabel(propsDict.lookup("sampleFrequency")));
// outputMode: compositeFile, filePerTrack
//word outputmode(propsDict.lookup("outputMode"))
label maxPositions(readLabel(propsDict.lookup("maxPositions")));
// outputFormat: raw, vtk
//word outputFormat(propsDict.lookup("outputFormat"));
word setFormat(propsDict.lookupOrDefault<word>("setFormat", "vtk"));

View File

@ -0,0 +1,25 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object particleTrackProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
cloudName reactingCloud1;
sampleFrequency 1;
maxPositions 1000000;
setFormat vtk;
// ************************************************************************* //

View File

@ -0,0 +1,254 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
particleTracks
Description
Generates a VTK file of particle tracks for cases that were computed using
a tracked-parcel-type cloud.
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Cloud.H"
#include "IOdictionary.H"
#include "fvMesh.H"
#include "Time.H"
#include "timeSelector.H"
#include "OFstream.H"
#include "passiveParticleCloud.H"
#include "writer.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
timeSelector::addOptions();
#include "addRegionOption.H"
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createNamedMesh.H"
#include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
fileName vtkPath(runTime.path()/"VTK");
mkDir(vtkPath);
Info<< "Scanning times to determine track data for cloud " << cloudName
<< nl << endl;
labelList maxIds(Pstream::nProcs(), -1);
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Time = " << runTime.timeName() << endl;
Info<< " Reading particle positions" << endl;
passiveParticleCloud myCloud(mesh, cloudName);
Info<< " Read " << returnReduce(myCloud.size(), sumOp<label>())
<< " particles" << endl;
forAllConstIter(passiveParticleCloud, myCloud, iter)
{
label origId = iter().origId();
label origProc = iter().origProc();
if (origProc >= maxIds.size())
{
// Expand size
maxIds.setSize(origProc+1, -1);
}
maxIds[origProc] = max(maxIds[origProc], origId);
}
}
label maxNProcs = returnReduce(maxIds.size(), maxOp<label>());
Info<< "Detected particles originating from " << maxNProcs
<< " processors." << nl << endl;
maxIds.setSize(maxNProcs, -1);
Pstream::listCombineGather(maxIds, maxEqOp<label>());
Pstream::listCombineScatter(maxIds);
labelList numIds = maxIds + 1;
Info<< nl << "Particle statistics:" << endl;
forAll(maxIds, proci)
{
Info<< " Found " << numIds[proci] << " particles originating"
<< " from processor " << proci << endl;
}
Info<< nl << endl;
// calc starting ids for particles on each processor
List<label> startIds(numIds.size(), 0);
for (label i = 0; i < numIds.size()-1; i++)
{
startIds[i+1] += startIds[i] + numIds[i];
}
label nParticle = startIds.last() + numIds[startIds.size()-1];
// number of tracks to generate
label nTracks = nParticle/sampleFrequency;
// storage for all particle tracks
List<DynamicList<vector>> allTracks(nTracks);
Info<< "\nGenerating " << nTracks << " particle tracks for cloud "
<< cloudName << nl << endl;
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Time = " << runTime.timeName() << endl;
List<pointField> allPositions(Pstream::nProcs());
List<labelField> allOrigIds(Pstream::nProcs());
List<labelField> allOrigProcs(Pstream::nProcs());
// Read particles. Will be size 0 if no particles.
Info<< " Reading particle positions" << endl;
passiveParticleCloud myCloud(mesh, cloudName);
// collect the track data on all processors that have positions
allPositions[Pstream::myProcNo()].setSize
(
myCloud.size(),
point::zero
);
allOrigIds[Pstream::myProcNo()].setSize(myCloud.size(), 0);
allOrigProcs[Pstream::myProcNo()].setSize(myCloud.size(), 0);
label i = 0;
forAllConstIter(passiveParticleCloud, myCloud, iter)
{
allPositions[Pstream::myProcNo()][i] = iter().position();
allOrigIds[Pstream::myProcNo()][i] = iter().origId();
allOrigProcs[Pstream::myProcNo()][i] = iter().origProc();
i++;
}
// collect the track data on the master processor
Pstream::gatherList(allPositions);
Pstream::gatherList(allOrigIds);
Pstream::gatherList(allOrigProcs);
Info<< " Constructing tracks" << nl << endl;
if (Pstream::master())
{
forAll(allPositions, proci)
{
forAll(allPositions[proci], i)
{
label globalId =
startIds[allOrigProcs[proci][i]]
+ allOrigIds[proci][i];
if (globalId % sampleFrequency == 0)
{
label trackId = globalId/sampleFrequency;
if (allTracks[trackId].size() < maxPositions)
{
allTracks[trackId].append
(
allPositions[proci][i]
);
}
}
}
}
}
}
if (Pstream::master())
{
PtrList<coordSet> tracks(allTracks.size());
forAll(allTracks, trackI)
{
tracks.set
(
trackI,
new coordSet
(
"track" + Foam::name(trackI),
"distance"
)
);
tracks[trackI].transfer(allTracks[trackI]);
}
autoPtr<writer<scalar>> scalarFormatterPtr = writer<scalar>::New
(
setFormat
);
//OFstream vtkTracks(vtkPath/"particleTracks.vtk");
fileName vtkFile
(
scalarFormatterPtr().getFileName
(
tracks[0],
wordList(0)
)
);
OFstream vtkTracks
(
vtkPath
/ "particleTracks." + vtkFile.ext()
);
Info<< "\nWriting particle tracks in " << setFormat
<< " format to " << vtkTracks.name()
<< nl << endl;
scalarFormatterPtr().write
(
true, // writeTracks
tracks,
wordList(0),
List<List<scalarField>>(0),
vtkTracks
);
}
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,3 @@
steadyParticleTracks.C
EXE = $(FOAM_APPBIN)/steadyParticleTracks

View File

@ -0,0 +1,9 @@
EXE_INC = \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-llagrangian \
-lmeshTools \
-lfiniteVolume

View File

@ -0,0 +1,9 @@
const word dictName("particleTrackDict");
#include "setConstantMeshDictionaryIO.H"
IOdictionary propsDict(dictIO);
word cloudName(propsDict.lookup("cloudName"));
List<word> userFields(propsDict.lookup("fields"));

View File

@ -0,0 +1,21 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object particleTrackDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
cloudName reactingCloud1Tracks;
fields ( d U T );
// ************************************************************************* //

View File

@ -0,0 +1,340 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
steadyParticleTracks
Description
Generates a VTK file of particle tracks for cases that were computed using
a steady-state cloud
NOTE: case must be re-constructed (if running in parallel) before use
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Cloud.H"
#include "IOdictionary.H"
#include "fvMesh.H"
#include "Time.H"
#include "timeSelector.H"
#include "OFstream.H"
#include "passiveParticleCloud.H"
#include "SortableList.H"
#include "IOobjectList.H"
#include "PtrList.H"
#include "Field.H"
#include "steadyParticleTracksTemplates.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
using namespace Foam;
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
label validateFields
(
const List<word>& userFields,
const IOobjectList& cloudObjs
)
{
List<bool> ok(userFields.size(), false);
forAll(userFields, i)
{
ok[i] = ok[i] || fieldOk<label>(cloudObjs, userFields[i]);
ok[i] = ok[i] || fieldOk<scalar>(cloudObjs, userFields[i]);
ok[i] = ok[i] || fieldOk<vector>(cloudObjs, userFields[i]);
ok[i] = ok[i] || fieldOk<sphericalTensor>(cloudObjs, userFields[i]);
ok[i] = ok[i] || fieldOk<symmTensor>(cloudObjs, userFields[i]);
ok[i] = ok[i] || fieldOk<tensor>(cloudObjs, userFields[i]);
}
label nOk = 0;
forAll(ok, i)
{
if (ok[i])
{
nOk++;
}
else
{
Info << "\n*** Warning: user specified field '" << userFields[i]
<< "' unavailable" << endl;
}
}
return nOk;
}
template<>
void writeVTK(OFstream& os, const label& value)
{
os << value;
}
template<>
void writeVTK(OFstream& os, const scalar& value)
{
os << value;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::noParallel();
timeSelector::addOptions();
#include "addRegionOption.H"
#include "addDictOption.H"
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createNamedMesh.H"
#include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
fileName vtkPath(runTime.path()/"VTK");
mkDir(vtkPath);
typedef HashTable<label, labelPair, labelPair::Hash<>> trackTableType;
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Time = " << runTime.timeName() << endl;
fileName vtkTimePath(runTime.path()/"VTK"/runTime.timeName());
mkDir(vtkTimePath);
Info<< " Reading particle positions" << endl;
PtrList<passiveParticle> particles(0);
// transfer particles to (more convenient) list
{
passiveParticleCloud ppc(mesh, cloudName);
Info<< "\n Read " << returnReduce(ppc.size(), sumOp<label>())
<< " particles" << endl;
particles.setSize(ppc.size());
label i = 0;
forAllIter(passiveParticleCloud, ppc, iter)
{
particles.set(i++, ppc.remove(&iter()));
}
// myCloud should now be empty
}
List<label> particleToTrack(particles.size());
label nTracks = 0;
{
trackTableType trackTable;
forAll(particles, i)
{
const label origProc = particles[i].origProc();
const label origId = particles[i].origId();
const trackTableType::const_iterator& iter =
trackTable.find(labelPair(origProc, origId));
if (iter == trackTable.end())
{
particleToTrack[i] = nTracks;
trackTable.insert(labelPair(origProc, origId), nTracks);
nTracks++;
}
else
{
particleToTrack[i] = iter();
}
}
}
if (nTracks == 0)
{
Info<< "\n No track data" << endl;
}
else
{
Info<< "\n Generating " << nTracks << " tracks" << endl;
// determine length of each track
labelList trackLengths(nTracks, 0);
forAll(particleToTrack, i)
{
const label trackI = particleToTrack[i];
trackLengths[trackI]++;
}
// particle "age" property used to sort the tracks
List<SortableList<scalar>> agePerTrack(nTracks);
List<List<label>> particleMap(nTracks);
forAll(trackLengths, i)
{
const label length = trackLengths[i];
agePerTrack[i].setSize(length);
particleMap[i].setSize(length);
}
// store the particle age per track
IOobjectList cloudObjs
(
mesh,
runTime.timeName(),
cloud::prefix/cloudName
);
// TODO: gather age across all procs
{
tmp<scalarField> tage =
readParticleField<scalar>("age", cloudObjs);
const scalarField& age = tage();
List<label> trackSamples(nTracks, 0);
forAll(particleToTrack, i)
{
const label trackI = particleToTrack[i];
const label sampleI = trackSamples[trackI];
agePerTrack[trackI][sampleI] = age[i];
particleMap[trackI][sampleI] = i;
trackSamples[trackI]++;
}
tage.clear();
}
if (Pstream::master())
{
OFstream os(vtkTimePath/"particleTracks.vtk");
Info<< "\n Writing particle tracks to " << os.name() << endl;
label nPoints = sum(trackLengths);
os << "# vtk DataFile Version 2.0" << nl
<< "particleTracks" << nl
<< "ASCII" << nl
<< "DATASET POLYDATA" << nl
<< "POINTS " << nPoints << " float" << nl;
Info<< "\n Writing points" << endl;
{
forAll(agePerTrack, i)
{
agePerTrack[i].sort();
const labelList& ids = agePerTrack[i].indices();
labelList& particleIds = particleMap[i];
{
// update addressing
List<label> sortedIds(ids);
forAll(sortedIds, j)
{
sortedIds[j] = particleIds[ids[j]];
}
particleIds = sortedIds;
}
forAll(ids, j)
{
const label localId = particleIds[j];
const vector& pos = particles[localId].position();
os << pos.x() << ' ' << pos.y() << ' ' << pos.z()
<< nl;
}
}
}
// write track (line) connectivity to file
Info<< "\n Writing track lines" << endl;
os << "\nLINES " << nTracks << ' ' << nPoints + nTracks << nl;
// Write ids of track points to file
{
label globalPtI = 0;
forAll(particleMap, i)
{
os << particleMap[i].size() << nl;
forAll(particleMap[i], j)
{
os << ' ' << globalPtI++;
if (((j + 1) % 10 == 0) && (j != 0))
{
os << nl;
}
}
os << nl;
}
}
const label nFields = validateFields(userFields, cloudObjs);
os << "POINT_DATA " << nPoints << nl
<< "FIELD attributes " << nFields << nl;
Info<< "\n Processing fields" << nl << endl;
processFields<label>(os, particleMap, userFields, cloudObjs);
processFields<scalar>(os, particleMap, userFields, cloudObjs);
processFields<vector>(os, particleMap, userFields, cloudObjs);
processFields<sphericalTensor>
(os, particleMap, userFields, cloudObjs);
processFields<symmTensor>
(os, particleMap, userFields, cloudObjs);
processFields<tensor>(os, particleMap, userFields, cloudObjs);
}
}
Info<< endl;
}
Info<< "\ndone" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,190 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "steadyParticleTracksTemplates.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
bool fieldOk(const IOobjectList& cloudObjs, const word& name)
{
IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
return (objects.lookup(name) != NULL);
}
template<class Type>
tmp<Field<Type>> readParticleField
(
const word& name,
const IOobjectList cloudObjs
)
{
IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
const IOobject* obj = objects.lookup(name);
if (obj != NULL)
{
IOField<Type> newField(*obj);
return tmp<Field<Type>>(new Field<Type>(newField.xfer()));
}
FatalErrorInFunction
<< "error: cloud field name " << name << " not found"
<< abort(FatalError);
return Field<Type>::null();
}
template<class Type>
void readFields
(
PtrList<List<Type>>& values,
const List<word>& fieldNames,
const IOobjectList& cloudObjs
)
{
IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
forAll(fieldNames, j)
{
const IOobject* obj = objects.lookup(fieldNames[j]);
if (obj != NULL)
{
Info<< " reading field " << fieldNames[j] << endl;
IOField<Type> newField(*obj);
values.set(j, new List<Type>(newField.xfer()));
}
else
{
FatalErrorInFunction
<< "Unable to read field " << fieldNames[j]
<< abort(FatalError);
}
}
}
template<class Type>
void writeVTK(OFstream& os, const Type& value)
{
os << value.component(0);
for (label i=1; i<pTraits<Type>::nComponents; i++)
{
os << ' ' << value.component(i);
}
}
template<class Type>
void writeVTKFields
(
OFstream& os,
const PtrList<List<Type>>& values,
const List<List<label>>& addr,
const List<word>& fieldNames
)
{
label step = max(floor(8/pTraits<Type>::nComponents), 1);
forAll(values, fieldi)
{
Info<< " writing field " << fieldNames[fieldi] << endl;
os << nl << fieldNames[fieldi] << ' ' << pTraits<Type>::nComponents
<< ' ' << values[fieldi].size() << " float" << nl;
label offset = 0;
forAll(addr, trackI)
{
const List<label> ids(addr[trackI]);
List<Type> data(UIndirectList<Type>(values[fieldi], ids));
label nData = data.size() - 1;
forAll(data, i)
{
writeVTK<Type>(os, data[i]);
if (((i + 1) % step == 0) || (i == nData))
{
os << nl;
}
else
{
os << ' ';
}
}
offset += ids.size();
}
}
}
template<class Type>
void processFields
(
OFstream& os,
const List<List<label>>& addr,
const List<word>& userFieldNames,
const IOobjectList& cloudObjs
)
{
IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
if (objects.size())
{
DynamicList<word> fieldNames(objects.size());
forAll(userFieldNames, i)
{
IOobject* obj = objects.lookup(userFieldNames[i]);
if (obj != NULL)
{
fieldNames.append(obj->name());
}
}
fieldNames.shrink();
PtrList<List<Type>> values(fieldNames.size());
readFields<Type>(values, fieldNames, cloudObjs);
writeVTKFields<Type>
(
os,
values,
addr,
fieldNames
);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,93 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#ifndef steadyParticleTracksTemplates_H
#define steadyParticleTracksTemplates_H
#include "OFstream.H"
#include "SortableList.H"
#include "IOobjectList.H"
#include "PtrList.H"
#include "Field.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
bool fieldOk(const IOobjectList& cloudObjs, const word& name);
template<class Type>
tmp<Field<Type>> readParticleField
(
const word& name,
const IOobjectList cloudObjs
);
template<class Type>
void readFields
(
PtrList<List<Type>>& values,
const List<word>& fields,
const IOobjectList& cloudObjs
);
template<class Type>
void writeVTK(OFstream& os, const Type& value);
template<class Type>
void writeVTKFields
(
OFstream& os,
const PtrList<List<Type>>& values,
const List<SortableList<scalar>>& agePerTrack,
const List<word>& fieldNames
);
void processFields
(
OFstream& os,
const List<SortableList<scalar>>& agePerTrack,
const List<word>& userFieldNames,
const IOobjectList& cloudObjs
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "steadyParticleTracksTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,3 @@
streamFunction.C
EXE = $(FOAM_APPBIN)/streamFunction

View File

@ -0,0 +1,8 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lgenericPatchFields \
-lfiniteVolume \
-lmeshTools

View File

@ -0,0 +1,491 @@
/*---------------------------------------------------------------------------* \
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
streamFunction
Description
Calculates and writes the stream function of velocity field U at each
time.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "pointFields.H"
#include "emptyPolyPatch.H"
#include "symmetryPlanePolyPatch.H"
#include "symmetryPolyPatch.H"
#include "wedgePolyPatch.H"
#include "OSspecific.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
timeSelector::addOptions();
#include "addRegionOption.H"
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createNamedMesh.H"
label nD = mesh.nGeometricD();
if (nD != 2)
{
FatalErrorInFunction
<< "Case is not 2D, stream-function cannot be computed"
<< exit(FatalError);
}
Vector<label> slabNormal((Vector<label>::one - mesh.geometricD())/2);
const direction slabDir
(
slabNormal
& Vector<label>(Vector<label>::X, Vector<label>::Y, Vector<label>::Z)
);
scalar thickness = vector(slabNormal) & mesh.bounds().span();
const pointMesh& pMesh = pointMesh::New(mesh);
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
Info<< nl << "Time: " << runTime.timeName() << endl;
IOobject phiHeader
(
"phi",
runTime.timeName(),
mesh,
IOobject::NO_READ
);
if (phiHeader.headerOk())
{
mesh.readUpdate();
Info<< nl << "Reading field phi" << endl;
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
);
pointScalarField streamFunction
(
IOobject
(
"streamFunction",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
pMesh,
dimensionedScalar("zero", phi.dimensions(), 0.0)
);
labelList visitedPoint(mesh.nPoints());
forAll(visitedPoint, pointi)
{
visitedPoint[pointi] = 0;
}
label nVisited = 0;
label nVisitedOld = 0;
const faceUList& faces = mesh.faces();
const pointField& points = mesh.points();
label nInternalFaces = mesh.nInternalFaces();
vectorField unitAreas(mesh.faceAreas());
unitAreas /= mag(unitAreas);
const polyPatchList& patches = mesh.boundaryMesh();
bool finished = true;
// Find the boundary face with zero flux. set the stream function
// to zero on that face
bool found = false;
do
{
found = false;
forAll(patches, patchi)
{
const primitivePatch& bouFaces = patches[patchi];
if (!isType<emptyPolyPatch>(patches[patchi]))
{
forAll(bouFaces, facei)
{
if
(
magSqr(phi.boundaryField()[patchi][facei])
< SMALL
)
{
const labelList& zeroPoints = bouFaces[facei];
// Zero flux face found
found = true;
forAll(zeroPoints, pointi)
{
if (visitedPoint[zeroPoints[pointi]] == 1)
{
found = false;
break;
}
}
if (found)
{
Info<< "Zero face: patch: " << patchi
<< " face: " << facei << endl;
forAll(zeroPoints, pointi)
{
streamFunction[zeroPoints[pointi]] = 0;
visitedPoint[zeroPoints[pointi]] = 1;
nVisited++;
}
break;
}
}
}
}
if (found) break;
}
if (!found)
{
Info<< "zero flux boundary face not found. "
<< "Using cell as a reference."
<< endl;
const cellList& c = mesh.cells();
forAll(c, cI)
{
labelList zeroPoints = c[cI].labels(mesh.faces());
bool found = true;
forAll(zeroPoints, pointi)
{
if (visitedPoint[zeroPoints[pointi]] == 1)
{
found = false;
break;
}
}
if (found)
{
forAll(zeroPoints, pointi)
{
streamFunction[zeroPoints[pointi]] = 0.0;
visitedPoint[zeroPoints[pointi]] = 1;
nVisited++;
}
break;
}
else
{
FatalErrorInFunction
<< "Cannot find initialisation face or a cell."
<< abort(FatalError);
}
}
}
// Loop through all faces. If one of the points on
// the face has the streamfunction value different
// from -1, all points with -1 ont that face have the
// streamfunction value equal to the face flux in
// that point plus the value in the visited point
do
{
finished = true;
for
(
label facei = nInternalFaces;
facei<faces.size();
facei++
)
{
const labelList& curBPoints = faces[facei];
bool bPointFound = false;
scalar currentBStream = 0.0;
vector currentBStreamPoint(0, 0, 0);
forAll(curBPoints, pointi)
{
// Check if the point has been visited
if (visitedPoint[curBPoints[pointi]] == 1)
{
// The point has been visited
currentBStream =
streamFunction[curBPoints[pointi]];
currentBStreamPoint =
points[curBPoints[pointi]];
bPointFound = true;
break;
}
}
if (bPointFound)
{
// Sort out other points on the face
forAll(curBPoints, pointi)
{
// Check if the point has been visited
if (visitedPoint[curBPoints[pointi]] == 0)
{
label patchNo =
mesh.boundaryMesh().whichPatch(facei);
if
(
!isType<emptyPolyPatch>
(patches[patchNo])
&& !isType<symmetryPlanePolyPatch>
(patches[patchNo])
&& !isType<symmetryPolyPatch>
(patches[patchNo])
&& !isType<wedgePolyPatch>
(patches[patchNo])
)
{
label faceNo =
mesh.boundaryMesh()[patchNo]
.whichFace(facei);
vector edgeHat =
points[curBPoints[pointi]]
- currentBStreamPoint;
edgeHat.replace(slabDir, 0);
edgeHat /= mag(edgeHat);
vector nHat = unitAreas[facei];
if (edgeHat.y() > VSMALL)
{
visitedPoint[curBPoints[pointi]] =
1;
nVisited++;
streamFunction[curBPoints[pointi]]
=
currentBStream
+ phi.boundaryField()
[patchNo][faceNo]
*sign(nHat.x());
}
else if (edgeHat.y() < -VSMALL)
{
visitedPoint[curBPoints[pointi]] =
1;
nVisited++;
streamFunction[curBPoints[pointi]]
=
currentBStream
- phi.boundaryField()
[patchNo][faceNo]
*sign(nHat.x());
}
else
{
if (edgeHat.x() > VSMALL)
{
visitedPoint
[curBPoints[pointi]] = 1;
nVisited++;
streamFunction
[curBPoints[pointi]] =
currentBStream
+ phi.boundaryField()
[patchNo][faceNo]
*sign(nHat.y());
}
else if (edgeHat.x() < -VSMALL)
{
visitedPoint
[curBPoints[pointi]] = 1;
nVisited++;
streamFunction
[curBPoints[pointi]] =
currentBStream
- phi.boundaryField()
[patchNo][faceNo]
*sign(nHat.y());
}
}
}
}
}
}
else
{
finished = false;
}
}
for (label facei=0; facei<nInternalFaces; facei++)
{
// Get the list of point labels for the face
const labelList& curPoints = faces[facei];
bool pointFound = false;
scalar currentStream = 0.0;
point currentStreamPoint(0, 0, 0);
forAll(curPoints, pointi)
{
// Check if the point has been visited
if (visitedPoint[curPoints[pointi]] == 1)
{
// The point has been visited
currentStream =
streamFunction[curPoints[pointi]];
currentStreamPoint =
points[curPoints[pointi]];
pointFound = true;
break;
}
}
if (pointFound)
{
// Sort out other points on the face
forAll(curPoints, pointi)
{
// Check if the point has been visited
if (visitedPoint[curPoints[pointi]] == 0)
{
vector edgeHat =
points[curPoints[pointi]]
- currentStreamPoint;
edgeHat.replace(slabDir, 0);
edgeHat /= mag(edgeHat);
vector nHat = unitAreas[facei];
if (edgeHat.y() > VSMALL)
{
visitedPoint[curPoints[pointi]] = 1;
nVisited++;
streamFunction[curPoints[pointi]] =
currentStream
+ phi[facei]*sign(nHat.x());
}
else if (edgeHat.y() < -VSMALL)
{
visitedPoint[curPoints[pointi]] = 1;
nVisited++;
streamFunction[curPoints[pointi]] =
currentStream
- phi[facei]*sign(nHat.x());
}
}
}
}
else
{
finished = false;
}
}
Info<< ".";
if (nVisited == nVisitedOld)
{
// Find new seed. This must be a
// multiply connected domain
Info<< nl << "Exhausted a seed. Looking for new seed "
<< "(this is correct for multiply connected "
<< "domains).";
break;
}
else
{
nVisitedOld = nVisited;
}
} while (!finished);
Info<< endl;
} while (!finished);
// Normalise the stream-function by the 2D mesh thickness
streamFunction /= thickness;
streamFunction.boundaryFieldRef() = 0.0;
streamFunction.write();
}
else
{
WarningInFunction
<< "Flux field does not exist."
<< " Stream function not calculated" << endl;
}
}
Info<< "\nEnd\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,3 @@
wallHeatFlux.C
EXE = $(FOAM_APPBIN)/wallHeatFlux

View File

@ -0,0 +1,23 @@
EXE_INC = \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-lreactionThermophysicalModels \
-lgenericPatchFields \
-lspecie \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lsolidThermo \
-lfiniteVolume \
-lfvOptions \
-lmeshTools

View File

@ -0,0 +1,72 @@
autoPtr<basicThermo> thermo
(
basicThermo::New(mesh)
);
const volScalarField& h = thermo->he();
// Register copy of thermo density
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh
),
thermo->rho()
);
// Construct turbulence model (if fluid)
autoPtr<volVectorField> UPtr;
autoPtr<surfaceScalarField> phiPtr;
autoPtr<compressible::turbulenceModel> turbulence;
if (isA<fluidThermo>(thermo()))
{
UPtr.reset
(
new volVectorField
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
)
);
const volVectorField& U = UPtr();
#include "compressibleCreatePhi.H"
// Copy phi to autoPtr. Rename to make sure copy is now registered as 'phi'.
phi.rename("phiFluid");
phiPtr.reset(new surfaceScalarField("phi", phi));
turbulence = compressible::turbulenceModel::New
(
rho,
U,
phiPtr(),
refCast<const fluidThermo>(thermo())
);
}
// Read radiative heat-flux if available
volScalarField Qr
(
IOobject
(
"Qr",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
);

View File

@ -0,0 +1,157 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
wallHeatFlux
Description
Calculates and writes the heat flux for all patches as the boundary field
of a volScalarField and also prints the integrated flux for all wall
patches.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "turbulentFluidThermoModel.H"
#include "solidThermo.H"
#include "wallFvPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
timeSelector::addOptions();
#include "addRegionOption.H"
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createNamedMesh.H"
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Time = " << runTime.timeName() << endl;
mesh.readUpdate();
#include "createFields.H"
surfaceScalarField heatFlux
(
fvc::interpolate
(
(
turbulence.valid()
? turbulence->alphaEff()()
: thermo->alpha()
)
)*fvc::snGrad(h)
);
const surfaceScalarField::Boundary& patchHeatFlux =
heatFlux.boundaryField();
const volScalarField::Boundary& patchRadHeatFlux =
Qr.boundaryField();
const surfaceScalarField::Boundary& magSf =
mesh.magSf().boundaryField();
Info<< "\nWall heat fluxes [W]" << endl;
forAll(patchHeatFlux, patchi)
{
if (isA<wallFvPatch>(mesh.boundary()[patchi]))
{
scalar convFlux = gSum(magSf[patchi]*patchHeatFlux[patchi]);
scalar radFlux = -gSum(magSf[patchi]*patchRadHeatFlux[patchi]);
Info<< mesh.boundary()[patchi].name() << endl
<< " convective: " << convFlux << endl
<< " radiative: " << radFlux << endl
<< " total: " << convFlux + radFlux << endl;
}
}
Info<< endl;
volScalarField wallHeatFlux
(
IOobject
(
"wallHeatFlux",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("wallHeatFlux", heatFlux.dimensions(), 0.0)
);
volScalarField::Boundary& wallHeatFluxBf =
wallHeatFlux.boundaryFieldRef();
forAll(wallHeatFluxBf, patchi)
{
wallHeatFluxBf[patchi] = patchHeatFlux[patchi];
}
wallHeatFlux.write();
// Write the total heat-flux including the radiative contribution
// if available
if (Qr.headerOk())
{
volScalarField totalWallHeatFlux
(
IOobject
(
"totalWallHeatFlux",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar
(
"totalWallHeatFlux",
heatFlux.dimensions(),
0.0
)
);
volScalarField::Boundary& totalWallHeatFluxBf =
totalWallHeatFlux.boundaryFieldRef();
forAll(totalWallHeatFluxBf, patchi)
{
totalWallHeatFluxBf[patchi] =
patchHeatFlux[patchi] - patchRadHeatFlux[patchi];
}
totalWallHeatFlux.write();
}
}
Info<< "End" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,3 @@
wdot.C
EXE = $(FOAM_APPBIN)/wdot

View File

@ -0,0 +1,8 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lgenericPatchFields \
-lfiniteVolume \
-lmeshTools

View File

@ -0,0 +1,123 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
wdot
Description
Calculates and writes wdot for each time.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
timeSelector::addOptions();
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createMesh.H"
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
mesh.readUpdate();
volScalarField mgb
(
IOobject
(
"mgb",
runTime.timeName(),
mesh,
IOobject::MUST_READ
),
mesh
);
volScalarField Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh,
IOobject::MUST_READ
),
mesh
);
volScalarField Xi
(
IOobject
(
"Xi",
runTime.timeName(),
mesh,
IOobject::MUST_READ
),
mesh
);
volScalarField St
(
IOobject
(
"St",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
Xi*Su
);
St.write();
volScalarField wdot
(
IOobject
(
"wdot",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
St*mgb
);
wdot.write();
}
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,3 @@
writeCellCentres.C
EXE = $(FOAM_APPBIN)/writeCellCentres

View File

@ -0,0 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lfiniteVolume

View File

@ -0,0 +1,140 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
Write the three components of the cell centres as volScalarFields so
they can be used in postprocessing in thresholding.
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "timeSelector.H"
#include "Time.H"
#include "fvMesh.H"
#include "vectorIOField.H"
#include "volFields.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
timeSelector::addOptions();
#include "addRegionOption.H"
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createNamedMesh.H"
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Time = " << runTime.timeName() << endl;
// Check for new mesh
mesh.readUpdate();
//volVectorField cc
//(
// IOobject
// (
// "cellCentres",
// runTime.timeName(),
// mesh,
// IOobject::NO_READ,
// IOobject::AUTO_WRITE
// ),
// 1.0*mesh.C()
//);
//
//Info<< "Writing cellCentre positions to " << cc.name() << " in "
// << runTime.timeName() << endl;
//cc.write();
Info<< "Writing components of cellCentre positions to volScalarFields"
<< " ccx, ccy, ccz in " << runTime.timeName() << endl;
for (direction i=0; i<vector::nComponents; i++)
{
volScalarField cci
(
IOobject
(
"cc" + word(vector::componentNames[i]),
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh.C().component(i)
);
cci.write();
}
volScalarField V
(
IOobject
(
"V",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
mesh,
dimensionedScalar("V", mesh.V().dimensions(), 0.0),
calculatedFvPatchField<scalar>::typeName
);
V.ref() = mesh.V();
volScalarField::Boundary& Vbf = V.boundaryFieldRef();
forAll(Vbf, patchi)
{
Vbf[patchi] = Vbf[patchi].patch().magSf();
}
Info<< "Writing cellVolumes and patch faceAreas to " << V.name()
<< " in " << runTime.timeName() << endl;
V.write();
}
Info<< "\nEnd\n" << endl;
return 0;
}
// ************************************************************************* //