From 42a9078dad14ff0e7394eb9676c14852604496c2 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 20 Oct 2021 14:51:38 +0100
Subject: [PATCH] ENH: decomposePar: cyclicAMI in constraints. Fixes #2242

---
 .../decompositionConstraint.C                 | 93 ++++++++++++++++-
 .../decompositionConstraint.H                 |  9 ++
 .../preserveFaceZonesConstraint.C             | 99 +++++++++----------
 .../preservePatchesConstraint.C               | 76 +++++++-------
 4 files changed, 182 insertions(+), 95 deletions(-)

diff --git a/src/parallel/decompose/decompositionMethods/decompositionConstraints/decompositionConstraint/decompositionConstraint.C b/src/parallel/decompose/decompositionMethods/decompositionConstraints/decompositionConstraint/decompositionConstraint.C
index 7c24d3dec8..12f8a316c6 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionConstraints/decompositionConstraint/decompositionConstraint.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionConstraints/decompositionConstraint/decompositionConstraint.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2019 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,6 +27,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "decompositionConstraint.H"
+#include "syncTools.H"
+#include "cyclicAMIPolyPatch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -36,6 +38,95 @@ namespace Foam
     defineRunTimeSelectionTable(decompositionConstraint, dictionary);
 }
 
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::decompositionConstraint::getMinBoundaryValue
+(
+    const polyMesh& mesh,
+    const labelList& decomposition,
+    labelList& destProc
+) const
+{
+    destProc.setSize(mesh.nBoundaryFaces());
+    destProc = labelMax;
+
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    for (const polyPatch& pp : pbm)
+    {
+        const labelUList& faceCells = pp.faceCells();
+
+        forAll(faceCells, i)
+        {
+            label bFacei = pp.offset()+i;
+            destProc[bFacei] = decomposition[faceCells[i]];
+        }
+    }
+
+    // Take minimum of coupled faces (over all patches!)
+    syncTools::syncBoundaryFaceList(mesh, destProc, minEqOp<label>());
+
+    // Do cyclicAMI ourselves
+    for (const auto& pp : pbm)
+    {
+        const auto* ppp = isA<cyclicAMIPolyPatch>(pp);
+        if (ppp)
+        {
+            const auto& cycPp = *ppp;
+            const auto& nbrPp = cycPp.neighbPatch();
+            const labelList nbrDecomp(decomposition, nbrPp.faceCells());
+            labelList thisDecomp(decomposition, cycPp.faceCells());
+
+            if (cycPp.owner())
+            {
+                cycPp.AMI().interpolateToSource
+                (
+                    nbrDecomp,
+                    []
+                    (
+                        label& res,
+                        const label facei,
+                        const label& fld,
+                        const scalar& w
+                    )
+                    {
+                        res = min(res, fld);
+                    },
+                    thisDecomp,
+                    thisDecomp      // used in case of low-weight-corr
+                );
+            }
+            else
+            {
+                nbrPp.AMI().interpolateToTarget
+                (
+                    nbrDecomp,
+                    []
+                    (
+                        label& res,
+                        const label facei,
+                        const label& fld,
+                        const scalar& w
+                    )
+                    {
+                        res = min(res, fld);
+                    },
+                    thisDecomp,
+                    thisDecomp      // used in case of low-weight-corr
+                );
+            }
+
+            forAll(thisDecomp, i)
+            {
+                label& proc = destProc[cycPp.offset()+i];
+                proc = min(proc, thisDecomp[i]);
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::decompositionConstraint::decompositionConstraint
diff --git a/src/parallel/decompose/decompositionMethods/decompositionConstraints/decompositionConstraint/decompositionConstraint.H b/src/parallel/decompose/decompositionMethods/decompositionConstraints/decompositionConstraint/decompositionConstraint.H
index 98913f4c7c..70ca0c1145 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionConstraints/decompositionConstraint/decompositionConstraint.H
+++ b/src/parallel/decompose/decompositionMethods/decompositionConstraints/decompositionConstraint/decompositionConstraint.H
@@ -68,6 +68,15 @@ protected:
 
     // Protected Member Functions
 
+        //- Get minimum label across coupled boundary faces
+        void getMinBoundaryValue
+        (
+            const polyMesh& mesh,
+            const labelList& decomposition,
+            labelList& destProc
+        ) const;
+
+
         //- No copy construct
         decompositionConstraint(const decompositionConstraint&) = delete;
 
diff --git a/src/parallel/decompose/decompositionMethods/decompositionConstraints/preserveFaceZones/preserveFaceZonesConstraint.C b/src/parallel/decompose/decompositionMethods/decompositionConstraints/preserveFaceZones/preserveFaceZonesConstraint.C
index 3f0a052778..5952e305a2 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionConstraints/preserveFaceZones/preserveFaceZonesConstraint.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionConstraints/preserveFaceZones/preserveFaceZonesConstraint.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018,2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -140,75 +140,64 @@ void Foam::decompositionConstraints::preserveFaceZones::apply
     // If the decomposition has not enforced the constraint do it over
     // here.
 
+    label nChanged;
 
-    // Synchronise decomposition on boundary
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
-
-    labelList destProc(mesh.nBoundaryFaces(), labelMax);
-
-    for (const polyPatch& pp : pbm)
+    do
     {
-        label bFacei = pp.start() - mesh.nInternalFaces();
+        // Extract min coupled boundary data
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        const labelUList& faceCells = pp.faceCells();
+        labelList destProc;
+        getMinBoundaryValue(mesh, decomposition, destProc);
 
-        for (const label celli : faceCells)
+
+        // Override if differing
+        // ~~~~~~~~~~~~~~~~~~~~~
+
+        const faceZoneMesh& fZones = mesh.faceZones();
+
+        const labelList zoneIDs(zones_.matching(fZones.names()));
+
+        nChanged = 0;
+        for (const label zonei : zoneIDs)
         {
-            destProc[bFacei] = decomposition[celli];
+            const faceZone& fz = fZones[zonei];
 
-            ++bFacei;
-        }
-    }
-
-    syncTools::syncBoundaryFaceList(mesh, destProc, minEqOp<label>());
-
-
-    // Override if differing
-    // ~~~~~~~~~~~~~~~~~~~~~
-
-    const faceZoneMesh& fZones = mesh.faceZones();
-
-    const labelList zoneIDs(zones_.matching(fZones.names()));
-
-    label nChanged = 0;
-
-    for (const label zonei : zoneIDs)
-    {
-        const faceZone& fz = fZones[zonei];
-
-        for (const label facei : fz)
-        {
-            const label own = mesh.faceOwner()[facei];
-
-            if (mesh.isInternalFace(facei))
+            for (const label facei : fz)
             {
-                const label nei = mesh.faceNeighbour()[facei];
-                if (decomposition[own] != decomposition[nei])
+                const label own = mesh.faceOwner()[facei];
+
+                if (mesh.isInternalFace(facei))
                 {
-                    decomposition[nei] = decomposition[own];
-                    ++nChanged;
+                    const label nei = mesh.faceNeighbour()[facei];
+                    if (decomposition[nei] < decomposition[own])
+                    {
+                        decomposition[own] = decomposition[nei];
+                        ++nChanged;
+                    }
                 }
-            }
-            else
-            {
-                const label bFaceI = facei-mesh.nInternalFaces();
-                if (decomposition[own] != destProc[bFaceI])
+                else
                 {
-                    decomposition[own] = destProc[bFaceI];
-                    ++nChanged;
+                    const label bFaceI = facei-mesh.nInternalFaces();
+                    if (destProc[bFaceI] < decomposition[own])
+                    {
+                        decomposition[own] = destProc[bFaceI];
+                        ++nChanged;
+                    }
                 }
             }
         }
-    }
 
-    if (decompositionConstraint::debug & 2)
-    {
         reduce(nChanged, sumOp<label>());
-        Info<< type() << " : changed decomposition on " << nChanged
-            << " cells" << endl;
-    }
+
+        if (decompositionConstraint::debug & 2)
+        {
+            reduce(nChanged, sumOp<label>());
+            Info<< type() << " : changed decomposition on " << nChanged
+                << " cells" << endl;
+        }
+
+    } while (nChanged > 0);
 }
 
 
diff --git a/src/parallel/decompose/decompositionMethods/decompositionConstraints/preservePatches/preservePatchesConstraint.C b/src/parallel/decompose/decompositionMethods/decompositionConstraints/preservePatches/preservePatchesConstraint.C
index 2c0f6b1382..10f65341aa 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionConstraints/preservePatches/preservePatchesConstraint.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionConstraints/preservePatches/preservePatchesConstraint.C
@@ -63,7 +63,8 @@ Foam::decompositionConstraints::preservePatches::preservePatches
         Info<< type()
             << " : adding constraints to keep owner and (coupled) neighbour"
             << " of faces in patches " << patches_
-            << " on same processor. This only makes sense for cyclics." << endl;
+            << " on same processor. This only makes sense for cyclics"
+            << " and cyclicAMI." << endl;
     }
 }
 
@@ -81,7 +82,8 @@ Foam::decompositionConstraints::preservePatches::preservePatches
         Info<< type()
             << " : adding constraints to keep owner and (coupled) neighbour"
             << " of faces in patches " << patches_
-            << " on same processor. This only makes sense for cyclics." << endl;
+            << " on same processor. This only makes sense for cyclics"
+            << " and cyclicAMI." << endl;
     }
 }
 
@@ -143,58 +145,54 @@ void Foam::decompositionConstraints::preservePatches::apply
 {
     // If the decomposition has not enforced the constraint, do it over here.
 
-    // Synchronise decomposition on patchIDs
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
-    labelList destProc(mesh.nBoundaryFaces(), labelMax);
-
-    for (const polyPatch& pp : pbm)
-    {
-        const labelUList& faceCells = pp.faceCells();
-
-        forAll(faceCells, i)
-        {
-            label bFaceI = pp.start()+i-mesh.nInternalFaces();
-            destProc[bFaceI] = decomposition[faceCells[i]];
-        }
-    }
-
-    syncTools::syncBoundaryFaceList(mesh, destProc, minEqOp<label>());
-
-
-    // Override if differing
-    // ~~~~~~~~~~~~~~~~~~~~~
-
     const labelList patchIDs(pbm.patchSet(patches_).sortedToc());
 
-    label nChanged = 0;
+    // Synchronise decomposition on patchIDs
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    for (const label patchi : patchIDs)
+    label nChanged;
+
+    do
     {
-        const polyPatch& pp = pbm[patchi];
+        // Extract min coupled boundary data
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        const labelUList& faceCells = pp.faceCells();
+        labelList destProc;
+        getMinBoundaryValue(mesh, decomposition, destProc);
 
-        forAll(faceCells, i)
+
+        // Override (patchIDs only) if differing
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        nChanged = 0;
+        for (const label patchi : patchIDs)
         {
-            const label bFaceI = pp.start()+i-mesh.nInternalFaces();
+            const polyPatch& pp = pbm[patchi];
 
-            if (decomposition[faceCells[i]] != destProc[bFaceI])
+            const labelUList& faceCells = pp.faceCells();
+
+            forAll(faceCells, i)
             {
-                decomposition[faceCells[i]] = destProc[bFaceI];
-                ++nChanged;
+                const label bFacei = pp.offset()+i;
+                if (destProc[bFacei] < decomposition[faceCells[i]])
+                {
+                    decomposition[faceCells[i]] = destProc[bFacei];
+                    ++nChanged;
+                }
             }
         }
-    }
 
-    if (decompositionConstraint::debug & 2)
-    {
         reduce(nChanged, sumOp<label>());
-        Info<< type() << " : changed decomposition on " << nChanged
-            << " cells" << endl;
-    }
+
+        if (decompositionConstraint::debug & 2)
+        {
+            Info<< type() << " : changed decomposition on " << nChanged
+                << " cells" << endl;
+        }
+
+    } while (nChanged > 0);
 }