zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: Updating particle member functions

Browse Source 
- Dealing with detA < 0 tracking issues
- Modified locate function
This commit is contained in:
Sergio Ferraris
2020-08-28 08:19:06 -07:00
committed by Andrew Heather
parent 8427eccd00
commit f7dc0d8edb
4 changed files with 278 additions and 103 deletions
Download Patch File Download Diff File Expand all files Collapse all files

253
src/lagrangian/basic/particle/particle.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2018 OpenCFD Ltd.
Copyright (C) 2018-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -40,7 +40,7 @@ namespace Foam
defineTypeNameAndDebug(particle, 0);
}
const Foam::scalar Foam::particle::negativeSpaceDisplacementFactor = 1.01;
const Foam::label Foam::particle::maxNBehind_ = 10;
Foam::label Foam::particle::particleCount_ = 0;
@ -365,7 +365,6 @@ void Foam::particle::changeCell()
const label ownerCellI = mesh_.faceOwner()[tetFacei_];
const bool isOwner = celli_ == ownerCellI;
celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
// Reflect to account for the change of triangle orientation in the new cell
reflect();
}
@ -414,78 +413,84 @@ void Foam::particle::locate
const string& boundaryMsg
)
{
if (debug)
{
Info << "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
}
celli_ = celli;
// Find the cell, if it has not been given
if (celli_ < 0)
if (celli < 0)
{
celli_ = mesh_.cellTree().findInside(position);
}
if (celli_ < 0)
if (celli < 0)
{
FatalErrorInFunction
<< "Cell not found for particle position " << position << "."
<< exit(FatalError);
}
// Put the particle at (almost) the cell centre and in a random tet.
// Note perturbing the cell centre to make sure we find at least one
// tet containing it. With start point exactly at the cell centre very
// occasionally it would not get found in any of the tets
coordinates_ = barycentric(1-3*SMALL, SMALL, SMALL, SMALL);
tetFacei_ = mesh_.cells()[celli_][0];
tetPti_ = 1;
const vector displacement = position - mesh_.cellCentres()[celli_];
// Loop all cell tets to find the one containing the position. Track
// through each tet from the cell centre. If a tet contains the position
// then the track will end with a single trackToTri.
const class cell& c = mesh_.cells()[celli_];
scalar minF = VGREAT;
label minTetFacei = -1, minTetPti = -1;
forAll(c, cellTetFacei)
{
const class face& f = mesh_.faces()[c[cellTetFacei]];
for (label tetPti = 1; tetPti < f.size() - 1; ++tetPti)
{
coordinates_ = barycentric(1, 0, 0, 0);
tetFacei_ = c[cellTetFacei];
tetPti_ = tetPti;
facei_ = -1;
// Track to the injection point
track(position - this->position(), 0);
label tetTriI = -1;
const scalar f = trackToTri(displacement, 0, tetTriI);
if (tetTriI == -1)
{
return;
}
if (f < minF)
{
minF = f;
minTetFacei = tetFacei_;
minTetPti = tetPti_;
}
}
}
// The particle must be (hopefully only slightly) outside the cell. Track
// into the tet which got the furthest.
coordinates_ = barycentric(1, 0, 0, 0);
tetFacei_ = minTetFacei;
tetPti_ = minTetPti;
facei_ = -1;
track(displacement, 0);
if (!onFace())
{
return;
}
// We hit a boundary ...
// If we are here then we hit a boundary
if (boundaryFail)
{
FatalErrorInFunction << boundaryMsg
<< " when tracking from centre " << mesh_.cellCentres()[celli_]
<< " of cell " << celli_ << " to position " << position
<< exit(FatalError);
FatalErrorInFunction << boundaryMsg << exit(FatalError);
}
else
{
// Re-do the track, but this time do the bit tangential to the
// direction/patch first. This gets us as close as possible to the
// original path/position.
if (direction == nullptr)
{
const polyPatch& p = mesh_.boundaryMesh()[patch()];
direction = &p.faceNormals()[p.whichFace(facei_)];
}
const vector n = *direction/mag(*direction);
const vector s = position - mesh_.cellCentres()[celli_];
const vector sN = (s & n)*n;
const vector sT = s - sN;
coordinates_ = barycentric(1, 0, 0, 0);
tetFacei_ = mesh_.cells()[celli_][0];
tetPti_ = 1;
facei_ = -1;
track(sT, 0);
track(sN, 0);
static label nWarnings = 0;
static const label maxNWarnings = 100;
if (nWarnings < maxNWarnings)
if ((nWarnings < maxNWarnings) && boundaryFail)
{
WarningInFunction << boundaryMsg.c_str()
<< " when tracking from centre " << mesh_.cellCentres()[celli_]
<< " of cell " << celli_ << " to position " << position
<< endl;
WarningInFunction << boundaryMsg.c_str() << endl;
++ nWarnings;
}
if (nWarnings == maxNWarnings)
@ -496,6 +501,7 @@ void Foam::particle::locate
++ nWarnings;
}
}
}
@ -516,7 +522,9 @@ Foam::particle::particle
tetFacei_(tetFacei),
tetPti_(tetPti),
facei_(-1),
stepFraction_(0.0),
stepFraction_(1.0),
behind_(0.0),
nBehind_(0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{}
@ -535,7 +543,9 @@ Foam::particle::particle
tetFacei_(-1),
tetPti_(-1),
facei_(-1),
stepFraction_(0.0),
stepFraction_(1.0),
behind_(0.0),
nBehind_(0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{
@ -566,7 +576,9 @@ Foam::particle::particle
tetFacei_(tetFacei),
tetPti_(tetPti),
facei_(-1),
stepFraction_(0.0),
stepFraction_(1.0),
behind_(0.0),
nBehind_(0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{
@ -593,6 +605,8 @@ Foam::particle::particle(const particle& p)
tetPti_(p.tetPti_),
facei_(p.facei_),
stepFraction_(p.stepFraction_),
behind_(p.behind_),
nBehind_(p.nBehind_),
origProc_(p.origProc_),
origId_(p.origId_)
{}
@ -607,6 +621,8 @@ Foam::particle::particle(const particle& p, const polyMesh& mesh)
tetPti_(p.tetPti_),
facei_(p.facei_),
stepFraction_(p.stepFraction_),
behind_(p.behind_),
nBehind_(p.nBehind_),
origProc_(p.origProc_),
origId_(p.origId_)
{}
@ -645,7 +661,8 @@ Foam::scalar Foam::particle::trackToFace
facei_ = -1;
while (true)
// Loop the tets in the current cell
while (nBehind_ < maxNBehind_)
{
f *= trackToTri(f*displacement, f*fraction, tetTriI);
@ -666,6 +683,25 @@ Foam::scalar Foam::particle::trackToFace
changeTet(tetTriI);
}
}
// Warn if stuck, and incorrectly advance the step fraction to completion
static label stuckID = -1, stuckProc = -1;
if (origId_ != stuckID && origProc_ != stuckProc)
{
WarningInFunction
<< "Particle #" << origId_ << " got stuck at " << position()
<< endl;
}
stuckID = origId_;
stuckProc = origProc_;
stepFraction_ += f*fraction;
behind_ = 0;
nBehind_ = 0;
return 0;
}
@ -702,11 +738,8 @@ Foam::scalar Foam::particle::trackToStationaryTri
<< "Start local coordinates = " << y0 << endl;
}
// Get the factor by which the displacement is increased
const scalar f = detA >= 0 ? 1 : negativeSpaceDisplacementFactor;
// Calculate the local tracking displacement
barycentric Tx1(f*x1 & T);
barycentric Tx1(x1 & T);
if (debug)
{
@ -718,7 +751,7 @@ Foam::scalar Foam::particle::trackToStationaryTri
scalar muH = std::isnormal(detA) && detA <= 0 ? VGREAT : 1/detA;
for (label i = 0; i < 4; ++ i)
{
if (std::isnormal(Tx1[i]) && Tx1[i] < 0)
if (Tx1[i] < - detA*SMALL)
{
scalar mu = - y0[i]/Tx1[i];
@ -776,6 +809,30 @@ Foam::scalar Foam::particle::trackToStationaryTri
// Set the proportion of the track that has been completed
stepFraction_ += fraction*muH*detA;
if (debug)
{
Info << "Step Fraction : " << stepFraction_ << endl;
Info << "fraction*muH*detA : " << fraction*muH*detA << endl;
Info << "static muH : " << muH << endl;
Info << "origId() : " << origId() << endl;
}
// Accumulate displacement behind
if (detA <= 0 || nBehind_ > 0)
{
behind_ += muH*detA*mag(displacement);
if (behind_ > 0)
{
behind_ = 0;
nBehind_ = 0;
}
else
{
++ nBehind_;
}
}
return iH != -1 ? 1 - muH*detA : 0;
}
@ -803,12 +860,20 @@ Foam::scalar Foam::particle::trackToMovingTri
FixedList<barycentricTensor, 3> T;
movingTetReverseTransform(fraction, centre, detA, T);
// Get the factor by which the displacement is increased
const scalar f = detA[0] >= 0 ? 1 : negativeSpaceDisplacementFactor;
if (debug)
{
Pair<vector> o, b, v1, v2;
movingTetGeometry(fraction, o, b, v1, v2);
Info<< "Tet points o=" << o[0] << ", b=" << b[0]
<< ", v1=" << v1[0] << ", v2=" << v2[0] << endl
<< "Tet determinant = " << detA[0] << endl
<< "Start local coordinates = " << y0[0] << endl;
}
// Get the relative global position
const vector x0Rel = x0 - centre[0];
const vector x1Rel = f*x1 - centre[1];
const vector x1Rel = x1 - centre[1];
// Form the determinant and hit equations
cubicEqn detAEqn(sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
@ -826,6 +891,16 @@ Foam::scalar Foam::particle::trackToMovingTri
hitEqn[i] = cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
}
if (debug)
{
for (label i = 0; i < 4; ++ i)
{
Info<< (i ? " " : "Hit equation ") << i << " = "
<< hitEqn[i] << endl;
}
Info<< " DetA equation = " << detA << endl;
}
// Calculate the hit fraction
label iH = -1;
scalar muH = std::isnormal(detA[0]) && detA[0] <= 0 ? VGREAT : 1/detA[0];
@ -835,8 +910,33 @@ Foam::scalar Foam::particle::trackToMovingTri
for (label j = 0; j < 3; ++j)
{
if (mu.type(j) == roots::real && hitEqn[i].derivative(mu[j]) < 0)
if
(
mu.type(j) == roots::real
&& hitEqn[i].derivative(mu[j]) < - detA[0]*SMALL
)
{
if (debug)
{
const barycentric yH
(
hitEqn[0].value(mu[j]),
hitEqn[1].value(mu[j]),
hitEqn[2].value(mu[j]),
hitEqn[3].value(mu[j])
);
const scalar detAH = detAEqn.value(mu[j]);
Info<< "Hit on tet face " << i << " at local coordinate "
<< (std::isnormal(detAH) ? name(yH/detAH) : "???")
<< ", " << mu[j]*detA[0]*100 << "% of the "
<< "way along the track" << endl;
Info<< "derivative : " << hitEqn[i].derivative(mu[j]) << nl
<< " coord " << j << " mu[j]: " << mu[j] << nl
<< " hitEq " << i << endl;
}
if (0 <= mu[j] && mu[j] < muH)
{
iH = i;
@ -887,8 +987,26 @@ Foam::scalar Foam::particle::trackToMovingTri
coordinates_ = yH;
tetTriI = iH;
scalar advance = muH*detA[0];
// Set the proportion of the track that has been completed
stepFraction_ += fraction*muH*detA[0];
stepFraction_ += fraction*advance;
// Accumulate displacement behind
if (detA[0] <= 0 || nBehind_ > 0)
{
behind_ += muH*detA[0]*mag(displacement);
if (behind_ > 0)
{
behind_ = 0;
nBehind_ = 0;
}
else
{
++ nBehind_;
}
}
if (debug)
{
@ -900,10 +1018,15 @@ Foam::scalar Foam::particle::trackToMovingTri
{
Pout<< "Track hit no tet faces" << endl;
}
Pout<< "End local coordinates = " << yH << endl
<< "End global coordinates = " << position() << endl;
// Pout<< "End local coordinates = " << yH << endl
// << "End global coordinates = " << position() << endl
// << "Tracking displacement = " << position() - x0 << endl
// << muH*detA[0]*100 << "% of the step from " << stepFraction_
// << " to " << stepFraction_ + fraction << " completed" << endl
// << endl;
}
return iH != -1 ? 1 - muH*detA[0] : 0;
}
@ -915,7 +1038,7 @@ Foam::scalar Foam::particle::trackToTri
label& tetTriI
)
{
if (mesh_.moving())
if ((mesh_.moving() && (stepFraction_ != 1 || fraction != 0)))
{
return trackToMovingTri(displacement, fraction, tetTriI);
}
@ -1049,7 +1172,7 @@ void Foam::particle::correctAfterInteractionListReferral(const label celli)
// so this approximate topology is good enough. By using the nearby cell we
// minimise the error associated with the incorrect topology.
coordinates_ = barycentric(1, 0, 0, 0);
if (mesh_.moving())
if (mesh_.moving() && stepFraction_ != 1)
{
Pair<vector> centre;
FixedList<scalar, 4> detA;

44
src/lagrangian/basic/particle/particle.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2017-2019 OpenCFD Ltd.
Copyright (C) 2017-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -78,7 +78,7 @@ class particle
:
public IDLList<particle>::link
{
// Private member data
// Private Data
//- Size in bytes of the position data
static const std::size_t sizeofPosition;
@ -86,15 +86,9 @@ class particle
//- Size in bytes of the fields
static const std::size_t sizeofFields;
//- The factor by which the displacement is increased when passing
// through negative space. This should be slightly bigger than one.
// This is needed as a straight trajectory can form a closed loop
// through regions of overlapping positive and negative space, leading
// to a track which never ends. By increasing the rate of displacement
// through negative regions, the change in track fraction over this
// part of the loop no longer exactly cancels the change over the
// positive part, and the track therefore terminates.
static const scalar negativeSpaceDisplacementFactor;
//- The value of nBehind_ at which tracking is abandoned. See the
// description of nBehind_.
static const label maxNBehind_;
public:
@ -103,7 +97,7 @@ public:
{
public:
// Public data
// Public Data
//- Flag to switch processor
bool switchProcessor;
@ -135,7 +129,7 @@ public:
private:
// Private data
// Private Data
//- Reference to the polyMesh database
const polyMesh& mesh_;
@ -147,12 +141,12 @@ private:
label celli_;
//- Index of the face that owns the decomposed tet that the
// particle is in
//- particle is in
label tetFacei_;
//- Index of the point on the face that defines the decomposed
// tet that the particle is in. Relative to the face base
// point.
//- tet that the particle is in. Relative to the face base
//- point.
label tetPti_;
//- Face index if the particle is on a face otherwise -1
@ -161,6 +155,18 @@ private:
//- Fraction of time-step completed
scalar stepFraction_;
//- The distance behind the maximum distance reached so far
scalar behind_;
//- The number of tracks carried out that ended in a distance behind the
//- maximum distance reached so far. Once this reaches maxNBehind_,
// tracking is abandoned for the current step. This is needed because
// when tetrahedra are inverted a straight trajectory can form a closed
// loop through regions of overlapping positive and negative space.
// Without this break clause, such loops can result in a valid track
// which never ends.
label nBehind_;
//- Originating processor id
label origProc_;
@ -331,7 +337,7 @@ protected:
public:
// Static data members
// Static Data Members
//- Runtime type information
TypeName("particle");
@ -526,6 +532,9 @@ public:
//- Return current particle position
inline vector position() const;
//- Reset particle data
inline void reset();
// Track
@ -664,6 +673,7 @@ public:
void relocate(const point& position, const label celli = -1);
// I-O
//- Write the name representation to stream

21
src/lagrangian/basic/particle/particleI.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2018 OpenFOAM Foundation
Copyright (C) 2011-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -68,6 +69,7 @@ inline void Foam::particle::movingTetGeometry
) const
{
const triFace triIs(currentTetIndices().faceTriIs(mesh_));
const pointField& ptsOld = mesh_.oldPoints();
const pointField& ptsNew = mesh_.points();
@ -75,8 +77,11 @@ inline void Foam::particle::movingTetGeometry
// we need to put a mesh_.oldCellCentres() method in for this to work. The
// values obtained from the mesh and those obtained from the cell do not
// necessarily match. See mantis #1993.
const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
//const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
//const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
const vector ccOld = mesh_.oldCellCentres()[celli_];
const vector ccNew = mesh_.cellCentres()[celli_];
// Old and new points and cell centres are not sub-cycled. If we are sub-
// cycling, then we have to account for the timestep change here by
@ -265,7 +270,7 @@ inline Foam::tetIndices Foam::particle::currentTetIndices() const
inline Foam::barycentricTensor Foam::particle::currentTetTransform() const
{
if (mesh_.moving())
if (mesh_.moving() && stepFraction_ != 1)
{
return movingTetTransform(0)[0];
}
@ -312,6 +317,14 @@ inline Foam::vector Foam::particle::position() const
}
inline void Foam::particle::reset()
{
stepFraction_ = 0;
nBehind_ = 0;
behind_ = 0;
}
void Foam::particle::patchData(vector& n, vector& U) const
{
if (!onBoundaryFace())
@ -321,7 +334,7 @@ void Foam::particle::patchData(vector& n, vector& U) const
<< exit(FatalError);
}
if (mesh_.moving())
if ((mesh_.moving() && stepFraction_ != 1))
{
Pair<vector> centre, base, vertex1, vertex2;
movingTetGeometry(1, centre, base, vertex1, vertex2);

55
src/lagrangian/basic/particle/particleTemplates.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2016-2019 OpenCFD Ltd.
Copyright (C) 2016-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -293,7 +293,7 @@ void Foam::particle::writeObjects(const CloudType& c, objectRegistry& obr)
template<class TrackCloudType>
void Foam::particle::hitFace
(
const vector& direction,
const vector& displacement,
TrackCloudType& cloud,
trackingData& td
)
@ -337,11 +337,11 @@ void Foam::particle::hitFace
}
else if (isA<cyclicACMIPolyPatch>(patch))
{
p.hitCyclicACMIPatch(cloud, ttd, direction);
p.hitCyclicACMIPatch(cloud, ttd, displacement);
}
else if (isA<cyclicAMIPolyPatch>(patch))
{
p.hitCyclicAMIPatch(cloud, ttd, direction);
p.hitCyclicAMIPatch(cloud, ttd, displacement);
}
else if (isA<processorPolyPatch>(patch))
{
@ -459,7 +459,7 @@ void Foam::particle::hitCyclicAMIPatch
(
TrackCloudType&,
trackingData& td,
const vector& direction
const vector& displacement
)
{
vector pos = position();
@ -468,7 +468,14 @@ void Foam::particle::hitCyclicAMIPatch
static_cast<const cyclicAMIPolyPatch&>(mesh_.boundaryMesh()[patch()]);
const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch();
const label sendFacei = cpp.whichFace(facei_);
const label receiveFacei = cpp.pointFace(sendFacei, direction, pos);
const label receiveFacei = cpp.pointFace(sendFacei, displacement, pos);
if (false)
{
Info<< "My new pos : " << pos << endl;
Info<< "Particle " << origId() << " crossing AMI from " << cpp.name()
<< " to " << receiveCpp.name() << endl << endl;
}
if (receiveFacei < 0)
{
@ -485,12 +492,13 @@ void Foam::particle::hitCyclicAMIPatch
facei_ = tetFacei_ = receiveFacei + receiveCpp.start();
// Locate the particle on the receiving side
vector directionT = direction;
cpp.reverseTransformDirection(directionT, sendFacei);
vector displacementT = displacement;
cpp.reverseTransformDirection(displacementT, sendFacei);
locate
(
pos,
&directionT,
&displacementT,
mesh_.faceOwner()[facei_],
false,
"Particle crossed between " + cyclicAMIPolyPatch::typeName +
@ -523,6 +531,27 @@ void Foam::particle::hitCyclicAMIPatch
);
transformProperties(-s);
}
//if (onBoundaryFace())
{
//DebugVar("On boudanry")
// vector receiveNormal, receiveDisplacement;
// patchData(receiveNormal, receiveDisplacement);
//
// if (((displacementT - fraction*receiveDisplacement)&receiveNormal) > 0)
// {
// td.keepParticle = false;
// WarningInFunction
// << "Particle transfer from " << cyclicAMIPolyPatch::typeName
// << " patches " << cpp.name() << " to " << receiveCpp.name()
// << " failed at position " << pos << " and with displacement "
// << (displacementT - fraction*receiveDisplacement) << nl
// << " The displacement points into both the source and "
// << "receiving faces, so the tracking cannot proceed" << nl
// << " The particle has been removed" << nl << endl;
// return;
// }
}
}
@ -531,7 +560,7 @@ void Foam::particle::hitCyclicACMIPatch
(
TrackCloudType& cloud,
trackingData& td,
const vector& direction
const vector& displacement
)
{
const cyclicACMIPolyPatch& cpp =
@ -551,20 +580,20 @@ void Foam::particle::hitCyclicACMIPatch
if (!couple && !nonOverlap)
{
vector pos = position();
couple = cpp.pointFace(localFacei, direction, pos) >= 0;
couple = cpp.pointFace(localFacei, displacement, pos) >= 0;
nonOverlap = !couple;
}
if (couple)
{
hitCyclicAMIPatch(cloud, td, direction);
hitCyclicAMIPatch(cloud, td, displacement);
}
else
{
// Move to the face associated with the non-overlap patch and redo the
// face interaction.
tetFacei_ = facei_ = cpp.nonOverlapPatch().start() + localFacei;
hitFace(direction, cloud, td);
hitFace(displacement, cloud, td);
}
}