ENH: Support AMI for multi-world operation. Fixes #2099

Multi-world operation now supports AMI:

    // What to sample:
    sampleMode      nearestPatchFaceAMI;

src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.C

@@ -67,6 +67,7 @@ void Foam::mappedPatchFieldBase<Type>::storeField
     for (label domain = 0; domain < nProcs; domain++)
     {
         const labelList& map = procToMap[domain];
+
         if (map.size())
         {
             const Field<T> subFld(fld, map);
@@ -100,7 +101,7 @@ void Foam::mappedPatchFieldBase<Type>::storeField
 
 template<class Type>
 template<class T>
-void Foam::mappedPatchFieldBase<Type>::retrieveField
+bool Foam::mappedPatchFieldBase<Type>::retrieveField
 (
     const bool allowUnset,
     const objectRegistry& obr,
@@ -112,9 +113,10 @@ void Foam::mappedPatchFieldBase<Type>::retrieveField
 ) const
 {
     // Store my data onto database
-    //const label myRank = Pstream::myProcNo(0);  // comm_
     const label nProcs = Pstream::nProcs(0);    // comm_
 
+    bool ok = true;
+
     for (label domain = 0; domain < nProcs; domain++)
     {
         const labelList& map = procToMap[domain];
@@ -128,39 +130,68 @@ void Foam::mappedPatchFieldBase<Type>::retrieveField
               / patch
             );
 
-            //const IOField<T>& subFld = subObr.lookupObject<IOField<T>>
-            //(
-            //    fieldName
-            //);
             const IOField<T>* subFldPtr = subObr.getObjectPtr<IOField<T>>
             (
                 fieldName
             );
             if (subFldPtr)
             {
-                UIndirectList<T>(fld, map) = *subFldPtr;
-
-                if (fvPatchField<Type>::debug)
+                if (subFldPtr->size() != map.size())
                 {
-                    Pout<< "*** RETRIEVED :"
-                        << " field:" << fieldName
-                        << " values:" << flatOutput(fld)
-                        << " from:" << subObr.objectPath() << endl;
+                    // This is the dummy value inserted at start-up since the
+                    // map is always non-zero size (checked above)
+                    //Pout<< "*** RETRIEVED DUMMY :"
+                    //    << " field:" << fieldName
+                    //    << " subFldPtr:" << subFldPtr->size()
+                    //    << " map:" << map.size() << endl;
+
+                    ok = false;
+                }
+                else
+                {
+                    UIndirectList<T>(fld, map) = *subFldPtr;
+
+                    if (fvPatchField<Type>::debug)
+                    {
+                        Pout<< "*** RETRIEVED :"
+                            << " field:" << fieldName
+                            << " values:" << flatOutput(fld)
+                            << " from:" << subObr.objectPath() << endl;
+                    }
                 }
             }
             else if (allowUnset)
             {
-                WarningInFunction << "Not found"
-                    << " field:" << fieldName
-                    << " in:" << subObr.objectPath() << endl;
+                if (fvPatchField<Type>::debug)
+                {
+                    WarningInFunction << "Not found"
+                        << " field:" << fieldName
+                        << " in:" << subObr.objectPath() << endl;
+                }
+
+                // Store dummy value so the database has something on it.
+                // Note that size 0 should never occur naturally so we can
+                // detect it if necessary.
+                const Field<T> dummyFld(0);
+
+                mappedPatchBase::storeField
+                (
+                    const_cast<objectRegistry&>(subObr),
+                    fieldName,
+                    dummyFld
+                );
+
+                ok = false;
             }
             else
             {
                 // Not found. Make it fail
                 (void)subObr.lookupObject<IOField<T>>(fieldName);
+                ok = false;
             }
         }
     }
+    return ok;
 }
 
 
@@ -171,18 +202,17 @@ void Foam::mappedPatchFieldBase<Type>::initRetrieveField
     const objectRegistry& obr,
     const word& region,
     const word& patch,
-    const mapDistribute& map,
+    const labelListList& map,
     const word& fieldName,
     const Field<T>& fld
 ) const
 {
     // Store my data onto database
-    //const label myRank = Pstream::myProcNo(0);  // comm_
     const label nProcs = Pstream::nProcs(0);    // comm_
 
     for (label domain = 0; domain < nProcs; domain++)
     {
-        const labelList& constructMap = map.constructMap()[domain];
+        const labelList& constructMap = map[domain];
         if (constructMap.size())
         {
             const objectRegistry& subObr = mappedPatchBase::subRegistry
@@ -216,6 +246,71 @@ void Foam::mappedPatchFieldBase<Type>::initRetrieveField
 }
 
 
+template<class Type>
+template<class T>
+bool Foam::mappedPatchFieldBase<Type>::storeAndRetrieveField
+(
+    const word& fieldName,
+    const labelListList& subMap,
+    const label constructSize,
+    const labelListList& constructMap,
+    const labelListList& address,
+    const scalarListList& weights,
+    Field<T>& fld
+) const
+{
+    storeField
+    (
+        patchField_.internalField().time(),
+        patchField_.patch().boundaryMesh().mesh().name(),
+        patchField_.patch().name(),
+        subMap,
+        fieldName,
+        fld
+    );
+
+    Field<T> work(constructSize);
+    const bool ok = retrieveField
+    (
+        true,                           // allow unset
+        patchField_.internalField().time(),
+        mapper_.sampleRegion(),
+        mapper_.samplePatch(),
+        constructMap,
+        fieldName,
+        work
+    );
+
+    if (ok)
+    {
+        // Do interpolation
+
+        fld.setSize(address.size());
+        fld = Zero;
+
+        const plusEqOp<T> cop;
+        const multiplyWeightedOp<T, plusEqOp<T>> mop(cop);
+
+        forAll(address, facei)
+        {
+            const labelList& slots = address[facei];
+            const scalarList& w = weights[facei];
+
+            forAll(slots, i)
+            {
+                mop(fld[facei], facei, work[slots[i]], w[i]);
+            }
+        }
+    }
+    else
+    {
+        // Leave fld intact
+    }
+
+    return ok;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
@@ -263,16 +358,24 @@ Foam::mappedPatchFieldBase<Type>::mappedPatchFieldBase
     if
     (
         mapper_.sampleDatabase()
-     && mapper_.mode() != mappedPatchBase::NEARESTPATCHFACE
+     && (
+            mapper_.mode() != mappedPatchBase::NEARESTPATCHFACE
+         && mapper_.mode() != mappedPatchBase::NEARESTPATCHFACEAMI
+        )
     )
     {
         FatalErrorInFunction
             << "Mapping using the database only supported for "
-            << "sampleMode "
+            << "sampleModes "
             <<  mappedPatchBase::sampleModeNames_
                 [
                     mappedPatchBase::NEARESTPATCHFACE
                 ]
+            << " and "
+            <<  mappedPatchBase::sampleModeNames_
+                [
+                    mappedPatchBase::NEARESTPATCHFACEAMI
+                ]
             << exit(FatalError);
     }
 
@@ -297,22 +400,29 @@ Foam::mappedPatchFieldBase<Type>::mappedPatchFieldBase
 :
     mappedPatchFieldBase<Type>::mappedPatchFieldBase(mapper, patchField, dict)
 {
-    if
-    (
-        mapper_.mode() == mappedPatchBase::NEARESTPATCHFACE
-     && mapper_.sampleDatabase()
-    )
+    if (mapper_.sampleDatabase())
     {
-        // Store my data on receive buffers so we have some initial data
-        initRetrieveField
-        (
-            patchField_.internalField().time(),
-            patchField_.patch().boundaryMesh().mesh().name(),
-            patchField_.patch().name(),
-            mapper_.map(),
-            patchField_.internalField().name(),
-            patchField_
-        );
+        if (mapper_.mode() == mappedPatchBase::NEARESTPATCHFACE)
+        {
+            // Store my data on receive buffers so we have some initial data
+            initRetrieveField
+            (
+                patchField_.internalField().time(),
+                //patchField_.patch().boundaryMesh().mesh().name(),
+                mapper_.sampleRegion(),
+                //patchField_.patch().name(),
+                mapper_.samplePatch(),
+                mapper_.map().constructMap(),
+                patchField_.internalField().name(),
+                patchField_
+            );
+        }
+        else if (mapper_.mode() == mappedPatchBase::NEARESTPATCHFACEAMI)
+        {
+            // Depend on fall-back (sorting dummy field) in retrieveField
+            // since it would be too hard to determine the field that gives
+            // the wanted result after interpolation
+        }
     }
 }
 
@@ -410,37 +520,80 @@ template<class T>
 void Foam::mappedPatchFieldBase<Type>::distribute
 (
     const word& fieldName,
-    Field<T>& newValues
+    Field<T>& fld
 ) const
 {
     if (mapper_.sampleDatabase())
     {
-        // Store my data on send buffers
-        storeField
-        (
-            patchField_.internalField().time(),
-            patchField_.patch().boundaryMesh().mesh().name(),
-            patchField_.patch().name(),
-            mapper_.map().subMap(),
-            fieldName,
-            newValues
-        );
-        // Construct my data from receive buffers
-        newValues.setSize(mapper_.map().constructSize());
-        retrieveField
-        (
-            true,                           // allow unset
-            patchField_.internalField().time(),
-            mapper_.sampleRegion(),
-            mapper_.samplePatch(),
-            mapper_.map().constructMap(),
-            fieldName,
-            newValues
-        );
+        if (mapper_.mode() != mappedPatchBase::NEARESTPATCHFACEAMI)
+        {
+            // Store my data on send buffers
+            storeField
+            (
+                patchField_.internalField().time(),
+                patchField_.patch().boundaryMesh().mesh().name(),
+                patchField_.patch().name(),
+                mapper_.map().subMap(),
+                fieldName,
+                fld
+            );
+            // Construct my data from receive buffers
+            fld.setSize(mapper_.map().constructSize());
+            retrieveField
+            (
+                true,                           // allow unset
+                patchField_.internalField().time(),
+                mapper_.sampleRegion(),
+                mapper_.samplePatch(),
+                mapper_.map().constructMap(),
+                fieldName,
+                fld
+            );
+        }
+        else
+        {
+            const AMIPatchToPatchInterpolation& AMI = mapper_.AMI();
+
+            // The AMI does an interpolateToSource/ToTarget. This is a
+            // mapDistribute (so using subMap/constructMap) and then a
+            // weighted sum. We'll store the sent data as before and
+            // do the weighted summation after the retrieveField
+
+            if (mapper_.masterWorld())
+            {
+                // See AMIInterpolation::interpolateToSource. Use tgtMap,
+                // srcAddress, srcWeights
+                storeAndRetrieveField
+                (
+                    fieldName,
+                    AMI.srcMap().subMap(),
+                    AMI.tgtMap().constructSize(),
+                    AMI.tgtMap().constructMap(),
+                    AMI.srcAddress(),
+                    AMI.srcWeights(),
+                    fld
+                );
+            }
+            else
+            {
+                // See AMIInterpolation::interpolateToTarget.
+                // Use srcMap, tgtAddress, tgtWeights
+                storeAndRetrieveField
+                (
+                    fieldName,
+                    AMI.tgtMap().subMap(),
+                    AMI.srcMap().constructSize(),
+                    AMI.srcMap().constructMap(),
+                    AMI.tgtAddress(),
+                    AMI.tgtWeights(),
+                    fld
+                );
+            }
+        }
     }
     else
     {
-        mapper_.distribute(newValues);
+        mapper_.distribute(fld);
     }
 }
 
@@ -825,15 +978,18 @@ void Foam::mappedPatchFieldBase<Type>::initRetrieveField
     {
         // Store my data on receive buffers (reverse of storeField;
         // i.e. retrieveField will obtain patchField)
-        initRetrieveField
-        (
-            patchField_.internalField().time(),
-            mapper_.sampleRegion(),
-            mapper_.samplePatch(),
-            mapper_.map(),
-            fieldName,
-            fld
-        );
+        if (mapper_.mode() == mappedPatchBase::NEARESTPATCHFACE)
+        {
+            initRetrieveField
+            (
+                patchField_.internalField().time(),
+                mapper_.sampleRegion(),
+                mapper_.samplePatch(),
+                mapper_.map().constructMap(),
+                fieldName,
+                fld
+            );
+        }
     }
 }
 

src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -116,11 +117,26 @@ protected:
             const objectRegistry& obr,
             const word& region,
             const word& patch,
-            const mapDistribute& map,
+            const labelListList& map,
             const word& fieldName,
             const Field<T>& fld
         ) const;
 
+        //- Helper : storeField and retrieveField and interpolate. Leaves fld
+        //  unchanged (and returns false) if new values cannot be retrieved.
+        //  Returns true otherwise.
+        template<class T>
+        bool storeAndRetrieveField
+        (
+            const word& fieldName,
+            const labelListList& subMap,
+            const label constructSize,
+            const labelListList& constructMap,
+            const labelListList& address,
+            const scalarListList& weights,
+            Field<T>& fld
+        ) const;
+
 
 public:
 
@@ -244,9 +260,10 @@ public:
                 const Field<T>& fld
             ) const;
 
-            //- Construct field from registered elements
+            //- Construct field from registered elements. Return true if
+            //  successful
             template<class T>
-            void retrieveField
+            bool retrieveField
             (
                 const bool allowUnset,
                 const objectRegistry& obr,

src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C

@@ -1105,14 +1105,14 @@ void Foam::mappedPatchBase::calcAMI() const
     }
 
     // Check if running locally
-    if (sampleWorld_.empty())
+    if (sampleWorld_.empty() || sameWorld())
     {
         const polyPatch& nbr = samplePolyPatch();
 
         // Transform neighbour patch to local system
-        pointField nbrPoints(samplePoints(nbr.localPoints()));
+        const pointField nbrPoints(samplePoints(nbr.localPoints()));
 
-        primitivePatch nbrPatch0
+        const primitivePatch nbrPatch0
         (
             SubList<face>
             (
@@ -1135,43 +1135,56 @@ void Foam::mappedPatchBase::calcAMI() const
             meshTools::writeOBJ(osO, patch_.localFaces(), patch_.localPoints());
         }
 
-        // Construct/apply AMI interpolation to determine addressing and weights
+        // Construct/apply AMI interpolation to determine addressing and
+        // weights. Make sure to use optional inter-world communicator.
+
+        const label oldWorldComm = Pstream::worldComm;
+        Pstream::worldComm = comm_;
+
+        const label oldComm(Pstream::warnComm);
+        Pstream::warnComm = UPstream::worldComm;
+
         AMIPtr_->calculate(patch_, nbrPatch0, surfPtr());
+
+        Pstream::warnComm = oldComm;
+        Pstream::worldComm = oldWorldComm;
     }
     else
     {
-        faceList nbrFaces;
-        pointField nbrPoints;
-        if (sampleWorld() == UPstream::myWorld())
-        {
-            const polyPatch& nbr = samplePolyPatch();
-            nbrFaces = nbr.localFaces();
-            nbrPoints = samplePoints(nbr.localPoints());
-        }
-        else
-        {
-            // Leave empty
-        }
-
-        primitivePatch nbrPatch0
+        faceList dummyFaces;
+        pointField dummyPoints;
+        const primitivePatch dummyPatch
         (
             SubList<face>
             (
-                nbrFaces,
-                nbrFaces.size()
+                dummyFaces
             ),
-            nbrPoints
+            dummyPoints
         );
 
-        // Change to use ALL processors communicator
+        // Change to use inter-world communicator
         const label oldWorldComm = Pstream::worldComm;
         Pstream::worldComm = comm_;
 
         const label oldComm(Pstream::warnComm);
         Pstream::warnComm = UPstream::worldComm;
 
-        // Construct/apply AMI interpolation to determine addressing and weights
-        AMIPtr_->calculate(patch_, nbrPatch0, surfPtr());
+        if (masterWorld())
+        {
+            // Construct/apply AMI interpolation to determine addressing
+            // and weights. Have patch_ for src faces, 0 faces for the
+            // target side
+            AMIPtr_->calculate(patch_, dummyPatch, surfPtr());
+        }
+        else
+        {
+            // Construct/apply AMI interpolation to determine addressing
+            // and weights. Have 0 faces for src side, patch_ for the tgt
+            // side
+            AMIPtr_->calculate(dummyPatch, patch_, surfPtr());
+        }
+        // Now the AMI addressing/weights will be from src side (on masterWorld
+        // processors) to tgt side (on other processors)
 
         Pstream::warnComm = oldComm;
         Pstream::worldComm = oldWorldComm;

src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H

@@ -492,9 +492,12 @@ public:
             //- Communicator
             inline label comm() const;
 
-            //- Is world the local world
+            //- Is sample world the local world?
             inline bool sameWorld() const;
 
+            //- Is my world ordered before the sampleWorld?
+            inline bool masterWorld() const;
+
             //- Cached sampleRegion != mesh.name()
             inline bool sameRegion() const;
 

src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBaseI.H

@@ -158,6 +158,23 @@ inline bool Foam::mappedPatchBase::sameWorld() const
 }
 
 
+inline bool Foam::mappedPatchBase::masterWorld() const
+{
+    if (sameWorld())
+    {
+        return true;
+    }
+    else
+    {
+        // Use ordering in allWorlds
+        const label myWorld = UPstream::myWorldID();
+        const label mySampleWorld =
+            UPstream::allWorlds().find(sampleWorld_);
+        return myWorld < mySampleWorld;
+    }
+}
+
+
 inline bool Foam::mappedPatchBase::sameRegion() const
 {