ENH: Added new extractEulerianParticles function object

Generates discrete particle data from multiphase calculations by
interrogating the phase fraction field at a faceZone.

Data is written in raw form, i.e. per particle collected, with
as an optional binned distribution
This commit is contained in:
Andrew Heather
2016-11-29 14:56:48 +00:00
parent 47439e4917
commit 271c8c8c6e
8 changed files with 1654 additions and 5 deletions

View File

@ -72,6 +72,8 @@ externalCoupled/externalCoupled.C
externalCoupled/externalCoupledMixed/externalCoupledMixedFvPatchFields.C externalCoupled/externalCoupledMixed/externalCoupledMixedFvPatchFields.C
externalCoupled/externalCoupledTemperatureMixed/externalCoupledTemperatureMixedFvPatchScalarField.C externalCoupled/externalCoupledTemperatureMixed/externalCoupledTemperatureMixedFvPatchScalarField.C
extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C
ddt2/ddt2.C ddt2/ddt2.C
zeroGradient/zeroGradient.C zeroGradient/zeroGradient.C

View File

@ -2,6 +2,7 @@ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \ -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude \ -I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \ -I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \ -I$(LIB_SRC)/surfMesh/lnInclude \
@ -19,18 +20,22 @@ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \ -I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude -I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/fvAgglomerationMethods/pairPatchAgglomeration/lnInclude
LIB_LIBS = \ LIB_LIBS = \
-lfiniteVolume \ -lfiniteVolume \
-lmeshTools \
-llagrangian \
-ldistributionModels \
-lsampling \
-lsurfMesh \
-lfluidThermophysicalModels \ -lfluidThermophysicalModels \
-lincompressibleTransportModels \ -lincompressibleTransportModels \
-lturbulenceModels \ -lturbulenceModels \
-lcompressibleTransportModels \ -lcompressibleTransportModels \
-lincompressibleTurbulenceModels \ -lincompressibleTurbulenceModels \
-lcompressibleTurbulenceModels \ -lcompressibleTurbulenceModels \
-lmeshTools \
-lsampling \
-lsurfMesh \
-lchemistryModel \ -lchemistryModel \
-lreactionThermophysicalModels -lreactionThermophysicalModels \
-lpairPatchAgglomeration

View File

@ -0,0 +1,104 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ 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 "eulerianParticle.H"
#include "mathematicalConstants.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
Foam::functionObjects::eulerianParticle::eulerianParticle()
:
globalFaceIHit(-1),
VC(vector::zero),
VU(vector::zero),
V(0),
time(0),
timeIndex(0)
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<(Ostream& os, const eulerianParticle& p)
{
os << p.globalFaceIHit << token::SPACE
<< p.VC << token::SPACE
<< p.VU << token::SPACE
<< p.V << token::SPACE
<< p.time << token::SPACE
<< p.timeIndex;
return os;
}
Foam::Istream& Foam::operator>>(Istream& is, eulerianParticle& p)
{
is >> p.globalFaceIHit
>> p.VC
>> p.VU
>> p.V
>> p.time
>> p.timeIndex;
return is;
}
void Foam::functionObjects::eulerianParticle::write(Ostream& os) const
{
scalar pDiameter = cbrt(6*V/constant::mathematical::pi);
vector U = VU/(V + ROOTVSMALL);
vector C = VC/(V + ROOTVSMALL);
os << time << token::SPACE
<< globalFaceIHit << token::SPACE
<< C << token::SPACE
<< pDiameter << token::SPACE
<< U << token::SPACE
<< endl;
}
Foam::dictionary Foam::functionObjects::eulerianParticle::writeDict() const
{
scalar pDiameter = cbrt(6*V/constant::mathematical::pi);
vector U = VU/(V + ROOTVSMALL);
vector C = VC/(V + ROOTVSMALL);
dictionary dict;
dict.add("time", time);
dict.add("meshFace", globalFaceIHit);
dict.add("position", C);
dict.add("diameter", pDiameter);
dict.add("U", U);
return dict;
}
// ************************************************************************* //

View File

@ -0,0 +1,155 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ 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::eulerianParticle
Description
Lightweight class to store particle data derived from VOF calculations,
with special handling for input, output and parallel reduction.
SourceFiles
eulerianParticle.H
eulerianParticle.C
eulerianParticleTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef functionObjects_eulerianParticle_H
#define functionObjects_eulerianParticle_H
#include "label.H"
#include "scalar.H"
#include "vector.H"
#include "dictionary.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class Istream;
class Ostream;
namespace functionObjects
{
class eulerianParticle;
}
// Forward declaration of friend functions and operators
Istream& operator>>(Istream&, functionObjects::eulerianParticle&);
Ostream& operator<<(Ostream&, const functionObjects::eulerianParticle&);
namespace functionObjects
{
/*---------------------------------------------------------------------------*\
Class eulerianParticle Declaration
\*---------------------------------------------------------------------------*/
class eulerianParticle
{
public:
// Public data
//- Index of face in faceZone that this particle hits. Also used to
// identify the index of the coarse face of the surface agglomeration
// Note: value of -1 used to indicate that the particle has not
// been initialised
label globalFaceIHit;
//- Volume multiplied by face centres [m4]
vector VC;
//- Volume multiplied by velocity [m4/s]
vector VU;
//- Volume [m3]
scalar V;
//- Injection time - set at collection [s]
scalar time;
//- Index of last output time
label timeIndex;
//- Constructor
eulerianParticle();
// Public Member Functions
//- Write to stream
void write(Ostream& os) const;
//- Write to dictionary
Foam::dictionary writeDict() const;
// Operators
friend bool operator==
(
const eulerianParticle& a,
const eulerianParticle& b
)
{
return
a.globalFaceIHit == b.globalFaceIHit
&& a.VC == b.VC
&& a.VU == b.VU
&& a.V == b.V
&& a.time == b.time
&& a.timeIndex == b.timeIndex;
}
friend bool operator!=
(
const eulerianParticle& a,
const eulerianParticle& b
)
{
return !(a == b);
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace functionObjects
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "eulerianParticleTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,88 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ 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/>.
\*---------------------------------------------------------------------------*/
namespace Foam
{
namespace functionObjects
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class VOFParticle>
class sumParticleOp
{
public:
eulerianParticle operator()
(
const eulerianParticle& p0,
const eulerianParticle& p1
) const
{
if ((p0.globalFaceIHit != -1) && (p1.globalFaceIHit == -1))
{
return p0;
}
else if ((p0.globalFaceIHit == -1) && (p1.globalFaceIHit != -1))
{
return p1;
}
else if ((p0.globalFaceIHit != -1) && (p1.globalFaceIHit != -1))
{
// Choose particle with the largest collected volume and
// accumulate total volume
if (p0.V > p1.V)
{
eulerianParticle p = p0;
p.V = p0.V + p1.V;
p.VC = p0.VC + p1.VC;
p.VU = p0.VU + p1.VU;
return p;
}
else
{
eulerianParticle p = p1;
p.V = p0.V + p1.V;
p.VC = p0.VC + p1.VC;
p.VU = p0.VU + p1.VU;
return p;
}
}
else
{
eulerianParticle p;
return p;
}
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace functionObjects
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -0,0 +1,893 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ 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 "extractEulerianParticles.H"
#include "regionSplit2D.H"
#include "mathematicalConstants.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "surfaceInterpolate.H"
#include "pairPatchAgglomeration.H"
#include "emptyPolyPatch.H"
#include "coupledPolyPatch.H"
#include "binned.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace functionObjects
{
defineTypeNameAndDebug(extractEulerianParticles, 0);
addToRunTimeSelectionTable
(
functionObject,
extractEulerianParticles,
dictionary
);
}
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::fileName
Foam::functionObjects::extractEulerianParticles::dictBaseFileDir() const
{
fileName baseDir(".."); // = obr_.time().path();
if (Pstream::parRun())
{
// Put in undecomposed case (Note: gives problems for
// distributed data running)
baseDir = baseDir/".."/"postProcessing";
}
else
{
baseDir = baseDir/"postProcessing";
}
return baseDir;
}
void Foam::functionObjects::extractEulerianParticles::checkFaceZone()
{
DebugInFunction << endl;
zoneID_ = mesh_.faceZones().findZoneID(faceZoneName_);
if (zoneID_ == -1)
{
FatalErrorInFunction
<< "Unable to find faceZone " << faceZoneName_
<< ". Available faceZones are: "
<< exit(FatalError);
}
const faceZone& fz = mesh_.faceZones()[zoneID_];
const label nFaces = fz.size();
const label allFaces = returnReduce(nFaces, sumOp<label>());
if (allFaces < nInjectorLocations_)
{
FatalErrorInFunction
<< "faceZone " << faceZoneName_
<< ": Number of faceZone faces (" << allFaces
<< ") is less than the number of requested locations ("
<< nInjectorLocations_ << ")."
<< exit(FatalError);
}
Log << type() << " " << name() << " output:" << nl
<< " faceZone : " << faceZoneName_ << nl
<< " faces : " << allFaces << nl
<< endl;
// Initialise old iteration blocked faces
// Note: for restart, this info needs to be written/read
regions0_.setSize(fz.size(), -1);
// Create global addressing for faceZone
globalFaces_ = globalIndex(fz.size());
}
void Foam::functionObjects::extractEulerianParticles::initialiseBins()
{
DebugInFunction << endl;
const faceZone& fz = mesh_.faceZones()[zoneID_];
// Agglomerate faceZone faces into nInjectorLocations_ global locations
const indirectPrimitivePatch patch
(
IndirectList<face>(mesh_.faces(), fz),
mesh_.points()
);
const label nFaces = fz.size();
label nLocations = nInjectorLocations_;
if (Pstream::parRun())
{
label nGlobalFaces = returnReduce(nFaces, sumOp<label>());
scalar fraction = scalar(nFaces)/scalar(nGlobalFaces);
nLocations = ceil(fraction*nInjectorLocations_);
if (debug)
{
Pout<< "nFaces:" << nFaces
<< ", nGlobalFaces:" << nGlobalFaces
<< ", fraction:" << fraction
<< ", nLocations:" << nLocations
<< endl;
}
}
pairPatchAgglomeration ppa(patch, 10, 50, nLocations, labelMax, 180);
ppa.agglomerate();
label nCoarseFaces = 0;
if (nFaces != 0)
{
fineToCoarseAddr_ = ppa.restrictTopBottomAddressing();
nCoarseFaces = max(fineToCoarseAddr_) + 1;
coarseToFineAddr_ = invertOneToMany(nCoarseFaces, fineToCoarseAddr_);
// Set coarse face centres as area average of fine face centres
const vectorField& faceCentres = mesh_.faceCentres();
const vectorField& faceAreas = mesh_.faceAreas();
coarsePosition_.setSize(coarseToFineAddr_.size());
forAll(coarseToFineAddr_, coarsei)
{
const labelList& fineFaces = coarseToFineAddr_[coarsei];
scalar sumArea = 0;
vector averagePosition(vector::zero);
forAll(fineFaces, i)
{
label facei = fz[fineFaces[i]];
scalar magSf = mag(faceAreas[facei]);
sumArea += magSf;
averagePosition += magSf*faceCentres[facei];
}
coarsePosition_[coarsei] = averagePosition/sumArea;
}
}
// Create global addressing for coarse face addressing
globalCoarseFaces_ = globalIndex(nCoarseFaces);
Info<< "Created " << returnReduce(nCoarseFaces, sumOp<label>())
<< " coarse faces" << endl;
}
Foam::tmp<Foam::surfaceScalarField>
Foam::functionObjects::extractEulerianParticles::phiU() const
{
DebugInFunction << endl;
const surfaceScalarField& phi
(
mesh_.lookupObject<surfaceScalarField>(phiName_)
);
if (phi.dimensions() == dimMass/dimTime)
{
const volScalarField& rho =
obr_.lookupObject<volScalarField>(rhoName_);
return phi/fvc::interpolate(rho);
}
return phi;
}
void Foam::functionObjects::extractEulerianParticles::setBlockedFaces
(
const surfaceScalarField& alphaf,
const faceZone& fz,
boolList& blockedFaces
)
{
DebugInFunction << endl;
// Initialise storage for patch and patch-face indices where faceZone
// intersects mesh patch(es)
patchIDs_.setSize(fz.size(), -1);
patchFaceIDs_.setSize(fz.size(), -1);
label nBlockedFaces = 0;
forAll(fz, localFacei)
{
const label facei = fz[localFacei];
if (mesh_.isInternalFace(facei))
{
if (alphaf[facei] > alphaThreshold_)
{
blockedFaces[localFacei] = true;
}
}
else
{
label patchi = mesh_.boundaryMesh().whichPatch(facei);
label patchFacei = -1;
const polyPatch& pp = mesh_.boundaryMesh()[patchi];
const scalarField& alphafp = alphaf.boundaryField()[patchi];
if (isA<coupledPolyPatch>(pp))
{
const coupledPolyPatch& cpp =
refCast<const coupledPolyPatch>(pp);
if (cpp.owner())
{
patchFacei = cpp.whichFace(facei);
}
}
else if (!isA<emptyPolyPatch>(pp))
{
patchFacei = pp.whichFace(facei);
}
if (patchFacei == -1)
{
patchi = -1;
}
else if (alphafp[patchFacei] > alphaThreshold_)
{
blockedFaces[localFacei] = true;
}
patchIDs_[localFacei] = patchi;
patchFaceIDs_[localFacei] = patchFacei;
}
}
DebugInFunction << "Number of blocked faces: " << nBlockedFaces << endl;
}
void Foam::functionObjects::extractEulerianParticles::collectParticles
(
const scalar time,
const boolList& collectParticle
)
{
DebugInFunction << "collectParticle: " << collectParticle << endl;
// Collect particles on local processors that have passed through faceZone
forAll(collectParticle, regioni)
{
if (!collectParticle[regioni])
{
continue;
}
Map<label>::const_iterator iter = regionToParticleMap_.find(regioni);
eulerianParticle p;
if (iter != Map<label>::end())
{
// Particle on local processor
p = particles_[iter()];
}
reduce(p, sumParticleOp<eulerianParticle>());
const scalar pDiameter = cbrt(6.0*p.V/constant::mathematical::pi);
if ((pDiameter > minDiameter_) && (pDiameter < maxDiameter_))
{
p.time = time;
p.timeIndex = outputTimes_.size() - 1;
if (Pstream::master() && writeRawData_)
{
rawParticlesDictPtr_->add
(
word("particle" + Foam::name(nCollectedParticles_)),
p.writeDict()
);
}
if (globalFaces_.isLocal(p.globalFaceIHit))
{
collectedParticles_.append(p);
}
nCollectedParticles_++;
}
else
{
// Discard particles over/under diameter threshold
nDiscardedParticles_++;
discardedVolume_ += p.V;
}
}
}
void Foam::functionObjects::extractEulerianParticles::calculateAddressing
(
const label nRegionsOld,
const label nRegionsNew,
const scalar time,
labelList& regionFaceIDs
)
{
DebugInFunction << endl;
// New region can only point to one old region
// Old region can only point to one new region. If old region intersects
// multiple new regions, select max of new region indices.
labelList oldToNewRegion(nRegionsOld, -1);
labelList newToOldRegion(nRegionsNew, -1);
forAll(regionFaceIDs, facei)
{
label newRegioni = regionFaceIDs[facei];
label oldRegioni = regions0_[facei];
if (newRegioni != -1)
{
newToOldRegion[newRegioni] = oldRegioni;
if (oldRegioni != -1)
{
// New region linked to old (so can accumulate particle data)
// Old region might already be mapped to a new region
oldToNewRegion[oldRegioni] =
max(newRegioni, oldToNewRegion[oldRegioni]);
}
}
}
// Need to re-number the new region indices based on knowledge of which
// old region they intersect. After operations, there should be a
// one-to-one match between the old and new regions.
// Ensure all old regions point to the same new regions (if present)
Pstream::listCombineGather(oldToNewRegion, maxEqOp<label>());
Pstream::listCombineScatter(oldToNewRegion);
// Any new region that points to an old region should be renumbered to the
// new region specified by the oldToNewRegion index
if (oldToNewRegion.size())
{
// Create corrected new to new addressing
labelList newToNewRegionCorr(newToOldRegion.size(), -1);
forAll(newToOldRegion, newRegioni)
{
label oldRegioni = newToOldRegion[newRegioni];
if (oldRegioni != -1)
{
label newRegionICorr = oldToNewRegion[oldRegioni];
newToNewRegionCorr[newRegioni] = newRegionICorr;
newToOldRegion[newRegioni] = oldRegioni;
}
}
// Renumber the new (current) face region IDs
forAll(regionFaceIDs, facei)
{
label newRegioni = regionFaceIDs[facei];
if (newRegioni != -1)
{
label newRegionICorr = newToNewRegionCorr[newRegioni];
// Note: newRegionICorr can be -1 if the region is new, since
// the address corrections are based on inverting the
// old-to-new addressing
if (newRegionICorr != -1)
{
regionFaceIDs[facei] = newRegionICorr;
}
}
}
boolList collectParticleFlag(nRegionsOld, true);
forAll(oldToNewRegion, oldRegioni)
{
label newRegioni = oldToNewRegion[oldRegioni];
if (newRegioni != -1)
{
collectParticleFlag[oldRegioni] = false;
}
}
// Collect particles whose IDs are no longer active
collectParticles(time, collectParticleFlag);
}
// Re-order collection bins
Map<label> newRegionToParticleMap(nRegionsNew);
List<eulerianParticle> newParticles(nRegionsNew);
label particlei = 0;
forAll(newToOldRegion, newRegioni)
{
label oldRegioni = newToOldRegion[newRegioni];
if (oldRegioni == -1)
{
// No mapping from old to new - this is a new particle
newRegionToParticleMap.insert(newRegioni, particlei);
particlei++;
}
else
{
// Update addressing for old to new regions
label oldParticlei = regionToParticleMap_[oldRegioni];
if (newRegionToParticleMap.insert(newRegioni, particlei))
{
// First time this particle has been seen
newParticles[particlei] = particles_[oldParticlei];
particlei++;
}
else
{
// Combine with existing contribution
label newParticlei = newRegionToParticleMap[newRegioni];
newParticles[newParticlei] =
sumParticleOp<eulerianParticle>()
(
newParticles[newParticlei],
particles_[oldParticlei]
);
}
}
}
// Reset the particles list and addressing for latest available info
particles_.transfer(newParticles);
regionToParticleMap_ = newRegionToParticleMap;
}
void Foam::functionObjects::extractEulerianParticles::accumulateParticleInfo
(
const surfaceScalarField& alphaf,
const surfaceScalarField& phi,
const labelList& regionFaceIDs,
const faceZone& fz
)
{
DebugInFunction << endl;
const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
const surfaceVectorField Uf(fvc::interpolate(U));
const scalar deltaT = mesh_.time().deltaTValue();
const pointField& faceCentres = mesh_.faceCentres();
forAll(regionFaceIDs, localFacei)
{
const label newRegioni = regionFaceIDs[localFacei];
if (newRegioni != -1)
{
const label particlei = regionToParticleMap_[newRegioni];
const label meshFacei = fz[localFacei];
eulerianParticle& p = particles_[particlei];
if (p.globalFaceIHit < 0)
{
// New particle - does not exist in particles_ list
p.globalFaceIHit = globalFaces_.toGlobal(localFacei);
p.V = 0;
p.VC = vector::zero;
p.VU = vector::zero;
}
// Accumulate particle properties
scalar magPhii = mag(faceValue(phi, localFacei, meshFacei));
vector Ufi = faceValue(Uf, localFacei, meshFacei);
scalar dV = magPhii*deltaT;
p.V += dV;
p.VC += dV*faceCentres[meshFacei];
p.VU += dV*Ufi;
}
}
}
void Foam::functionObjects::extractEulerianParticles::writeBinnedParticleData()
{
DebugInFunction << endl;
if (!createDistribution_)
{
// No need to store collected particles if not creating a distribution
collectedParticles_.clear();
return;
}
// Gather particles ready for collection from all procs
List<List<eulerianParticle> > allProcParticles(Pstream::nProcs());
allProcParticles[Pstream::myProcNo()] = collectedParticles_;
Pstream::gatherList(allProcParticles);
Pstream::scatterList(allProcParticles);
List<eulerianParticle> allParticles =
ListListOps::combine<List<eulerianParticle> >
(
allProcParticles,
accessOp<List<eulerianParticle> >()
);
// Determine coarse face index (global) and position for each particle
label nCoarseFaces = globalCoarseFaces_.size();
List<label> particleCoarseFacei(allParticles.size(), -1);
List<point> particleCoarseFacePosition(nCoarseFaces, point::min);
forAll(allParticles, particlei)
{
const eulerianParticle& p = allParticles[particlei];
label globalFaceHiti = p.globalFaceIHit;
if (globalFaces_.isLocal(globalFaceHiti))
{
label localFacei = globalFaces_.toLocal(globalFaceHiti);
label coarseFacei = fineToCoarseAddr_[localFacei];
label globalCoarseFacei = globalCoarseFaces_.toGlobal(coarseFacei);
particleCoarseFacei[particlei] = globalCoarseFacei;
particleCoarseFacePosition[globalCoarseFacei] =
coarsePosition_[coarseFacei];
}
}
Pstream::listCombineGather(particleCoarseFacei, maxEqOp<label>());
Pstream::listCombineGather(particleCoarseFacePosition, maxEqOp<point>());
// Write the agglomerated particle data to file
DynamicList<label> processedCoarseFaces;
if (Pstream::master())
{
fileName baseDir(dictBaseFileDir()/name());
IOdictionary dict
(
IOobject
(
"particleDistribution",
obr_.time().timeName(),
baseDir,
obr_,
IOobject::NO_READ,
IOobject::NO_WRITE
)
);
labelListList coarseFaceToParticle =
invertOneToMany(nCoarseFaces, particleCoarseFacei);
// Process the allParticles per coarse face
forAll(coarseFaceToParticle, globalCoarseFacei)
{
const List<label>& particleIDs =
coarseFaceToParticle[globalCoarseFacei];
const label nParticle = particleIDs.size();
if (nParticle == 0)
{
continue;
}
Field<scalar> pd(particleIDs.size());
scalar sumV = 0;
vector sumVU = vector::zero;
scalar startTime = GREAT;
scalar endTime = -GREAT;
forAll(particleIDs, i)
{
const label particlei = particleIDs[i];
const eulerianParticle& p = allParticles[particlei];
scalar pDiameter = cbrt(6*p.V/constant::mathematical::pi);
pd[i] = pDiameter;
sumV += p.V;
sumVU += p.VU;
startTime = min(startTime, outputTimes_[p.timeIndex]);
endTime = max(endTime, outputTimes_[p.timeIndex + 1]);
}
if (sumV < ROOTVSMALL)
{
// Started collecting particle info, but not accumulated any
// volume yet
continue;
}
distributionModels::binned binnedDiameters
(
pd,
distributionBinWidth_,
rndGen_
);
// Velocity info hard-coded to volume average
vector Uave = sumVU/sumV;
dictionary particleDict;
particleDict.add("startTime", startTime);
particleDict.add("endTime", endTime);
particleDict.add("nParticle", nParticle);
particleDict.add
(
"position",
particleCoarseFacePosition[globalCoarseFacei]
);
particleDict.add("volume", sumV);
particleDict.add("U", Uave);
particleDict.add
(
"binnedDistribution",
binnedDiameters.writeDict("distribution")
);
dict.add
(
word("sample" + Foam::name(globalCoarseFacei)),
particleDict
);
processedCoarseFaces.append(globalCoarseFacei);
}
dict.regIOobject::write();
}
if (resetDistributionOnWrite_)
{
// Remove processed coarse faces from collectedParticles_
Pstream::scatter(processedCoarseFaces);
labelHashSet processedFaces(processedCoarseFaces);
DynamicList<eulerianParticle> nonProcessedParticles;
forAll(collectedParticles_, particlei)
{
const eulerianParticle& p = collectedParticles_[particlei];
label localFacei = globalFaces_.toLocal(p.globalFaceIHit);
label coarseFacei = fineToCoarseAddr_[localFacei];
label globalCoarseFacei = globalCoarseFaces_.toGlobal(coarseFacei);
if (!processedFaces.found(globalCoarseFacei))
{
nonProcessedParticles.append(p);
}
}
collectedParticles_.transfer(nonProcessedParticles);
}
}
// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
Foam::functionObjects::extractEulerianParticles::extractEulerianParticles
(
const word& name,
const Time& runTime,
const dictionary& dict
)
:
fvMeshFunctionObject(name, runTime, dict),
writeFile(runTime, name),
writeRawData_(false),
rawParticlesDictPtr_(NULL),
outputTimes_(100),
faceZoneName_(word::null),
zoneID_(-1),
patchIDs_(),
patchFaceIDs_(),
alphaName_("alpha"),
alphaThreshold_(0.1),
UName_("U"),
rhoName_("rho"),
phiName_("phi"),
nInjectorLocations_(0),
fineToCoarseAddr_(),
coarseToFineAddr_(),
coarsePosition_(),
globalCoarseFaces_(),
regions0_(),
nRegions0_(0),
particles_(),
regionToParticleMap_(),
collectedParticles_(),
minDiameter_(ROOTVSMALL),
maxDiameter_(GREAT),
createDistribution_(false),
resetDistributionOnWrite_(false),
distributionBinWidth_(0),
rndGen_(1234, -1),
nCollectedParticles_(0),
nDiscardedParticles_(0),
discardedVolume_(0)
{
if (mesh_.nSolutionD() != 3)
{
FatalErrorInFunction
<< name << " function object only applicable to 3-D cases"
<< exit(FatalError);
}
read(dict);
if (Pstream::master() && writeRawData_)
{
fileName baseDir(dictBaseFileDir()/name);
rawParticlesDictPtr_.reset
(
new IOdictionary
(
IOobject
(
"rawParticlesDict",
mesh_.time().timeName(),
baseDir,
mesh_,
IOobject::NO_READ, // _IF_PRESENT,
IOobject::NO_WRITE
)
)
);
DebugVar(rawParticlesDictPtr_->objectPath());
DebugVar(rawParticlesDictPtr_->path());
}
outputTimes_.append(mesh_.time().value());
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::functionObjects::extractEulerianParticles::~extractEulerianParticles()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::functionObjects::extractEulerianParticles::read(const dictionary& dict)
{
DebugInFunction << endl;
if (fvMeshFunctionObject::read(dict) && writeFile::read(dict))
{
dict.lookup("writeRawData") >> writeRawData_;
dict.lookup("faceZone") >> faceZoneName_;
dict.lookup("nLocations") >> nInjectorLocations_;
dict.lookup("alphaName") >> alphaName_;
dict.readIfPresent("alphaThreshold", alphaThreshold_);
dict.lookup("UName") >> UName_;
dict.lookup("rhoName") >> rhoName_;
dict.lookup("phiName") >> phiName_;
dict.readIfPresent("minDiameter", minDiameter_);
dict.readIfPresent("maxDiameter", maxDiameter_);
dict.lookup("createDistribution") >> createDistribution_;
dict.lookup("distributionBinWidth") >> distributionBinWidth_;
dict.lookup("resetDistributionOnWrite") >> resetDistributionOnWrite_;
checkFaceZone();
initialiseBins();
return true;
}
return false;
}
bool Foam::functionObjects::extractEulerianParticles::execute()
{
DebugInFunction << endl;
Log << type() << " " << name() << " output:" << nl;
const volScalarField& alpha =
obr_.lookupObject<volScalarField>(alphaName_);
const surfaceScalarField alphaf
(
typeName + ":alphaf",
fvc::interpolate(alpha)
);
const faceZone& fz = mesh_.faceZones()[zoneID_];
const indirectPrimitivePatch patch
(
IndirectList<face>(mesh_.faces(), fz),
mesh_.points()
);
// Set the blocked faces, i.e. where alpha > alpha threshold value
boolList blockedFaces(fz.size(), false);
setBlockedFaces(alphaf, fz, blockedFaces);
// Split the faceZone according to the blockedFaces
// - Returns a list of (disconnected) region index per face zone face
regionSplit2D regionFaceIDs(mesh_, patch, blockedFaces);
const label nRegionsNew = regionFaceIDs.nRegions();
// Calculate the addressing between the old and new region information
// Also collects particles that have traversed the faceZone
calculateAddressing
(
nRegions0_,
nRegionsNew,
mesh_.time().value(),
regionFaceIDs
);
// Process latest region information
tmp<surfaceScalarField> tphi = phiU();
accumulateParticleInfo(alphaf, tphi(), regionFaceIDs, fz);
// Reset the blocked faces fot the next integration step
nRegions0_ = nRegionsNew;
regions0_ = regionFaceIDs;
Log << " Collected particles: " << nCollectedParticles_ << nl
<< " Discarded particles: " << nDiscardedParticles_ << nl
<< " Discarded volume: " << discardedVolume_ << nl
<< endl;
return true;
}
bool Foam::functionObjects::extractEulerianParticles::write()
{
DebugInFunction << endl;
outputTimes_.append(obr_.time().value());
writeBinnedParticleData();
if (Pstream::master() && writeRawData_)
{
rawParticlesDictPtr_->regIOobject::write();
}
return true;
}
// ************************************************************************* //

View File

@ -0,0 +1,346 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ 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::functionObjects::extractEulerianParticles
Group
grpFieldFunctionObjects
Description
Generates particle size information from Eulerian calculations, e.g. VoF.
Particle data is written to the directory:
\verbatim
$FOAM_CASE/postProcessing/<name>/rawParticlesDict
$FOAM_CASE/postProcessing/<name>/particleDistribution
\endverbatim
Usage
extractEulerianParticles1
{
type extractEulerianParticles;
libs ("libfieldFunctionObjects.so");
...
faceZone f0;
nLocations 10;
alphaName alpha.water;
UName U;
rhoName rho;
phiName phi;
createDistribution yes;
distributionBinWidth 1e-4;
resetDistributionOnWrite no;
writeRawData yes;
}
\endverbatim
where the entries comprise:
\table
Property | Description | Required | Default value
type | type name: extractEulerianParticles | yes |
faceZone | Name of faceZone used as collection surface | yes |
nLocations | Number of injection bins to generate | yes |
aplhaName | Name of phase indicator field | yes |
rhoName | Name of density field | yes |
phiNane | Name of flux field | yes |
createDistribution | Flag to create a binned distribution | yes |
distributionBinWidth | Binned distribution bin width| yes |
writeRawData | Flag to write raw particle data | yes |
\endtable
SourceFiles
extractEulerianParticles.C
\*---------------------------------------------------------------------------*/
#ifndef functionObjects_extractEulerianParticles_H
#define functionObjects_extractEulerianParticles_H
#include "fvMeshFunctionObject.H"
#include "writeFile.H"
#include "runTimeSelectionTables.H"
#include "polyMesh.H"
#include "surfaceFieldsFwd.H"
#include "vectorList.H"
#include "globalIndex.H"
#include "cachedRandom.H"
#include "eulerianParticle.H"
#include "IOdictionary.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace functionObjects
{
/*---------------------------------------------------------------------------*\
Class extractEulerianParticlesFunctionObject Declaration
\*---------------------------------------------------------------------------*/
class extractEulerianParticles
:
public fvMeshFunctionObject,
public writeFile
{
protected:
// Protected data
//- Flag to write the raw collected particle data
bool writeRawData_;
//- Raw particle output dictionary
autoPtr<IOdictionary> rawParticlesDictPtr_;
//- List to keep track of output times
DynamicList<scalar> outputTimes_;
// faceZone info
//- Name of faceZone to sample
word faceZoneName_;
//- Index of the faceZone
label zoneID_;
//- Patch indices where faceZone face intersect patch
labelList patchIDs_;
//- Patch face indices where faceZone face intersect patch
labelList patchFaceIDs_;
//- Global face addressing
globalIndex globalFaces_;
// Field names
//- Name of phase fraction field
word alphaName_;
//- Value of phase fraction used to identify particle boundaries
scalar alphaThreshold_;
//- Name of the velocity field, default = U
word UName_;
//- Name of the density field, defauls = rho
word rhoName_;
//- Name of the flux field, default ="rho"
word phiName_;
// Agglomeration
//- Number of sample locations to generate
label nInjectorLocations_;
//- Agglomeration addressing from fine to coarse
labelList fineToCoarseAddr_;
//- Agglomeration addressing from coarse to fine
labelListList coarseToFineAddr_;
//- Coarse face positions
vectorList coarsePosition_;
//- Global coarse face addressing
globalIndex globalCoarseFaces_;
// Particle collection info
//- Time of last write step
scalar lastOutputTime_;
//- Region indices in faceZone faces from last iteration
labelList regions0_;
//- Number of regions from last iteration
label nRegions0_;
//- Particle properties (partial, being accumulated)
List<eulerianParticle> particles_;
//- Map from region to index in particles_ list
Map<label> regionToParticleMap_;
// Collected particles
DynamicList<eulerianParticle> collectedParticles_;
//- Minimum diameter (optional)
// Can be used to filter out 'small' particles
scalar minDiameter_;
//- Maximum diameter (optional)
// Can be used to filter out 'large' particles
scalar maxDiameter_;
//- Flag to create a distribution
bool createDistribution_;
//- Flag to reset the distribution on each write
bool resetDistributionOnWrite_;
//- Diameter distribution bin width
scalar distributionBinWidth_;
//- Random class needed by distribution models
cachedRandom rndGen_;
// Statistics
//- Total number of collected particles
label nCollectedParticles_;
//- Total number of discarded particles
label nDiscardedParticles_;
//- Total discarded volume [m3]
scalar discardedVolume_;
// Protected Member Functions
//- Return the base directory for dictionary output
virtual fileName dictBaseFileDir() const;
//- Check that the faceZone is valid
virtual void checkFaceZone();
//- Initialise the particle collection bins
virtual void initialiseBins();
//- Return the volumetric flux
virtual tmp<surfaceScalarField> phiU() const;
//- Set the blocked faces, i.e. where alpha > alpha threshold value
virtual void setBlockedFaces
(
const surfaceScalarField& alphaf,
const faceZone& fz,
boolList& blockedFaces
);
//- Calculate the addressing between regions between iterations
virtual void calculateAddressing
(
const label nRegionsOld,
const label nRegionsNew,
const scalar time,
labelList& regionFaceIDs
);
//- Collect particles that have passed through the faceZone
virtual void collectParticles
(
const scalar time,
const boolList& collectParticle
);
//- Process latest region information
virtual void accumulateParticleInfo
(
const surfaceScalarField& alphaf,
const surfaceScalarField& phi,
const labelList& regionFaceIDs,
const faceZone& fz
);
//- Write agglomerated particle data to stream
virtual void writeBinnedParticleData();
template<class Type>
inline Type faceValue
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& field,
const label localFaceI,
const label globalFaceI
) const;
//- Disallow default bitwise copy construct
extractEulerianParticles(const extractEulerianParticles&);
//- Disallow default bitwise assignment
void operator=(const extractEulerianParticles&);
public:
// Static data members
//- Static data staticData
TypeName("extractEulerianParticles");
// Constructors
//- Construct from components
extractEulerianParticles
(
const word& name,
const Time& runTime,
const dictionary& dict
);
//- Destructor
virtual ~extractEulerianParticles();
// Member Functions
//- Read the field min/max data
virtual bool read(const dictionary&);
//- Execute
virtual bool execute();
//- Write
virtual bool write();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace functionObjects
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "extractEulerianParticlesTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,56 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenCFD Ltd.
\\/ 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 "fvMesh.H"
template<class Type>
Type Foam::functionObjects::extractEulerianParticles::faceValue
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& field,
const label localFacei,
const label globalFacei
) const
{
if (mesh_.isInternalFace(globalFacei))
{
return field[globalFacei];
}
else
{
label patchi = patchIDs_[localFacei];
label pFacei = patchFaceIDs_[localFacei];
if (patchi != 0)
{
return field.boundaryField()[patchi][pFacei];
}
else
{
return pTraits<Type>::zero;
}
}
}
// ************************************************************************* //