Bugfix: future patch processor ordering

src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C

@@ -126,6 +126,176 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
         List<SLList<label> > procPatchIndex(nProcs_);
 
+        // Rewrite:
+        // Handling of coupled patches is changed.  HJ, 11/Apr/2018
+
+        // Prepare collection of patch faces
+
+        // For all processors, set the size of start index and patch size
+        // lists to the number of patches in the mesh
+        forAll (procPatchSize_, procI)
+        {
+            procPatchSize_[procI].setSize(patches.size());
+            procPatchStartIndex_[procI].setSize(patches.size());
+        }
+
+        forAll (patches, patchI)
+        {
+            // Reset size for all processors
+            forAll (procPatchSize_, procI)
+            {
+                procPatchSize_[procI][patchI] = 0;
+            }
+        }
+
+
+        // Algorithm:
+        // When running the decomposition in parallel, it is assumed that
+        // the result will be used in dynamic load balancing
+        // For load balancing, the processor patches need to match in order for
+        // the patch-to-patch matching to work propely on mesh reconstruction
+        // Therefore, all processor patches need to be split into the matching
+        // and non-matching part be examining the cellToProc_ data on
+        // the neighbour side.  This has been prepared in patchNbrCellToProc_
+        // HJ, 11/Apr/2018
+
+        // The correct order of dumping the faces is to dump the
+        // slave processor patches first, then internal faces and then
+        // master processor patches
+        // HJ, 23/Apr/2018
+
+        // Dump current processor patch faces into new processor patches
+        // that will be created in decomposition when running in parallel
+        forAll (patches, patchI)
+        {
+            // Check the processor patch for which neighbour data exists
+            if
+            (
+                isA<processorFvPatch>(patches[patchI])
+             && !patchNbrCellToProc_[patchI].empty()
+            )
+            {
+                const processorFvPatch& procPatch =
+                    refCast<const processorFvPatch>(patches[patchI]);
+
+                // DO ONLY SLAVE SIDE
+                if (!procPatch.master())
+                {
+                    // Get patch start
+                    const label patchStart = patches[patchI].patch().start();
+
+                    // Get faceCells
+                    const labelList& fc = patches[patchI].faceCells();
+
+                    // Get neighbour cellToProc addressing across the interface
+                    const labelList& curNbrPtc = patchNbrCellToProc_[patchI];
+
+                    forAll (fc, patchFaceI)
+                    {
+                        // Local owner proc is looked up using faceCells
+                        const label ownerProc = cellToProc_[fc[patchFaceI]];
+
+                        // Neighbour proc is looked up directly
+                        const label neighbourProc = curNbrPtc[patchFaceI];
+
+                        // Check change in processor type across the processor
+                        // boundary
+                        // If ownerProc and neighbourProc are the same, the
+                        // processor face will be merged, meaning that it
+                        // remains in the (old) processor patch
+                        // If ownerProc and neighbourProc are different,
+                        // this will be a new processor boundary created from
+                        //  the existing processor face and added afterwards
+                        if (ownerProc != neighbourProc)
+                        {
+                            Pout<< "Slave: New proc from old proc["
+                                << patchFaceI
+                                << "]: " << ownerProc << " " << neighbourProc
+                                << " face: " << patchStart + patchFaceI
+                                << endl;
+
+                            // Search algorithm repeated in processor patches.
+                            // HJ, 11/Apr/2018
+                            SLList<label>::iterator curInterProcBdrsOwnIter =
+                                interProcBoundaries[ownerProc].begin();
+
+                            SLList<SLList<label> >::iterator
+                                curInterProcBFacesOwnIter =
+                                interProcBFaces[ownerProc].begin();
+
+                            bool interProcBouFound = false;
+
+                            // WARNING: Synchronous SLList iterators
+
+                            for
+                            (
+                                ;
+                                curInterProcBdrsOwnIter
+                             != interProcBoundaries[ownerProc].end()
+                             && curInterProcBFacesOwnIter
+                             != interProcBFaces[ownerProc].end();
+                                ++curInterProcBdrsOwnIter,
+                                ++curInterProcBFacesOwnIter
+                            )
+                            {
+                                if (curInterProcBdrsOwnIter() == neighbourProc)
+                                {
+                                    // Inter - processor boundary exists.
+                                    // Add the face
+                                    interProcBouFound = true;
+
+                                    // Add to owner only!
+                                    curInterProcBFacesOwnIter().append
+                                    (
+                                        patchStart + patchFaceI
+                                    );
+                                }
+
+                                if (interProcBouFound) break;
+                            }
+
+                            if (!interProcBouFound)
+                            {
+                                // inter - processor boundaries do not exist and
+                                // need to be created
+
+                                // set the new addressing information
+
+                                // Add to owner only!
+                                interProcBoundaries[ownerProc].append
+                                (
+                                    neighbourProc
+                                );
+                                interProcBFaces[ownerProc].append
+                                (
+                                    SLList<label>(patchStart + patchFaceI)
+                                );
+                            }
+                        }
+                        else
+                        {
+                            // Owner and neighbour processor index are the same
+                            // across the processor boundary.  This face will be
+                            // re-merged into internal faces in load balancing
+                            // and therefore remains in the processor patch
+                            Pout<< "Preserved proc[" << patchFaceI << "]: "
+                                << ownerProc << " " << neighbourProc << endl;
+
+                            // Add the face
+                            procFaceList[ownerProc].append
+                            (
+                                patchStart + patchFaceI
+                            );
+
+                            // Increment the number of faces for this patch
+                            procPatchSize_[ownerProc][patchI]++;
+                        }
+                    }
+                }
+            }
+        }
+
+        // Internal mesh faces
         forAll (neighbour, faceI)
         {
             if (cellToProc_[owner[faceI]] != cellToProc_[neighbour[faceI]])
@@ -199,15 +369,19 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                             if (neighbourFound) break;
                         }
 
-                        if (interProcBouFound && !neighbourFound)
+                        if (!neighbourFound)
                         { 
-                            FatalErrorIn
+                            // Owner found (slave proc boundary), but neighbour
+                            // not found: add it
+                            interProcBoundaries[neighbourProc].append
                             (
-                                "domainDecomposition::decomposeMesh()"
+                                ownerProc
-                            )   << "Inconsistency in inter - "
+                            );
-                                << "processor boundary lists for processors "
+
-                                << ownerProc << " and " << neighbourProc
+                            interProcBFaces[neighbourProc].append
-                                << abort(FatalError);
+                            (
+                                SLList<label>(faceI)
+                            );
                         }
                     }
 
@@ -221,11 +395,11 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
                     // set the new addressing information
 
-                    // owner
+                    // Add owner side
                     interProcBoundaries[ownerProc].append(neighbourProc);
                     interProcBFaces[ownerProc].append(SLList<label>(faceI));
 
-                    // neighbour
+                    // Add neighbour side
                     interProcBoundaries[neighbourProc].append(ownerProc);
                     interProcBFaces[neighbourProc].append
                     (
@@ -235,39 +409,6 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
             }
         }
 
-        // Rewrite:
-        // Handling of coupled patches is changed.  HJ, 11/Apr/2018
-
-        // Prepare collection of patch faces
-
-        // For all processors, set the size of start index and patch size
-        // lists to the number of patches in the mesh
-        forAll (procPatchSize_, procI)
-        {
-            procPatchSize_[procI].setSize(patches.size());
-            procPatchStartIndex_[procI].setSize(patches.size());
-        }
-
-        forAll (patches, patchI)
-        {
-            // Reset size for all processors
-            forAll (procPatchSize_, procI)
-            {
-                procPatchSize_[procI][patchI] = 0;
-            }
-        }
-
-
-        // Algorithm:
-        // When running the decomposition in parallel, it is assumed that
-        // the result will be used in dynamic load balancing
-        // For load balancing, the processor patches need to match in order for
-        // the patch-to-patch matching to work propely on mesh reconstruction
-        // Therefore, all processor patches need to be split into the matching
-        // and non-matching part be examining the cellToProc_ data on
-        // the neighbour side.  This has been prepared in patchNbrCellToProc_
-        // HJ, 11/Apr/2018
-
         // Dump current processor patch faces into new processor patches
         // that will be created in decomposition when running in parallel
         forAll (patches, patchI)
@@ -278,6 +419,12 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                 isA<processorFvPatch>(patches[patchI])
              && !patchNbrCellToProc_[patchI].empty()
             )
+            {
+                const processorFvPatch& procPatch =
+                    refCast<const processorFvPatch>(patches[patchI]);
+
+                // DO ONLY MASTER SIDE
+                if (procPatch.master())
                 {
                     // Get patch start
                     const label patchStart = patches[patchI].patch().start();
@@ -299,15 +446,18 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                         // Check change in processor type across the processor
                         // boundary
                         // If ownerProc and neighbourProc are the same, the
-                    // processor face will be merged, meaning that it remains
+                        // processor face will be merged, meaning that it
-                    // in the (old) processor patch
+                        // remains in the (old) processor patch
-                    // If ownerProc and neighbourProc are different, this will
+                        // If ownerProc and neighbourProc are different,
-                    // be a new processor boundary created from the existing
+                        // this will be a new processor boundary created from
-                    // processor face and added afterwards
+                        //  the existing processor face and added afterwards
                         if (ownerProc != neighbourProc)
                         {
-                        Pout<< "New proc from old proc[" << patchFaceI << "]: "
+                            Pout<< "Master: New proc from old proc["
-                            << ownerProc << " " << neighbourProc << endl;
+                                << patchFaceI << "]: "
+                                << ownerProc << " " << neighbourProc
+                                << " face: " << patchStart + patchFaceI << endl;
+
                             // Search algorithm repeated in processor patches.
                             // HJ, 11/Apr/2018
                             SLList<label>::iterator curInterProcBdrsOwnIter =
@@ -338,60 +488,11 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                                     // Add the face
                                     interProcBouFound = true;
 
+                                    // Add to owner only!
                                     curInterProcBFacesOwnIter().append
                                     (
                                         patchStart + patchFaceI
                                     );
-
-                                SLList<label>::iterator
-                                    curInterProcBdrsNeiIter =
-                                    interProcBoundaries[neighbourProc].begin();
-
-                                SLList<SLList<label> >::iterator
-                                    curInterProcBFacesNeiIter =
-                                    interProcBFaces[neighbourProc].begin();
-
-                                bool neighbourFound = false;
-
-                                // WARNING: Synchronous SLList iterators
-
-                                for
-                                (
-                                    ;
-                                    curInterProcBdrsNeiIter !=
-                                    interProcBoundaries[neighbourProc].end()
-                                 && curInterProcBFacesNeiIter !=
-                                    interProcBFaces[neighbourProc].end();
-                                    ++curInterProcBdrsNeiIter,
-                                    ++curInterProcBFacesNeiIter
-                                )
-                                {
-                                    if (curInterProcBdrsNeiIter() == ownerProc)
-                                    {
-                                        // Boundary found. Add the face
-                                        neighbourFound = true;
-
-                                        curInterProcBFacesNeiIter().append
-                                        (
-                                            patchStart + patchFaceI
-                                        );
-                                    }
-
-                                    if (neighbourFound) break;
-                                }
-
-                                if (interProcBouFound && !neighbourFound)
-                                {
-                                    FatalErrorIn
-                                    (
-                                        "domainDecomposition::decomposeMesh()"
-                                    )   << "Inconsistency in inter - "
-                                        << "processor boundary lists for "
-                                        << "processors "
-                                        << ownerProc << " and "
-                                        << neighbourProc
-                                        << abort(FatalError);
-                                }
                                 }
 
                                 if (interProcBouFound) break;
@@ -404,7 +505,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                                 
                                 // set the new addressing information
 
-                            // owner
+                                // Add to owner only!
                                 interProcBoundaries[ownerProc].append
                                 (
                                     neighbourProc
@@ -413,16 +514,6 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                                 (
                                     SLList<label>(patchStart + patchFaceI)
                                 );
-
-                            // neighbour
-                            interProcBoundaries[neighbourProc].append
-                            (
-                                ownerProc
-                            );
-                            interProcBFaces[neighbourProc].append
-                            (
-                                SLList<label>(patchStart + patchFaceI)
-                            );
                             }
                         }
                         else
@@ -434,15 +525,19 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                             Pout<< "Preserved proc[" << patchFaceI << "]: "
                                 << ownerProc << " " << neighbourProc << endl;
 
-                        // add the face
+                            // Add the face
-                        procFaceList[ownerProc].append(patchStart + patchFaceI);
+                            procFaceList[ownerProc].append
+                            (
+                                patchStart + patchFaceI
+                            );
 
-                        // increment the number of faces for this patch
+                            // Increment the number of faces for this patch
                             procPatchSize_[ownerProc][patchI]++;
                         }
                     }
                 }
             }
+        }
 
         // Loop through patches. For cyclic boundaries detect inter-processor
         // faces; for all other, add faces to the face list and remember start
@@ -462,11 +557,13 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
             // Do normal patches.  Note that processor patches
             // have already been done and need to be skipped
             // For all other patches patchNbrCellToProc_ will be empty
-            if
+            if (!patchNbrCellToProc_[patchI].empty())
-            (
+            {
-                !isA<cyclicFvPatch>(patches[patchI])
+                // Processor patch.  Skip it
-             && patchNbrCellToProc_[patchI].empty()
+                continue;
-            )
+            }
+
+            if (!isA<cyclicFvPatch>(patches[patchI]))
             {
                 // Normal patch. Add faces to processor where the cell
                 // next to the face lives

src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C

@@ -39,6 +39,7 @@ License
 
 defineTypeNameAndDebug(Foam::domainDecomposition, 0);
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::domainDecomposition::domainDecomposition
@@ -117,19 +118,19 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
 
     faceList procFaces(curFaceLabels.size());
 
-    forAll (curFaceLabels, facei)
+    forAll (curFaceLabels, faceI)
     {
         // Mark the original face as used
         // Remember to decrement the index by one (turning index)
         // HJ, 5/Dec/2001
-        label curF = mag(curFaceLabels[facei]) - 1;
+        label curF = mag(curFaceLabels[faceI]) - 1;
 
-        faceLookup[curF] = facei;
+        faceLookup[curF] = faceI;
 
         // get the original face
         labelList origFaceLabels;
 
-        if (curFaceLabels[facei] >= 0)
+        if (curFaceLabels[faceI] >= 0)
         {
             // face not turned
             origFaceLabels = meshFaces[curF];
@@ -140,7 +141,7 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
         }
 
         // translate face labels into local point list
-        face& procFaceLabels = procFaces[facei];
+        face& procFaceLabels = procFaces[faceI];
 
         procFaceLabels.setSize(origFaceLabels.size());
 
@@ -275,7 +276,7 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
             procPatches[nPatches] =
                 new passiveProcessorPolyPatch
                 (
-                    word("procBoundary") + Foam::name(procI)
+                    word("passiveProcBoundary") + Foam::name(procI)
                   + word("to")
                   + Foam::name(curNeighbourProcessors[procPatchI]),
                     curProcessorPatchSizes[procPatchI],

src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C

@@ -650,8 +650,8 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
     // Dump first valid mesh without checking
     {
-        Pout<< "Dump mesh 0" << endl;
         const label fvmId = firstValidMesh();
+        Pout<< "Dump mesh " << fvmId << endl;
 
         cellOffset[fvmId] = 0;
 
@@ -795,13 +795,14 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
     for (label procI = 1; procI < meshes_.size(); procI++)
     {
-        Pout<< "Dump mesh " << procI << endl;
         // Grab cell offset from previous offset and mesh size
         cellOffset[procI] =
             cellOffset[procI - 1] + meshes_[procI - 1].nCells();
 
         if (meshes_.set(procI))
         {
+            Pout<< "Dump mesh " << procI << endl;
+
             const polyMesh& curMesh = meshes_[procI];
             const polyBoundaryMesh& procPatches = curMesh.boundaryMesh();