cyclicACMIGAMGInterface.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2013-2016 OpenFOAM Foundation
    Copyright (C) 2019-2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "AMIInterpolation.H"
#include "cyclicACMIGAMGInterface.H"
#include "addToRunTimeSelectionTable.H"
#include "Map.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(cyclicACMIGAMGInterface, 0);
    addToRunTimeSelectionTable
    (
        GAMGInterface,
        cyclicACMIGAMGInterface,
        lduInterface
    );
    addToRunTimeSelectionTable
    (
        GAMGInterface,
        cyclicACMIGAMGInterface,
        Istream
    );
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::cyclicACMIGAMGInterface::cyclicACMIGAMGInterface
(
    const label index,
    const lduInterfacePtrsList& coarseInterfaces,
    const lduInterface& fineInterface,
    const labelField& localRestrictAddressing,
    const labelField& neighbourRestrictAddressing,
    const label fineLevelIndex,
    const label coarseComm
)
:
    GAMGInterface
    (
        index,
        coarseInterfaces
    ),
    neighbPatchID_
    (
        refCast<const cyclicACMILduInterface>(fineInterface).neighbPatchID()
    ),
    owner_
    (
        refCast<const cyclicACMILduInterface>(fineInterface).owner()
    ),
    forwardT_
    (
        refCast<const cyclicACMILduInterface>(fineInterface).forwardT()
    ),
    reverseT_
    (
        refCast<const cyclicACMILduInterface>(fineInterface).reverseT()
    ),
    myProcNo_(-1)
{
    const auto& fineCyclicACMIInterface =
        refCast<const cyclicACMILduInterface>(fineInterface);

    // Construct face agglomeration from cell agglomeration
    {
        // From coarse face to cell
        DynamicList<label> dynFaceCells(localRestrictAddressing.size());

        // From face to coarse face
        DynamicList<label> dynFaceRestrictAddressing
        (
            localRestrictAddressing.size()
        );

        Map<label> masterToCoarseFace(localRestrictAddressing.size());

        for (const label curMaster : localRestrictAddressing)
        {
            const auto iter = masterToCoarseFace.cfind(curMaster);

            if (iter.good())
            {
                // Already have coarse face
                dynFaceRestrictAddressing.append(iter.val());
            }
            else
            {
                // New coarse face
                const label coarseI = dynFaceCells.size();
                dynFaceRestrictAddressing.append(coarseI);
                dynFaceCells.append(curMaster);
                masterToCoarseFace.insert(curMaster, coarseI);
            }
        }

        faceCells_.transfer(dynFaceCells);
        faceRestrictAddressing_.transfer(dynFaceRestrictAddressing);
    }


    // On the owner side construct the AMI

    if (fineCyclicACMIInterface.owner())
    {
        // Construct the neighbour side agglomeration (as the neighbour would
        // do it so it the exact loop above using neighbourRestrictAddressing
        // instead of localRestrictAddressing)

        labelList nbrFaceRestrictAddressing;
        {
            // From face to coarse face
            DynamicList<label> dynNbrFaceRestrictAddressing
            (
                neighbourRestrictAddressing.size()
            );

            Map<label> masterToCoarseFace(neighbourRestrictAddressing.size());

            for (const label curMaster : neighbourRestrictAddressing)
            {
                const auto iter = masterToCoarseFace.cfind(curMaster);

                if (iter.good())
                {
                    // Already have coarse face
                    dynNbrFaceRestrictAddressing.append(iter.val());
                }
                else
                {
                    // New coarse face
                    const label coarseI = masterToCoarseFace.size();
                    dynNbrFaceRestrictAddressing.append(coarseI);
                    masterToCoarseFace.insert(curMaster, coarseI);
                }
            }

            nbrFaceRestrictAddressing.transfer(dynNbrFaceRestrictAddressing);
        }


        amiPtr_.reset
        (
            new AMIPatchToPatchInterpolation
            (
                fineCyclicACMIInterface.AMI(),
                faceRestrictAddressing_,
                nbrFaceRestrictAddressing
            )
        );
    }
}


Foam::cyclicACMIGAMGInterface::cyclicACMIGAMGInterface
(
    const label index,
    const lduInterfacePtrsList& coarseInterfaces,
    Istream& is
)
:
    GAMGInterface(index, coarseInterfaces, is),
    neighbPatchID_(readLabel(is)),
    owner_(readBool(is)),
    forwardT_(is),
    reverseT_(is),
    myProcNo_(-1)
{
    const bool hasAMI(readBool(is));

    if (hasAMI)
    {
        amiPtr_.reset(new AMIPatchToPatchInterpolation(is));

        // Store originating ranks locally - used when processor agglomerating
        // onto a processor that wasn't in the communicator originally (since
        // it had no faces)
        const label comm = AMI().comm();

        if (comm != -1)
        {
            is  >> myProcNo_;
        }
    }
}


Foam::cyclicACMIGAMGInterface::cyclicACMIGAMGInterface
(
    const label index,
    const lduInterfacePtrsList& coarseInterfaces,
    const lduInterface& fineInterface,
    const labelList& interfaceMap,
    const labelUList& faceCells,
    const labelUList& faceRestrictAddresssing,
    const labelUList& faceOffsets,
    const lduInterfacePtrsList& allInterfaces,
    const label coarseComm,
    const label myProcNo,
    const labelList& procAgglomMap
)
:
    GAMGInterface
    (
        index,
        coarseInterfaces,
        faceCells,
        faceRestrictAddresssing
    ),
    neighbPatchID_
    (
        interfaceMap.find
        (
            refCast
            <
                const cyclicACMILduInterface
            >(fineInterface).neighbPatchID()
        )
    ),
    owner_
    (
        refCast<const cyclicACMILduInterface>(fineInterface).owner()
    ),
    forwardT_
    (
        refCast<const cyclicACMILduInterface>(fineInterface).forwardT()
    ),
    reverseT_
    (
       refCast<const cyclicACMILduInterface>(fineInterface).reverseT()
    ),
    myProcNo_(-1)
{
    if (!owner_)
    {
        return;
    }


    // Get stats, sizes from the input interfaces. For the global settings
    // the problem is that the
    // local processor might not have any valid interfaces so here just
    // collect and do a global reduction afterwards.

    // Structure to pack all. First element is used to decide who has the
    // valid AMI.
    typedef
    Tuple2
    <
        label,
        Tuple2
        <
            Tuple2
            <
                FixedList<bool, 4>,
                scalar
            >,
            label
        >
    > AMIType;

    AMIType globalInfo;
    FixedList<bool, 4>& bools = globalInfo.second().first().first();

    // Define aliases to make our life easier
    label& firstValidAMI = globalInfo.first();
    bool& requireMatch = bools[0];
    bool& reverseTarget = bools[1];
    bool& srcHasFlip = bools[2];
    bool& tgtHasFlip = bools[3];
    scalar& lowWeightCorrection = globalInfo.second().first().second();
    label& singlePatchProc = globalInfo.second().second();

    // Initialise all global variables
    firstValidAMI = labelMax;
    requireMatch = false;
    reverseTarget = false;
    srcHasFlip = false;
    tgtHasFlip = false;
    lowWeightCorrection = -1;
    singlePatchProc = -1;

    // Initialise all local variables
    bool hasSrcMagSf = false;
    bool hasSrcCentroids = false;
    bool hasTgtMagSf = false;
    label nSrc = 0;
    label nTgt = 0;

    forAll(allInterfaces, inti)
    {
        if (allInterfaces.set(inti))
        {
            const auto& intf =
                refCast<const cyclicACMIGAMGInterface>(allInterfaces[inti]);

            if (!intf.amiPtr_)
            {
                continue;
            }

            if (firstValidAMI == labelMax)
            {
                firstValidAMI = inti;
            }

            const auto& AMI = intf.AMI();

            if (AMI.distributed() && AMI.comm() != -1)
            {
                singlePatchProc = -1;
                srcHasFlip =
                    srcHasFlip || AMI.srcMap().constructHasFlip();
                tgtHasFlip =
                    tgtHasFlip || AMI.tgtMap().constructHasFlip();
            }
            requireMatch = AMI.requireMatch();
            reverseTarget = AMI.reverseTarget();
            lowWeightCorrection = AMI.lowWeightCorrection();

            nSrc += AMI.srcAddress().size();
            nTgt += AMI.tgtAddress().size();

            if (AMI.srcMagSf().size())
            {
                hasSrcMagSf = true;
                if (AMI.srcMagSf().size() != AMI.srcAddress().size())
                {
                    FatalErrorInFunction
                        << "srcMagSf size:" << AMI.srcMagSf().size()
                        << "srcAddress size:" << AMI.srcAddress().size()
                        << exit(FatalError);
                }
            }
            if (AMI.srcCentroids().size())
            {
                hasSrcCentroids = true;
                if (AMI.srcCentroids().size() != AMI.srcAddress().size())
                {
                    FatalErrorInFunction
                        << "srcCentroids size:" << AMI.srcCentroids().size()
                        << "srcAddress size:" << AMI.srcAddress().size()
                        << exit(FatalError);
                }
            }
            if (AMI.tgtMagSf().size())
            {
                hasTgtMagSf = true;
                if (AMI.tgtMagSf().size() != AMI.tgtAddress().size())
                {
                    FatalErrorInFunction
                        << "tgtMagSf size:" << AMI.tgtMagSf().size()
                        << "tgtAddress size:" << AMI.tgtAddress().size()
                        << exit(FatalError);
                }
            }
        }
    }


    // Reduce global information in case one of the coarse ranks does not
    // have an input AMI to get data from. Could use minFirstEqOp from Tuple2
    // instead ...
    Pstream::combineReduce
    (
        globalInfo,
        [](AMIType& x, const AMIType& y)
        {
            if (y.first() < x.first())
            {
                x = y;
            }
        },
        Pstream::msgType(),
        coarseComm
    );

    DebugPout
        << "Input amis :"
        << " singlePatchProc:" << singlePatchProc
        << " srcHasFlip:" << srcHasFlip
        << " tgtHasFlip:" << tgtHasFlip
        << " requireMatch:" << requireMatch
        << " reverseTarget:" << reverseTarget
        << " lowWeightCorrection:" << lowWeightCorrection
        << " hasSrcMagSf:" << hasSrcMagSf
        << " hasSrcCentroids:" << hasSrcCentroids
        << " hasTgtMagSf:" << hasTgtMagSf
        << " nSrc:" << nSrc
        << " nTgt:" << nTgt
        << endl;


    labelListList srcAddress;
    scalarListList srcWeights;
    scalarList srcMagSf;
    // Needed?
    pointListList srcCentroids;

    labelListList tgtAddress;
    scalarListList tgtWeights;
    scalarList tgtMagSf;


    // Map to send src side data to tgt side
    autoPtr<mapDistribute> srcToTgtMap;

    // Map to send tgt side data to src side
    autoPtr<mapDistribute> tgtToSrcMap;

    if (singlePatchProc == -1)
    {
        // Find ranks that agglomerate together
        const label myAgglom = UPstream::myProcNo(coarseComm);

        // Per input map either -1 or the index in the maps that is local
        // data.
        labelList localRanks(allInterfaces.size(), -1);
        // From rank in coarse communicator back to rank in original (fine)
        // communicator.
        labelListList newToOldRanks;
        {
            // Pass 1: count number of valid maps
            label nOldRanks = 0;
            forAll(allInterfaces, inti)
            {
                if (allInterfaces.set(inti))
                {
                    const auto& intf = refCast<const cyclicACMIGAMGInterface>
                    (
                        allInterfaces[inti]
                    );

                    if (!intf.amiPtr_ || intf.AMI().comm() == -1)
                    {
                        continue;
                    }
                    nOldRanks++;
                }
            }

            // Pass 2: collect
            DynamicList<label> oldRanks(nOldRanks);
            forAll(allInterfaces, inti)
            {
                if (allInterfaces.set(inti))
                {
                    const auto& intf = refCast<const cyclicACMIGAMGInterface>
                    (
                        allInterfaces[inti]
                    );

                    if (!intf.amiPtr_ || intf.AMI().comm() == -1)
                    {
                        continue;
                    }

                    label fineRank = -1;
                    if (intf.myProcNo() == -1)
                    {
                        // The interface was already local so got never
                        // sent across so myProcNo_ is never set ...
                        fineRank = UPstream::myProcNo(intf.AMI().comm());
                    }
                    else
                    {
                        fineRank = intf.myProcNo();
                    }

                    oldRanks.append(fineRank);
                    localRanks[inti] = fineRank;
                }
            }

            // Pull individual parts together - this is the only communication
            // needed.
            newToOldRanks = Pstream::listGatherValues
            (
                labelList(std::move(oldRanks)),
                coarseComm
            );
            Pstream::broadcast(newToOldRanks, coarseComm);
        }


        // Create combined maps
        UPtrList<const mapDistribute> srcMaps(allInterfaces.size());
        UPtrList<const mapDistribute> tgtMaps(allInterfaces.size());
        forAll(allInterfaces, inti)
        {
            if (allInterfaces.set(inti))
            {
                const auto& intf = refCast<const cyclicACMIGAMGInterface>
                (
                    allInterfaces[inti]
                );

                if (!intf.amiPtr_)
                {
                    // Should not be in allInterfaces?
                    continue;
                }

                const auto& AMI = intf.AMI();

                if (AMI.comm() != -1)
                {
                    srcMaps.set(inti, &AMI.srcMap());
                    tgtMaps.set(inti, &AMI.tgtMap());
                }
            }
        }


        // Offsets for slots into results of srcToTgtMap
        labelList srcStartOfLocal;
        List<Map<label>> srcCompactMaps;

        srcToTgtMap.reset
        (
            new mapDistribute
            (
                srcMaps,
                localRanks,     // per src map which rank represents local data
                coarseComm,
                newToOldRanks,  // destination rank to source ranks
                srcStartOfLocal,
                srcCompactMaps
            )
        );


        // Assemble tgtAddress
        tgtAddress.setSize(nTgt);
        if (tgtAddress.size())
        {
            label alli = 0;
            forAll(allInterfaces, inti)
            {
                if (allInterfaces.set(inti))
                {
                    const auto& intf = refCast<const cyclicACMIGAMGInterface>
                    (
                        allInterfaces[inti]
                    );

                    if (!intf.amiPtr_)
                    {
                        continue;
                    }

                    const auto& AMI = intf.AMI();
                    const auto& tgtSlots = AMI.tgtAddress();
                    const label localSize =
                        srcStartOfLocal[inti+1]
                      - srcStartOfLocal[inti];

                    forAll(tgtSlots, tgti)
                    {
                        // Append old slots: copy old values and adapt
                        auto& newSlots = tgtAddress[alli++];
                        newSlots = tgtSlots[tgti];

                        // Renumber to new indices
                        mapDistributeBase::renumberMap
                        (
                            newSlots,
                            localSize,
                            srcStartOfLocal[inti],
                            srcCompactMaps[inti],
                            srcHasFlip //hasFlip
                        );

                        for (const label slot : newSlots)
                        {
                            if
                            (
                                slot < 0
                             || slot >= srcToTgtMap().constructSize()
                            )
                            {
                                FatalErrorInFunction << " newSlots:"
                                    << newSlots << exit(FatalError);
                            }
                        }
                    }
                }
            }

            if (nTgt != alli)
            {
                FatalErrorInFunction << "nTgt:" << nTgt
                    << " alli:" << alli << exit(FatalError);
            }
        }

        // Offsets for slots into results of tgtToSrcMap
        labelList tgtStartOfLocal;
        List<Map<label>> tgtCompactMaps;

        tgtToSrcMap.reset
        (
            new mapDistribute
            (
                tgtMaps,
                localRanks,
                coarseComm,
                newToOldRanks,
                tgtStartOfLocal,
                tgtCompactMaps
            )
        );


        // Assemble srcAddress
        srcAddress.setSize(nSrc);
        if (srcAddress.size())
        {
            label alli = 0;
            forAll(allInterfaces, inti)
            {
                if (allInterfaces.set(inti))
                {
                    const auto& intf = refCast<const cyclicACMIGAMGInterface>
                    (
                        allInterfaces[inti]
                    );

                    if (!intf.amiPtr_)
                    {
                        continue;
                    }

                    const auto& AMI = intf.AMI();
                    const auto& srcSlots = AMI.srcAddress();
                    const label localSize =
                        tgtStartOfLocal[inti+1]
                      - tgtStartOfLocal[inti];

                    forAll(srcSlots, srci)
                    {
                        // Append old slots: copy old values and adapt
                        auto& newSlots = srcAddress[alli++];
                        newSlots = srcSlots[srci];
                        // Renumber to new indices
                        mapDistributeBase::renumberMap
                        (
                            newSlots,
                            localSize,
                            tgtStartOfLocal[inti],
                            tgtCompactMaps[inti],
                            tgtHasFlip
                        );

                        for (const label slot : newSlots)
                        {
                            if
                            (
                                slot < 0
                             || slot >= tgtToSrcMap().constructSize()
                            )
                            {
                                FatalErrorInFunction << " newSlots:"
                                    << newSlots << exit(FatalError);
                            }
                        }
                    }
                }
            }

            if (nSrc != alli)
            {
                FatalErrorInFunction << "nSrc:" << nSrc
                    << " alli:" << alli << exit(FatalError);
            }
        }


        // Clean up: if no remote elements sent/received mark as
        // non-distributed. We could do this at the start but this
        // needs to take all the internal transport into account. Easier
        // (but less efficient) to do afterwards now all is compacted.
        {
            const auto& map = srcToTgtMap().subMap();

            bool usesRemote = false;
            forAll(map, proci)
            {
                if (proci != myAgglom)
                {
                    const auto& ss = srcToTgtMap().subMap()[proci];
                    const auto& sc = srcToTgtMap().constructMap()[proci];
                    const auto& ts = tgtToSrcMap().subMap()[proci];
                    const auto& tc = tgtToSrcMap().constructMap()[proci];

                    if (ss.size() || sc.size() || ts.size() || tc.size())
                    {
                        usesRemote = true;
                        break;
                    }
                }
            }

            // We can't have a single rank become fully-local since we
            // expect singlePatchProc to be synchronised. So make sure all
            // have become local

            if (!returnReduceOr(usesRemote, coarseComm))
            {
                DebugPout<< "** making fully local on new rank "
                    << myAgglom << " in comm:" << coarseComm << endl;
                singlePatchProc = myAgglom;
                srcToTgtMap.clear();
                tgtToSrcMap.clear();
            }
        }
    }
    else
    {
        // src/tgt address are straight indices

        srcAddress.setSize(nSrc);
        tgtAddress.setSize(nTgt);

        nSrc = 0;
        nTgt = 0;
        forAll(allInterfaces, inti)
        {
            if (allInterfaces.set(inti))
            {
                const auto& intf = refCast<const cyclicACMIGAMGInterface>
                (
                    allInterfaces[inti]
                );

                if (!intf.amiPtr_)
                {
                    continue;
                }

                const auto& AMI = intf.AMI();

                const auto& srcA = AMI.srcAddress();
                if (srcAddress.size())
                {
                    label srci = nSrc;
                    forAll(srcA, i)
                    {
                        srcAddress[srci++] = srcA[i]+nTgt;
                    }
                }

                const auto& tgtA = AMI.tgtAddress();
                if (tgtAddress.size())
                {
                    label tgti = nTgt;
                    forAll(tgtA, i)
                    {
                        tgtAddress[tgti++] = tgtA[i]+nSrc;
                    }
                }

                nSrc += srcA.size();
                nTgt += tgtA.size();
            }
        }
    }

    srcWeights.setSize(nSrc);
    if (hasSrcMagSf)
    {
        srcMagSf.setSize(nSrc);
    }
    if (hasSrcCentroids)
    {
        srcCentroids.setSize(nSrc);
    }
    tgtWeights.setSize(nTgt);
    if (hasTgtMagSf)
    {
        tgtMagSf.setSize(nTgt);
    }


    // Append individual data
    nSrc = 0;
    nTgt = 0;
    forAll(allInterfaces, inti)
    {
        if (allInterfaces.set(inti))
        {
            const auto& intf = refCast<const cyclicACMIGAMGInterface>
            (
                allInterfaces[inti]
            );

            if (!intf.amiPtr_)
            {
                continue;
            }

            const auto& AMI = intf.AMI();

            const auto& srcA = AMI.srcAddress();
            {
                // weights
                SubList<scalarList>(srcWeights, srcA.size(), nSrc) =
                    AMI.srcWeights();

                // magSf
                if (hasSrcMagSf)
                {
                    SubList<scalar>(srcMagSf, srcA.size(), nSrc) =
                        AMI.srcMagSf();
                }

                // centroids
                if (hasSrcCentroids)
                {
                    SubList<pointList>(srcCentroids, srcA.size(), nSrc) =
                        AMI.srcCentroids();
                }
            }

            const auto& tgtA = AMI.tgtAddress();
            {
                // weights
                SubList<scalarList>(tgtWeights, tgtA.size(), nTgt) =
                    AMI.tgtWeights();

                if (hasTgtMagSf)
                {
                    SubList<scalar>(tgtMagSf, tgtA.size(), nTgt) =
                        AMI.tgtMagSf();
                }
            }

            nSrc += srcA.size();
            nTgt += tgtA.size();
        }
    }


    // Construct with same arguments as original
    amiPtr_.reset
    (
        new AMIPatchToPatchInterpolation
        (
            requireMatch,
            reverseTarget,
            lowWeightCorrection
        )
    );
    amiPtr_().comm(coarseComm),
    amiPtr_().reset
    (
        std::move(srcToTgtMap),
        std::move(tgtToSrcMap),
        std::move(srcAddress),
        std::move(srcWeights),
        std::move(tgtAddress),
        std::move(tgtWeights),
        singlePatchProc
    );
    amiPtr_().srcMagSf() = std::move(srcMagSf);
    amiPtr_().srcCentroids() = std::move(srcCentroids);
    amiPtr_().tgtMagSf() = std::move(tgtMagSf);

    //Pout<< "** constructed new ami:"
    //    << " comm:" << amiPtr_().comm()
    //    << " srcMap.comm:" << amiPtr_().srcMap().comm()
    //    << " tgtMap.comm:" << amiPtr_().tgtMap().comm()
    //    << endl;
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::tmp<Foam::labelField> Foam::cyclicACMIGAMGInterface::internalFieldTransfer
(
    const Pstream::commsTypes commsType,
    const labelUList& iF
) const
{
    const cyclicACMIGAMGInterface& nbr =
        dynamic_cast<const cyclicACMIGAMGInterface&>(neighbPatch());
    const labelUList& nbrFaceCells = nbr.faceCells();

    auto tpnf = tmp<labelField>::New(nbrFaceCells.size());
    labelField& pnf = tpnf.ref();

    forAll(pnf, facei)
    {
        pnf[facei] = iF[nbrFaceCells[facei]];
    }

    return tpnf;
}


void Foam::cyclicACMIGAMGInterface::write(Ostream& os) const
{
    GAMGInterface::write(os);

    const bool hasAMI = bool(amiPtr_);

    os  << token::SPACE << neighbPatchID_
        << token::SPACE << owner_
        << token::SPACE << forwardT_
        << token::SPACE << reverseT_
        << token::SPACE << hasAMI;

    if (hasAMI)
    {
        os  << token::SPACE;
        AMI().writeData(os);

        // Write processors in communicator
        const label comm = AMI().comm();

        if (comm != -1)
        {
            os  << token::SPACE
                << UPstream::myProcNo(comm);
        }
    }
}


// ************************************************************************* //