Bugfix: Parallel Lagrangian (2x) & parallel point interpolation. Henrik Rusche

src/foam/fields/PointPatchFieldTemplates/constraint/processor/ProcessorPointPatchField.C

@@ -84,7 +84,7 @@ sendField
 
     //HJ: This needs complete rewrite:
     // - move communications into a patch
-    // - allow for various types of communication
+    // - allow for various types of communication - done HR, 12/6/2017
     // HJ, 15/Apr/2009
 
     if (commsType == Pstream::blocking || commsType == Pstream::scheduled)
@@ -101,6 +101,7 @@ sendField
     {
         resizeBuf(receiveBuf_, f.size()*sizeof(Type));
 
+        outstandingRecvRequest_ = Pstream::nRequests();
         IPstream::read
         (
             commsType,
@@ -112,6 +113,7 @@ sendField
         resizeBuf(sendBuf_, f.byteSize());
         memcpy(sendBuf_.begin(), f.begin(), f.byteSize());
 
+        outstandingSendRequest_ = Pstream::nRequests();
         OPstream::write
         (
             commsType,
@@ -127,22 +129,6 @@ sendField
             << exit(FatalError);
     }
 
-    // Not using non-blocking comms
-//     if (commsType == Pstream::nonBlocking)
-//     {
-//         FatalErrorIn("void ProcessorPointPatchField::sendField")
-//             << "Non-blocking comms not implemented"
-//             << abort(FatalError);
-//     }
-
-//     OPstream::write
-//     (
-//         commsType,
-//         procPatch_.neighbProcNo(),
-//         reinterpret_cast<const char*>(f.begin()),
-//         f.byteSize()
-//     );
-
     tf.clear();
 }
 
@@ -167,6 +153,27 @@ receivePointField
 {
     tmp<Field<Type2> > tf(new Field<Type2>(this->size()));
 
+    if (Pstream::parRun())
+    {
+        if (commsType == Pstream::nonBlocking)
+        {
+            // Receive into tf
+
+            if
+            (
+                outstandingRecvRequest_ >= 0
+             && outstandingRecvRequest_ < Pstream::nRequests()
+            )
+            {
+                Pstream::waitRequest(outstandingRecvRequest_);
+            }
+            outstandingSendRequest_ = -1;
+            outstandingRecvRequest_ = -1;
+
+            memcpy(tf().begin(), receiveBuf_.begin(), tf().byteSize());
+        }
+        else
+        {
             IPstream::read
             (
                 commsType,
@@ -174,6 +181,8 @@ receivePointField
                 reinterpret_cast<char*>(tf().begin()),
                 tf().byteSize()
             );
+        }
+    }
 
     return tf;
 }
@@ -202,6 +211,27 @@ receiveEdgeField
         new Field<Type2>(procPatch_.localEdgeIndices().size())
     );
 
+    if (Pstream::parRun())
+    {
+        if (commsType == Pstream::nonBlocking)
+        {
+            // Receive into tf
+
+            if
+            (
+                outstandingRecvRequest_ >= 0
+             && outstandingRecvRequest_ < Pstream::nRequests()
+            )
+            {
+                Pstream::waitRequest(outstandingRecvRequest_);
+            }
+            outstandingSendRequest_ = -1;
+            outstandingRecvRequest_ = -1;
+
+            memcpy(tf().begin(), receiveBuf_.begin(), tf().byteSize());
+        }
+        else
+        {
             IPstream::read
             (
                 commsType,
@@ -209,6 +239,8 @@ receiveEdgeField
                 reinterpret_cast<char*>(tf().begin()),
                 tf().byteSize()
             );
+        }
+    }
 
     return tf;
 }
@@ -470,7 +502,7 @@ initEvaluate
     {
         if (this->isPointField())
         {
-            initAddFieldTempl(Pstream::blocking, this->internalField());
+            initAddFieldTempl(commsType, this->internalField());
         }
     }
 }
@@ -551,7 +583,7 @@ ProcessorPointPatchField
 <PatchField, Mesh, PointPatch, ProcessorPointPatch, MatrixType, Type>::
 initAddField() const
 {
-    initAddFieldTempl(Pstream::blocking, this->internalField());
+    initAddFieldTempl(Pstream::defaultComms(), this->internalField());
 }
 
 
@@ -570,7 +602,7 @@ ProcessorPointPatchField
 <PatchField, Mesh, PointPatch, ProcessorPointPatch, MatrixType, Type>::
 addField(Field<Type>& f) const
 {
-    addFieldTempl(Pstream::blocking, f);
+    addFieldTempl(Pstream::defaultComms(), f);
 }
 
 
@@ -636,7 +668,7 @@ ProcessorPointPatchField
 <PatchField, Mesh, PointPatch, ProcessorPointPatch, MatrixType, Type>::
 initAddDiag(const scalarField& d) const
 {
-    initAddFieldTempl(Pstream::blocking, d);
+    initAddFieldTempl(Pstream::defaultComms(), d);
 }
 
 
@@ -655,7 +687,7 @@ ProcessorPointPatchField
 <PatchField, Mesh, PointPatch, ProcessorPointPatch, MatrixType, Type>::
 initAddSource(const scalarField& s) const
 {
-    initAddFieldTempl(Pstream::blocking, s);
+    initAddFieldTempl(Pstream::defaultComms(), s);
 }
 
 
@@ -674,7 +706,7 @@ ProcessorPointPatchField
 <PatchField, Mesh, PointPatch, ProcessorPointPatch, MatrixType, Type>::
 addDiag(scalarField& d) const
 {
-    addFieldTempl(Pstream::blocking, d);
+    addFieldTempl(Pstream::defaultComms(), d);
 }
 
 
@@ -693,7 +725,7 @@ ProcessorPointPatchField
 <PatchField, Mesh, PointPatch, ProcessorPointPatch, MatrixType, Type>::
 addSource(scalarField& s) const
 {
-    addFieldTempl(Pstream::blocking, s);
+    addFieldTempl(Pstream::defaultComms(), s);
 }
 
 

src/foam/fields/PointPatchFieldTemplates/constraint/processor/ProcessorPointPatchField.H

@@ -75,6 +75,12 @@ class ProcessorPointPatchField
         // Non-blocking parallel communications
         // Temporary: move to patch.  HJ, 15/Apr/2008
 
+            //- Outstanding request
+            mutable label outstandingSendRequest_;
+
+            //- Outstanding request
+            mutable label outstandingRecvRequest_;
+
             //- Send buffer.
             //  Only sized and used when compressed or non-blocking comms used.
             mutable List<char> sendBuf_;

src/foam/meshes/primitiveMesh/primitiveMeshFindCell.C

@@ -48,8 +48,6 @@ bool Foam::primitiveMesh::pointInCell(const point& p, label celli) const
     const vectorField& cf = faceCentres();
     const vectorField& Sf = faceAreas();
 
-    bool inCell = true;
-
     forAll(f, facei)
     {
         label nFace = f[facei];
@@ -59,10 +57,14 @@ bool Foam::primitiveMesh::pointInCell(const point& p, label celli) const
         {
             normal = -normal;
         }
-        inCell = inCell && ((normal & proj) <= 0);
+
+        if ((normal & proj) > 0)
+        {
+            return false;
+        }
     }
 
-    return inCell;
+    return true;
 }
 
 

src/lagrangian/basic/Particle/ParticleI.H

@@ -37,6 +37,7 @@ inline Foam::scalar Foam::Particle<ParticleType>::lambda
 ) const
 {
     const polyMesh& mesh = cloud_.polyMesh_;
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
     bool movingMesh = mesh.moving();
 
     if (movingMesh)
@@ -46,7 +47,11 @@ inline Foam::scalar Foam::Particle<ParticleType>::lambda
         vector Cf = mesh.faceCentres()[facei];
 
         // move reference point for wall
-        if (!cloud_.internalFace(facei))
+        if
+        (
+            !cloud_.internalFace(facei)
+         && isA<wallPolyPatch>(patches[cloud_.facePatch(facei)])
+        )
         {
             const vector& C = mesh.cellCentres()[celli_];
             scalar CCf = mag((C - Cf) & Sf);
@@ -184,13 +189,18 @@ inline Foam::scalar Foam::Particle<ParticleType>::lambda
 ) const
 {
     const polyMesh& mesh = cloud_.polyMesh_;
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
     vector Sf = mesh.faceAreas()[facei];
     Sf /= mag(Sf);
     vector Cf = mesh.faceCentres()[facei];
 
     // move reference point for wall
-    if (!cloud_.internalFace(facei))
+    if
+    (
+        !cloud_.internalFace(facei)
+     && isA<wallPolyPatch>(patches[cloud_.facePatch(facei)])
+    )
     {
         const vector& C = mesh.cellCentres()[celli_];
         scalar CCf = mag((C - Cf) & Sf);

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C

@@ -137,6 +137,8 @@ void Foam::InjectionModel<CloudType>::findCellAtPosition
 
     bool foundCell = false;
 
+    // Force creation of mesh.C() to avoid dead-lock in comms
+    owner_.mesh().C();
     cellI = owner_.mesh().findCell(position);
 
     if (cellI >= 0)