Commit f1780abe authored by Will Bainbridge's avatar Will Bainbridge Committed by Andrew Heather
Browse files

lagrangian: Optimised reduced-dimension constraints

The 4.x tracking enforces reduced dimensionality on the parcels by
moving them to the centre of the mesh at the start of each track,
without considering the topology. This can leave the parcel outside it's
associated tetrahedron.

The barycentric algorithm isn't tolerant to incorrect topology, so
instead of changing position, it was written to track to the mesh
centre. This worked, but effectively doubled the number of tracking
calls. This additional cost has now been removed by absorbing the
constraint displacement into the existing motion track, so that the same
number of tracking steps are performed as before.

Partially resolves bug report https://bugs.openfoam.org/view.php?id=2688
parent 9ab570ae
......@@ -53,17 +53,17 @@ bool Foam::DSMCParcel<ParcelType>::move
while (td.keepParticle && !td.switchProcessor && p.stepFraction() < 1)
{
// Apply correction to position for reduced-D cases
p.constrainToMeshCentre();
Utracking = U_;
// Apply correction to velocity to constrain tracking for
// reduced-D cases
meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking);
// Deviation from the mesh centre for reduced-D cases
const vector d = p.deviationFromMeshCentre();
const scalar f = 1 - p.stepFraction();
p.trackToAndHitFace(f*trackTime*Utracking, f, cloud, td);
p.trackToAndHitFace(f*trackTime*Utracking - d, f, cloud, td);
}
return td.keepParticle;
......
......@@ -984,25 +984,17 @@ Foam::scalar Foam::particle::trackToTri
}
void Foam::particle::constrainToMeshCentre()
Foam::vector Foam::particle::deviationFromMeshCentre() const
{
const Vector<label>& dirs = mesh_.geometricD();
bool isConstrained = false;
forAll(dirs, dirI)
if (cmptMin(mesh_.geometricD()) == -1)
{
if (dirs[dirI] == -1)
{
isConstrained = true;
break;
}
vector pos = position(), posC = pos;
meshTools::constrainToMeshCentre(mesh_, posC);
return pos - posC;
}
if (isConstrained)
else
{
vector pos = position();
meshTools::constrainToMeshCentre(mesh_, pos);
track(pos - position(), 0);
return vector::zero;
}
}
......
......@@ -598,9 +598,10 @@ public:
trackingData& td
);
//- Set the constrained components of the particle position to the
// mesh centre.
void constrainToMeshCentre();
//- Get the displacement from the mesh centre. Used to correct the
// particle position in cases with reduced dimensionality. Returns a
// zero vector for three-dimensional cases.
vector deviationFromMeshCentre() const;
// Patch data
......
......@@ -276,9 +276,6 @@ bool Foam::KinematicParcel<ParcelType>::move
while (ttd.keepParticle && !ttd.switchProcessor && p.stepFraction() < 1)
{
// Apply correction to position for reduced-D cases
p.constrainToMeshCentre();
// Cache the current position, cell and step-fraction
const point start = p.position();
const scalar sfrac = p.stepFraction();
......@@ -289,6 +286,9 @@ bool Foam::KinematicParcel<ParcelType>::move
// Cell length scale
const scalar l = cellLengthScale[p.cell()];
// Deviation from the mesh centre for reduced-D cases
const vector d = p.deviationFromMeshCentre();
// Fraction of the displacement to track in this loop. This is limited
// to ensure that the both the time and distance tracked is less than
// maxCo times the total value.
......@@ -298,7 +298,7 @@ bool Foam::KinematicParcel<ParcelType>::move
if (p.active())
{
// Track to the next face
p.trackToFace(f*s, f);
p.trackToFace(f*s - d, f);
}
else
{
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment