Regression in extractEulerianParticles

Summary

The extractEulerianParticles function object reports a much larger number of droplets/particles than actually present.

Steps to reproduce

Run the attached small test case on master (commit e651d635).

Example case

Box domain, two rows of droplets (3 + 2) initialized, constant velocity in -z direction (adaptation of a single-droplet case by Sptizenberger et al., 2020).

extractFuncDroplets.tar.gz

What is the current bug behaviour?

The extractEulerianParticles function object reports 160 particles instead of easily manually identifiable five. Consequently, the particle data in <time>/lagrangian/eulerianParticleCloud/* is incorrect, where * are files: coordinates d origId origProcId positions soi tag U.

What is the expected correct behavior?

The five identified particles with correct statistics in <time>/lagrangian/eulerianParticleCloud/*.

e.g the volumes of the 5 particles are 5(5.01041e-05 5.05945e-05 5.01041e-05 5.14354e-05 5.14354e-05)

Relevant logs and/or images

in the log.interFoam at the last timestep, the correct behaviour is:

extractEulerianParticles extractEulerianParticles output:
    Collected particles   : 5
    Collected volume      : 3.42031e-13
    Discarded particles   : 0
    Discarded volume      : 0
    Particles in progress : 0

current behaviour is:

extractEulerianParticles extractEulerianParticles output:
    Collected particles   : 160
    Collected volume      : 3.42031e-13
    Discarded particles   : 0
    Discarded volume      : 0
    Particles in progress : 0

Environment information

• OpenFOAM version : master (commit e651d635)
• Operating system : Ubuntu 22.04
• Compiler : gcc, Open MPI 4.1.2

Possible fixes

I tracked the issue down to a regression in regionSplit2D.C where the parchEdgeFaceRegion was replaced by the edgeTopoDistanceData, more specifically in the call to

PatchEdgeFaceWave
   <
       indirectPrimitivePatch,
       edgeTopoDistanceData<label>
   >
   (
       mesh,
       patch,
       changedEdges,
       changedRegions,
       allEdgeInfo,
       allFaceInfo,
       returnReduce(patch.nEdges(), sumOp<label>())
   );

which results in numerous member data updates, one of which is the following method call to update the face data

template
<
   class PrimitivePatchType,
   class Type,
   class TrackingData
>
bool Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
updateFace
(
   const label facei,
   const label neighbourEdgeI,
   const Type& neighbourInfo,
   Type& faceInfo
)
{
   nEvals_++;


   bool wasValid = faceInfo.valid(td_);


   bool propagate =
       faceInfo.updateFace
       (
           mesh_,
           patch_,
           facei,
           neighbourEdgeI,
           neighbourInfo,
           propagationTol_,
           td_
       );

which, in turn, calls the templated type’s faceInfo.updateFace(...) member.

The trouble is that parchEdgeFaceRegion and edgeTopoDistanceData implement different logic in the updateFace(...) member.

In edgeTopoDistanceData the region_ member is operated on:

template<class TrackingData>
inline bool Foam::patchEdgeFaceRegion::updateFace
(
   const polyMesh& mesh,
   const indirectPrimitivePatch& patch,
   const label facei,
   const label edgeI,
   const patchEdgeFaceRegion& edgeInfo,
   const scalar tol,
   TrackingData& td
)
{
   return update(edgeInfo, tol, td);
}

which calls

// Update this with w2 if w2 nearer to pt.
template<class TrackingData>
inline bool Foam::patchEdgeFaceRegion::update
(
   const patchEdgeFaceRegion& w2,
   const scalar tol,
   TrackingData& td
)
{
   if (!w2.valid(td))
   {
       FatalErrorInFunction
           << "problem." << abort(FatalError);
   }


   if (w2.region_ == -2 || region_ == -2)
   {
       // Blocked edge/face
       return false;
   }


   if (!valid(td))
   {
       // current not yet set so use any value
       operator=(w2);
       return true;
   }
   else
   {
       if (w2.region_ < region_)
       {
           operator=(w2);
           return true;
       }
       else
       {
           return false;
       }
   }
}

The equivalent of region_ in edgeTopoDistanceData would be data_, however the function

template<class Type, class PrimitivePatchType>
template<class TrackingData>
inline bool Foam::edgeTopoDistanceData<Type, PrimitivePatchType>::updateFace
(
   const polyMesh& mesh,
   const PrimitivePatchType& patch,
   const label facei,
   const label edgeI,
   const edgeTopoDistanceData<Type, PrimitivePatchType>& edgeInfo,
   const scalar tol,
   TrackingData& td
)
{
   // From edge to face
   if (distance_ == -1)
   {
       this->operator=(edgeInfo);
       return true;
   }


   return false;
}

only checks for distance_. I suppose this is the logic needed in fluxSummary, meshRefinementBaffles and other places where edgeTopoDistanceData is used. In contrast, patchEdgeFaceRegion doesn’t seem to be used anywhere else. The regression removed the only usage of it (while not deleting it from the repository entirely) which might suggest that it was unintentional (?).

When I made the change in regionSplit2D to once again use the patchEdgeFaceRegion, the issue was resolved. I have a branch with these changes, which I'd be happy to share with you if given push rights to the repo.