ENH: Added new particle collector cloud function object

This commit is contained in:
andy
2012-11-14 15:20:53 +00:00
parent fe18dc7e7f
commit 1bf00a47b8
4 changed files with 957 additions and 1 deletions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,6 +29,7 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "FacePostProcessing.H"
#include "ParticleCollector.H"
#include "ParticleErosion.H"
#include "ParticleTracks.H"
#include "ParticleTrap.H"
@ -42,6 +43,7 @@ License
makeCloudFunctionObject(CloudType); \
\
makeCloudFunctionObjectType(FacePostProcessing, CloudType); \
makeCloudFunctionObjectType(ParticleCollector, CloudType); \
makeCloudFunctionObjectType(ParticleErosion, CloudType); \
makeCloudFunctionObjectType(ParticleTracks, CloudType); \
makeCloudFunctionObjectType(ParticleTrap, CloudType); \

View File

@ -0,0 +1,668 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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 "ParticleCollector.H"
#include "Pstream.H"
#include "surfaceWriter.H"
#include "unitConversion.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class CloudType>
void Foam::ParticleCollector<CloudType>::makeLogFile
(
const faceList& faces,
const Field<point>& points,
const Field<scalar>& area
)
{
// Create the output file if not already created
if (log_)
{
if (debug)
{
Info<< "Creating output file" << endl;
}
if (Pstream::master())
{
const fileName logDir = outputDir_/this->owner().time().timeName();
// Create directory if does not exist
mkDir(logDir);
// Open new file at start up
outputFilePtr_.reset
(
new OFstream(logDir/(type() + ".dat"))
);
outputFilePtr_()
<< "# Source : " << type() << nl
<< "# Total area : " << sum(area) << nl
<< "# Time";
forAll(faces, i)
{
word id = Foam::name(i);
outputFilePtr_()
<< tab << "area[" << id << "]"
<< tab << "mass[" << id << "]"
<< tab << "massFlowRate[" << id << "]"
<< endl;
}
}
}
}
template<class CloudType>
void Foam::ParticleCollector<CloudType>::initPolygons()
{
mode_ = mtPolygon;
List<Field<point> > polygons(this->coeffDict().lookup("polygons"));
label nPoints = 0;
forAll(polygons, polyI)
{
label np = polygons[polyI].size();
if (np < 3)
{
FatalIOErrorIn
(
"Foam::ParticleCollector<CloudType>::initPolygons()",
this->coeffDict()
)
<< "polygons must consist of at least 3 points"
<< exit(FatalIOError);
}
nPoints += np;
}
label pointOffset = 0;
points_.setSize(nPoints);
faces_.setSize(polygons.size());
faceTris_.setSize(polygons.size());
area_.setSize(polygons.size());
forAll(faces_, faceI)
{
const Field<point>& polyPoints = polygons[faceI];
face f(identity(polyPoints.size()) + pointOffset);
UIndirectList<point>(points_, f) = polyPoints;
area_[faceI] = f.mag(points_);
DynamicList<face> tris;
f.triangles(points_, tris);
faceTris_[faceI].transfer(tris);
faces_[faceI].transfer(f);
pointOffset += polyPoints.size();
}
}
template<class CloudType>
void Foam::ParticleCollector<CloudType>::initConcentricCircles()
{
mode_ = mtConcentricCircle;
vector origin(this->coeffDict().lookup("origin"));
radius_ = this->coeffDict().lookup("radius");
nSector_ = readLabel(this->coeffDict().lookup("nSector"));
label nS = nSector_;
vector refDir;
if (nSector_ > 1)
{
refDir = this->coeffDict().lookup("refDir");
refDir -= normal_*(normal_ & refDir);
refDir /= mag(refDir);
}
else
{
// set 4 quadrants for single sector cases
nS = 4;
vector tangent = vector::zero;
scalar magTangent = 0.0;
Random rnd(1234);
while (magTangent < SMALL)
{
vector v = rnd.vector01();
tangent = v - (v & normal_)*normal_;
magTangent = mag(tangent);
}
refDir = tangent/magTangent;
}
scalar dTheta = 5.0;
scalar dThetaSector = 360.0/scalar(nS);
label intervalPerSector = max(1, ceil(dThetaSector/dTheta));
dTheta = dThetaSector/scalar(intervalPerSector);
label nPointPerSector = intervalPerSector + 1;
label nPointPerRadius = nS*(nPointPerSector - 1);
label nPoint = radius_.size()*nPointPerRadius;
label nFace = radius_.size()*nS;
// add origin
nPoint++;
points_.setSize(nPoint);
faces_.setSize(nFace);
area_.setSize(nFace);
coordSys_ = cylindricalCS("coordSys", origin, normal_, refDir, false);
List<label> ptIDs(identity(nPointPerRadius));
points_[0] = origin;
// points
forAll(radius_, radI)
{
label pointOffset = radI*nPointPerRadius + 1;
for (label i = 0; i < nPointPerRadius; i++)
{
label pI = i + pointOffset;
point pCyl(radius_[radI], degToRad(i*dTheta), 0.0);
points_[pI] = coordSys_.globalPosition(pCyl);
}
}
// faces
DynamicList<label> facePts(2*nPointPerSector);
forAll(radius_, radI)
{
if (radI == 0)
{
for (label secI = 0; secI < nS; secI++)
{
facePts.clear();
// append origin point
facePts.append(0);
for (label ptI = 0; ptI < nPointPerSector; ptI++)
{
label i = ptI + secI*(nPointPerSector - 1);
label id = ptIDs.fcIndex(i - 1) + 1;
facePts.append(id);
}
label faceI = secI + radI*nS;
faces_[faceI] = face(facePts);
area_[faceI] = faces_[faceI].mag(points_);
}
}
else
{
for (label secI = 0; secI < nS; secI++)
{
facePts.clear();
label offset = (radI - 1)*nPointPerRadius + 1;
for (label ptI = 0; ptI < nPointPerSector; ptI++)
{
label i = ptI + secI*(nPointPerSector - 1);
label id = offset + ptIDs.fcIndex(i - 1);
facePts.append(id);
}
for (label ptI = nPointPerSector-1; ptI >= 0; ptI--)
{
label i = ptI + secI*(nPointPerSector - 1);
label id = offset + nPointPerRadius + ptIDs.fcIndex(i - 1);
facePts.append(id);
}
label faceI = secI + radI*nS;
faces_[faceI] = face(facePts);
area_[faceI] = faces_[faceI].mag(points_);
}
}
}
}
template<class CloudType>
Foam::label Foam::ParticleCollector<CloudType>::collectParcelPolygon
(
const point& position,
const vector& U
) const
{
scalar dt = this->owner().db().time().deltaTValue();
point end(position + dt*U);
label dummyNearType = -1;
label dummyNearLabel = -1;
forAll(faces_, faceI)
{
const label facePoint0 = faces_[faceI][0];
const point p0 = points_[facePoint0];
const scalar d1 = normal_ & (position - p0);
const scalar d2 = normal_ & (end - p0);
if (sign(d1) == sign(d2))
{
// did not cross polygon plane
continue;
}
// intersection point
point pIntersect = position + (d1/(d1 - d2))*dt*U;
const List<face>& tris = faceTris_[faceI];
// identify if point is within poly bounds
forAll(tris, triI)
{
const face& tri = tris[triI];
triPointRef t
(
points_[tri[0]],
points_[tri[1]],
points_[tri[2]]
);
if (t.classify(pIntersect, dummyNearType, dummyNearLabel))
{
return faceI;
}
}
}
return -1;
}
template<class CloudType>
Foam::label Foam::ParticleCollector<CloudType>::collectParcelConcentricCircles
(
const point& position,
const vector& U
) const
{
label secI = -1;
scalar dt = this->owner().db().time().deltaTValue();
point end(position + dt*U);
const scalar d1 = normal_ & (position - coordSys_.origin());
const scalar d2 = normal_ & (end - coordSys_.origin());
if (sign(d1) == sign(d2))
{
// did not cross plane
return secI;
}
// intersection point in cylindrical co-ordinate system
point pCyl = coordSys_.localPosition(position + (d1/(d1 - d2))*dt*U);
scalar r = pCyl[0];
if (r < radius_.last())
{
label radI = 0;
while (r > radius_[radI])
{
radI++;
}
if (nSector_ == 1)
{
secI = 4*radI;
}
else
{
scalar theta = pCyl[1] + constant::mathematical::pi;
secI =
nSector_*radI
+ floor
(
scalar(nSector_)*theta/constant::mathematical::twoPi
);
}
}
return secI;
}
template<class CloudType>
void Foam::ParticleCollector<CloudType>::write()
{
const fvMesh& mesh = this->owner().mesh();
const Time& time = mesh.time();
scalar timeNew = time.value();
scalar timeElapsed = timeNew - timeOld_;
totalTime_ += timeElapsed;
const scalar alpha = (totalTime_ - timeElapsed)/totalTime_;
const scalar beta = timeElapsed/totalTime_;
forAll(faces_, faceI)
{
massFlowRate_[faceI] =
alpha*massFlowRate_[faceI] + beta*mass_[faceI]/timeElapsed;
massTotal_[faceI] += mass_[faceI];
}
const label procI = Pstream::myProcNo();
Info<< type() << " output:" << nl;
if (outputFilePtr_.valid())
{
outputFilePtr_() << time.timeName();
}
Field<scalar> faceMassTotal(mass_.size());
Field<scalar> faceMassFlowRate(massFlowRate_.size());
forAll(faces_, faceI)
{
scalarList allProcMass(Pstream::nProcs());
allProcMass[procI] = massTotal_[faceI];
Pstream::gatherList(allProcMass);
faceMassTotal[faceI] = sum(allProcMass);
scalarList allProcMassFlowRate(Pstream::nProcs());
allProcMassFlowRate[procI] = massFlowRate_[faceI];
Pstream::gatherList(allProcMassFlowRate);
faceMassFlowRate[faceI] = sum(allProcMassFlowRate);
Info<< " face " << faceI
<< ": total mass = " << faceMassTotal[faceI]
<< "; average mass flow rate = " << faceMassFlowRate[faceI]
<< nl;
if (outputFilePtr_.valid())
{
outputFilePtr_()
<< tab << area_[faceI]
<< tab << faceMassTotal[faceI]
<< tab << faceMassFlowRate[faceI]
<< endl;
}
}
Info<< endl;
if (surfaceFormat_ != "none")
{
if (Pstream::master())
{
autoPtr<surfaceWriter> writer(surfaceWriter::New(surfaceFormat_));
writer->write
(
outputDir_/time.timeName(),
"collector",
points_,
faces_,
"massTotal",
faceMassTotal,
false
);
writer->write
(
outputDir_/time.timeName(),
"collector",
points_,
faces_,
"massFlowRate",
faceMassFlowRate,
false
);
}
}
if (resetOnWrite_)
{
forAll(faces_, faceI)
{
massFlowRate_[faceI] = 0.0;
}
timeOld_ = timeNew;
totalTime_ = 0.0;
}
forAll(faces_, faceI)
{
mass_[faceI] = 0.0;
}
// writeProperties();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class CloudType>
Foam::ParticleCollector<CloudType>::ParticleCollector
(
const dictionary& dict,
CloudType& owner
)
:
CloudFunctionObject<CloudType>(dict, owner, typeName),
mode_(mtUnknown),
parcelType_(this->coeffDict().lookupOrDefault("parcelType", -1)),
points_(),
faces_(),
faceTris_(),
nSector_(0),
radius_(),
coordSys_(false),
normal_(this->coeffDict().lookup("normal")),
negateParcelsOppositeNormal_
(
readBool(this->coeffDict().lookup("negateParcelsOppositeNormal"))
),
surfaceFormat_(this->coeffDict().lookup("surfaceFormat")),
resetOnWrite_(this->coeffDict().lookup("resetOnWrite")),
totalTime_(0.0),
mass_(),
massTotal_(),
massFlowRate_(),
log_(this->coeffDict().lookup("log")),
outputFilePtr_(),
outputDir_(owner.mesh().time().path()),
timeOld_(owner.mesh().time().value())
{
if (Pstream::parRun())
{
// Put in undecomposed case (Note: gives problems for
// distributed data running)
outputDir_ =
outputDir_/".."/"postProcessing"/cloud::prefix/owner.name();
}
else
{
outputDir_ =
outputDir_/"postProcessing"/cloud::prefix/owner.name();
}
normal_ /= mag(normal_);
word mode(this->coeffDict().lookup("mode"));
if (mode == "polygon")
{
initPolygons();
}
else if (mode == "concentricCircle")
{
initConcentricCircles();
}
else
{
FatalErrorIn
(
"Foam::ParticleCollector<CloudType>::ParticleCollector"
"("
"const dictionary& dict,"
"CloudType& owner"
")"
)
<< "Unknown mode " << mode << ". Available options are "
<< "polygon and concentricCircle" << exit(FatalError);
}
mass_.setSize(faces_.size(), 0.0);
massTotal_.setSize(faces_.size(), 0.0);
massFlowRate_.setSize(faces_.size(), 0.0);
makeLogFile(faces_, points_, area_);
// readProperties(); AND initialise mass... fields
}
template<class CloudType>
Foam::ParticleCollector<CloudType>::ParticleCollector
(
const ParticleCollector<CloudType>& pc
)
:
CloudFunctionObject<CloudType>(pc),
parcelType_(pc.parcelType_),
points_(pc.points_),
faces_(pc.faces_),
faceTris_(pc.faceTris_),
normal_(pc.normal_),
negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
surfaceFormat_(pc.surfaceFormat_),
resetOnWrite_(pc.resetOnWrite_),
totalTime_(pc.totalTime_),
mass_(pc.mass_),
massTotal_(pc.massTotal_),
massFlowRate_(pc.massFlowRate_),
log_(pc.log_),
outputFilePtr_(),
outputDir_(pc.outputDir_),
timeOld_(0.0)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class CloudType>
Foam::ParticleCollector<CloudType>::~ParticleCollector()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CloudType>
void Foam::ParticleCollector<CloudType>::postMove
(
const parcelType& p,
const label cellI,
const scalar dt
// bool& keepParticle
)
{
if ((parcelType_ != -1) && (parcelType_ != p.typeId()))
{
return;
}
label faceI = -1;
switch (mode_)
{
case mtPolygon:
{
faceI = collectParcelPolygon(p.position(), p.U());
break;
}
case mtConcentricCircle:
{
faceI = collectParcelConcentricCircles(p.position(), p.U());
break;
}
default:
{
}
}
if (faceI != -1)
{
scalar m = p.nParticle()*p.mass();
if (negateParcelsOppositeNormal_)
{
vector Uhat = p.U();
Uhat /= mag(Uhat) + ROOTVSMALL;
if ((Uhat & normal_) < 0)
{
m *= -1.0;
}
}
// add mass contribution
mass_[faceI] += m;
if (nSector_ == 1)
{
mass_[faceI + 1] += m;
mass_[faceI + 2] += m;
mass_[faceI + 3] += m;
}
// remove particle
// keepParticle = false;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,252 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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/>.
Class
Foam::ParticleCollector
Description
Function object to collect the parcel mass- and mass flow rate over a
set of polygons. The polygons are defined as lists of points. If a
parcel is 'collected', it is subsequently flagged to be removed from the
domain.
SourceFiles
ParticleCollector.C
\*---------------------------------------------------------------------------*/
#ifndef ParticleCollector_H
#define ParticleCollector_H
#include "CloudFunctionObject.H"
#include "cylindricalCS.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class ParticleCollector Declaration
\*---------------------------------------------------------------------------*/
template<class CloudType>
class ParticleCollector
:
public CloudFunctionObject<CloudType>
{
public:
enum modeType
{
mtPolygon,
mtConcentricCircle,
mtUnknown
};
private:
// Private Data
// Typedefs
//- Convenience typedef for parcel type
typedef typename CloudType::parcelType parcelType;
//- Collector mode type
modeType mode_;
//- Index of parcel types to collect (-1 by default = all particles)
const label parcelType_;
//- List of points
Field<point> points_;
//- List of faces
List<face> faces_;
// Polygon collector
//- Triangulation of faces
List<List<face> > faceTris_;
// Concentric circles collector
//- Number of sectors per circle
label nSector_;
//- List of radii
List<scalar> radius_;
//- Cylindrical co-ordinate system
cylindricalCS coordSys_;
//- Face areas
Field<scalar> area_;
//- Polygon normal vector
vector normal_;
//- Remove mass of parcel travelling in opposite direction to normal_
bool negateParcelsOppositeNormal_;
//- Surface output format
const word surfaceFormat_;
//- Flag to indicate whether data should be reset/cleared on writing
Switch resetOnWrite_;
//- Total time
scalar totalTime_;
//- Mass storage
List<scalar> mass_;
//- Mass total storage
List<scalar> massTotal_;
//- Mass flow rate storage
List<scalar> massFlowRate_;
//- Flag to indicate whether data should be written to file
Switch log_;
//- Output file pointer
autoPtr<OFstream> outputFilePtr_;
//- Output directory
fileName outputDir_;
//- Last calculation time
scalar timeOld_;
// Private Member Functions
//- Helper function to create log files
void makeLogFile
(
const faceList& faces,
const Field<point>& points,
const Field<scalar>& area
);
//- Initialise polygon collectors
void initPolygons();
//- Initialise concentric circle collectors
void initConcentricCircles();
//- Collect parcels in polygon collectors
label collectParcelPolygon
(
const point& position,
const vector& U
) const;
//- Collect parcels in concentric circle collectors
label collectParcelConcentricCircles
(
const point& position,
const vector& U
) const;
protected:
// Protected Member Functions
//- Write post-processing info
void write();
public:
//- Runtime type information
TypeName("particleCollector");
// Constructors
//- Construct from dictionary
ParticleCollector(const dictionary& dict, CloudType& owner);
//- Construct copy
ParticleCollector(const ParticleCollector<CloudType>& pc);
//- Construct and return a clone
virtual autoPtr<CloudFunctionObject<CloudType> > clone() const
{
return autoPtr<CloudFunctionObject<CloudType> >
(
new ParticleCollector<CloudType>(*this)
);
}
//- Destructor
virtual ~ParticleCollector();
// Member Functions
// Access
//- Return const access to the reset on write flag
inline const Switch& resetOnWrite() const;
// Evaluation
//- Post-move hook
virtual void postMove
(
const parcelType& p,
const label cellI,
const scalar dt
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "ParticleCollectorI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "ParticleCollector.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,34 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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/>.
\*---------------------------------------------------------------------------*/
template<class CloudType>
inline const Foam::Switch&
Foam::ParticleCollector<CloudType>::resetOnWrite() const
{
return resetOnWrite_;
}
// ************************************************************************* //