Merge branch 'feature-distributedDIC' into 'develop'

Feature distributed DIC/DILU preconditioners

See merge request Development/openfoam!608

src/meshTools/Make/files

@@ -309,6 +309,11 @@ $(ACMICycPatches)/cyclicACMIPointPatchField/cyclicACMIPointPatchFields.C
 PeriodicAMICycPatches=$(AMI)/patches/cyclicPeriodicAMI
 $(PeriodicAMICycPatches)/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
 
+extraPreconditioners=matrices/lduMatrix/preconditioners/distributedDILUPreconditioner
+$(extraPreconditioners)/distributedDILUPreconditioner.C
+$(extraPreconditioners)/distributedDICPreconditioner.C
+$(extraPreconditioners)/processorColour.C
+
 multiWorld/multiWorldConnectionsObject.C
 
 mappedPatches/mappedPolyPatch/mappedPatchBase.C

src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDICPreconditioner.C

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "distributedDICPreconditioner.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(distributedDICPreconditioner, 0);
+
+    lduMatrix::preconditioner::
+        addsymMatrixConstructorToTable<distributedDICPreconditioner>
+        adddistributedDICPreconditionerSymMatrixConstructorToTable_;
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::distributedDICPreconditioner::forwardInternalDiag
+(
+    solveScalarField& rD,
+    const label colouri
+) const
+{
+    // Add forward constributions to diagonal
+
+    const auto& matrix = solver_.matrix();
+    const auto& lduAddr = matrix.lduAddr();
+
+    const label* const __restrict__ uPtr = lduAddr.upperAddr().begin();
+    const label* const __restrict__ lPtr = lduAddr.lowerAddr().begin();
+    const scalar* const __restrict__ upperPtr = matrix.upper().begin();
+
+    const label nFaces = matrix.upper().size();
+    if (cellColourPtr_.valid())
+    {
+        const auto& cellColour = cellColourPtr_();
+        for (label face=0; face<nFaces; face++)
+        {
+            const label cell = lPtr[face];
+            if (cellColour[cell] == colouri)
+            {
+                rD[uPtr[face]] -= upperPtr[face]*upperPtr[face]/rD[cell];
+            }
+        }
+    }
+    else
+    {
+        for (label face=0; face<nFaces; face++)
+        {
+            rD[uPtr[face]] -= upperPtr[face]*upperPtr[face]/rD[lPtr[face]];
+        }
+    }
+}
+
+
+void Foam::distributedDICPreconditioner::forwardInternal
+(
+    solveScalarField& wA,
+    const label colouri
+) const
+{
+    const auto& matrix = solver_.matrix();
+    const auto& lduAddr = matrix.lduAddr();
+
+    solveScalar* __restrict__ wAPtr = wA.begin();
+    const solveScalar* __restrict__ rDPtr = rD_.begin();
+
+    const label* const __restrict__ uPtr = lduAddr.upperAddr().begin();
+    const label* const __restrict__ lPtr = lduAddr.lowerAddr().begin();
+    const scalar* const __restrict__ upperPtr = matrix.upper().begin();
+
+    const label nFaces = matrix.upper().size();
+    if (cellColourPtr_.valid())
+    {
+        const auto& cellColour = cellColourPtr_();
+        for (label face=0; face<nFaces; face++)
+        {
+            const label cell = lPtr[face];
+            if (cellColour[cell] == colouri)
+            {
+                wAPtr[uPtr[face]] -=
+                    rDPtr[uPtr[face]]*upperPtr[face]*wAPtr[cell];
+            }
+        }
+    }
+    else
+    {
+        for (label face=0; face<nFaces; face++)
+        {
+            wAPtr[uPtr[face]] -=
+                rDPtr[uPtr[face]]*upperPtr[face]*wAPtr[lPtr[face]];
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::distributedDICPreconditioner::distributedDICPreconditioner
+(
+    const lduMatrix::solver& sol,
+    const dictionary& dict
+)
+:
+    distributedDILUPreconditioner(sol, dict)
+{}
+
+
+// ************************************************************************* //

src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDICPreconditioner.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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/>.
+
+Class
+    Foam::distributedDICPreconditioner
+
+Group
+    grpLduMatrixPreconditioners
+
+Description
+    Version of DICpreconditioner that uses preconditioning across processor
+    (and coupled) boundaries. Based on 'Parallel Preconditioners' chapter from
+    'Iterative Methods for Sparse Linear Systems' by Yousef Saad.
+
+    Leaves out the handling of boundary nodes after internal nodes since
+    probably not beneficial (so no overlap of comms and internal calculation)
+
+    By default preconditions across coupled boundaries (currently only
+    tested for cyclicAMI). This can be disabled with the 'coupled' setting
+
+        solver          PCG;
+        preconditioner
+        {
+            preconditioner  distributedDIC;
+            coupled         false;
+        }
+
+    The cyclicAMI boundary behaviour will only work if
+    - running non-parallel or
+    - different sides of cyclicAMI run on different processors i.e.
+      there is no processor which has cells on both owner and neighbour
+      of the patch pair.
+
+See Also
+    Foam::DICPreconditioner
+    Foam::distributedDILUPreconditioner
+
+SourceFiles
+    distributedDICPreconditioner.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef distributedDICPreconditioner_H
+#define distributedDICPreconditioner_H
+
+#include "distributedDILUPreconditioner.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                Class distributedDICPreconditioner Declaration
+\*---------------------------------------------------------------------------*/
+
+class distributedDICPreconditioner
+:
+    public distributedDILUPreconditioner
+{
+protected:
+
+    // Protected member functions
+
+        //- Update diagonal for colour
+        virtual void forwardInternalDiag
+        (
+            solveScalarField& rD,
+            const label colouri
+        ) const;
+
+        //- Update preconditioned variable walking forward on internal faces
+        virtual void forwardInternal
+        (
+            solveScalarField& wA,
+            const label colouri
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("distributedDIC");
+
+
+    // Constructors
+
+        //- Construct from matrix components and preconditioner solver controls
+        distributedDICPreconditioner
+        (
+            const lduMatrix::solver&,
+            const dictionary& solverControlsUnused
+        );
+
+
+    //- Destructor
+    virtual ~distributedDICPreconditioner() = default;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDILUPreconditioner.C

@@ -0,0 +1,897 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "distributedDILUPreconditioner.H"
+#include "processorLduInterface.H"
+#include "cyclicLduInterface.H"
+#include "cyclicAMILduInterface.H"
+#include "processorColour.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(distributedDILUPreconditioner, 0);
+
+    lduMatrix::preconditioner::
+        addasymMatrixConstructorToTable<distributedDILUPreconditioner>
+        adddistributedDILUPreconditionerAsymMatrixConstructorToTable_;
+
+    const lduMesh* distributedDILUPreconditioner::meshPtr_ = nullptr;
+
+    autoPtr<labelList> distributedDILUPreconditioner::procColoursPtr_(nullptr);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::distributedDILUPreconditioner::updateMatrixInterfaces
+(
+    const bool add,
+    const FieldField<Field, solveScalar>& coupleCoeffs,
+    const labelList& selectedInterfaces,
+    const solveScalarField& psiif,
+    solveScalarField& result,
+    const direction cmpt
+) const
+{
+    const auto& matrix = solver_.matrix();
+    const auto& lduAddr = matrix.lduAddr();
+    const auto& interfaces = solver_.interfaces();
+
+    const label startOfRequests = Pstream::nRequests();
+
+    // From lduMatrix::initMatrixInterfaces
+
+    // Avoid any conflicts with inter-processor comms
+    UPstream::msgType() += 321;
+
+    for (const label interfacei : selectedInterfaces)
+    {
+        interfaces[interfacei].initInterfaceMatrixUpdate
+        (
+            result,
+            add,
+            lduAddr,
+            interfacei,
+            psiif,
+            coupleCoeffs[interfacei],
+            cmpt,
+            UPstream::commsTypes::nonBlocking
+        );
+    }
+
+    UPstream::waitRequests(startOfRequests);
+
+    // Consume
+    for (const label interfacei : selectedInterfaces)
+    {
+        interfaces[interfacei].updateInterfaceMatrix
+        (
+            result,
+            add,
+            lduAddr,
+            interfacei,
+            psiif,
+            coupleCoeffs[interfacei],
+            cmpt,
+            UPstream::commsTypes::nonBlocking
+        );
+    }
+
+    UPstream::msgType() -= 321;
+}
+
+
+void Foam::distributedDILUPreconditioner::sendGlobal
+(
+    const labelList& selectedInterfaces,
+    solveScalarField& psi,
+    const label colouri
+) const
+{
+    const auto& interfaces = solver_.interfaces();
+
+    if (selectedInterfaces.size())
+    {
+        // Save old data
+        FieldField<Field, solveScalar> one(interfaces.size());
+        FieldField<Field, solveScalar> old(interfaces.size());
+        for (const label inti : selectedInterfaces)
+        {
+            const auto& intf = interfaces[inti].interface();
+            const auto& fc = intf.faceCells();
+            one.set(inti, new solveScalarField(fc.size(), 1.0));
+            old.set(inti, new solveScalarField(psi, intf.faceCells()));
+        }
+        updateMatrixInterfaces
+        (
+            false,              // add to psi
+            one,
+            selectedInterfaces,
+            psi,                // send data
+            psi,                // result
+            0                   // cmpt
+        );
+
+        if (!colourBufs_.set(colouri))
+        {
+            colourBufs_.set
+            (
+                colouri,
+                new FieldField<Field, solveScalar>
+                (
+                    interfaces.size()
+                )
+            );
+        }
+        auto& colourBuf = colourBufs_[colouri];
+        colourBuf.setSize(interfaces.size());
+        for (const label inti : selectedInterfaces)
+        {
+            const auto& intf = interfaces[inti].interface();
+            const auto& fc = intf.faceCells();
+            if (!colourBuf.set(inti))
+            {
+                colourBuf.set(inti, new Field<solveScalar>(fc.size()));
+            }
+            auto& cb = colourBuf[inti];
+            //cb.resize_nocopy(fc.size());
+            auto& oldValues = old[inti];
+
+            forAll(cb, face)
+            {
+                const label cell = fc[face];
+                // Store change in boundary value
+                cb[face] = psi[cell]-oldValues[face];
+                // Restore old value
+                std::swap(psi[cell], oldValues[face]);
+            }
+        }
+    }
+}
+
+
+void Foam::distributedDILUPreconditioner::receive
+(
+    const labelList& selectedInterfaces,
+    DynamicList<UPstream::Request>& requests
+) const
+{
+    const auto& interfaces = solver_.interfaces();
+    const auto& interfaceBouCoeffs = solver_.interfaceBouCoeffs();
+
+    // Start reads
+    for (const label inti : selectedInterfaces)
+    {
+        const auto& intf = interfaces[inti].interface();
+        const auto* ppp = isA<const processorLduInterface>(intf);
+
+        auto& recvBuf = recvBufs_[inti];
+        recvBuf.resize_nocopy(interfaceBouCoeffs[inti].size());
+
+        requests.push_back(UPstream::Request());
+        UIPstream::read
+        (
+            requests.back(),
+            ppp->neighbProcNo(),
+            recvBuf.data_bytes(),
+            recvBuf.size_bytes(),
+            ppp->tag()+70,          // random offset
+            ppp->comm()
+        );
+    }
+}
+
+
+void Foam::distributedDILUPreconditioner::send
+(
+    const labelList& selectedInterfaces,
+    const solveScalarField& psiInternal,
+    DynamicList<UPstream::Request>& requests
+) const
+{
+    const auto& interfaces = solver_.interfaces();
+
+    // Start writes
+    for (const label inti : selectedInterfaces)
+    {
+        const auto& intf = interfaces[inti].interface();
+        const auto* ppp = isA<const processorLduInterface>(intf);
+        const auto& faceCells = intf.faceCells();
+
+        auto& sendBuf = sendBufs_[inti];
+
+        sendBuf.resize_nocopy(faceCells.size());
+        forAll(faceCells, face)
+        {
+            sendBuf[face] = psiInternal[faceCells[face]];
+        }
+
+        requests.push_back(UPstream::Request());
+        UOPstream::write
+        (
+            requests.back(),
+            ppp->neighbProcNo(),
+            sendBuf.cdata_bytes(),
+            sendBuf.size_bytes(),
+            ppp->tag()+70,          // random offset
+            ppp->comm()
+        );
+    }
+}
+
+
+void Foam::distributedDILUPreconditioner::wait
+(
+    DynamicList<UPstream::Request>& requests
+) const
+{
+    UPstream::waitRequests(requests);
+    requests.clear();
+}
+
+
+// Diagonal calculation
+// ~~~~~~~~~~~~~~~~~~~~
+
+void Foam::distributedDILUPreconditioner::addInterfaceDiag
+(
+    solveScalarField& rD,
+    const label inti,
+    const Field<solveScalar>& recvBuf
+) const
+{
+    const auto& interfaces = solver_.interfaces();
+    const auto& interfaceBouCoeffs = solver_.interfaceBouCoeffs();
+    const auto& interfaceIntCoeffs = solver_.interfaceIntCoeffs();
+
+    const auto& intf = interfaces[inti].interface();
+    // TBD: do not use patch faceCells but passed-in addressing?
+    const auto& faceCells = intf.faceCells();
+    const auto& bc = interfaceBouCoeffs[inti];
+    const auto& ic = interfaceIntCoeffs[inti];
+
+    forAll(recvBuf, face)
+    {
+        // Note:interfaceBouCoeffs, interfaceIntCoeffs on the receiving
+        //      side are negated
+        rD[faceCells[face]] -= bc[face]*ic[face]/recvBuf[face];
+    }
+}
+
+
+void Foam::distributedDILUPreconditioner::forwardInternalDiag
+(
+    solveScalarField& rD,
+    const label colouri
+) const
+{
+    // Add forward constributions to diagonal
+
+    const auto& matrix = solver_.matrix();
+    const auto& lduAddr = matrix.lduAddr();
+
+    const label* const __restrict__ uPtr = lduAddr.upperAddr().begin();
+    const label* const __restrict__ lPtr = lduAddr.lowerAddr().begin();
+    const scalar* const __restrict__ upperPtr = matrix.upper().begin();
+    const scalar* const __restrict__ lowerPtr = matrix.lower().begin();
+
+    const label nFaces = matrix.upper().size();
+    if (cellColourPtr_.valid())
+    {
+        const auto& cellColour = cellColourPtr_();
+        for (label face=0; face<nFaces; face++)
+        {
+            const label cell = lPtr[face];
+            if (cellColour[cell] == colouri)
+            {
+                rD[uPtr[face]] -= upperPtr[face]*lowerPtr[face]/rD[cell];
+            }
+        }
+    }
+    else
+    {
+        for (label face=0; face<nFaces; face++)
+        {
+            rD[uPtr[face]] -= upperPtr[face]*lowerPtr[face]/rD[lPtr[face]];
+        }
+    }
+}
+
+
+void Foam::distributedDILUPreconditioner::addInterface
+(
+    solveScalarField& wA,
+    const label inti,
+    const Field<solveScalar>& recvBuf
+) const
+{
+    const auto& interfaces = solver_.interfaces();
+    const auto& interfaceBouCoeffs = solver_.interfaceBouCoeffs();
+
+    const auto& intf = interfaces[inti].interface();
+    const auto& faceCells = intf.faceCells();
+    const auto& bc = interfaceBouCoeffs[inti];
+    const solveScalar* __restrict__ rDPtr = rD_.begin();
+    solveScalar* __restrict__ wAPtr = wA.begin();
+
+    forAll(recvBuf, face)
+    {
+        // Note: interfaceBouCoeffs is -upperPtr
+        const label cell = faceCells[face];
+        wAPtr[cell] += rDPtr[cell]*bc[face]*recvBuf[face];
+    }
+}
+
+
+void Foam::distributedDILUPreconditioner::forwardInternal
+(
+    solveScalarField& wA,
+    const label colouri
+) const
+{
+    const auto& matrix = solver_.matrix();
+    const auto& lduAddr = matrix.lduAddr();
+
+    solveScalar* __restrict__ wAPtr = wA.begin();
+    const solveScalar* __restrict__ rDPtr = rD_.begin();
+
+    const label* const __restrict__ uPtr = lduAddr.upperAddr().begin();
+    const label* const __restrict__ lPtr = lduAddr.lowerAddr().begin();
+    const label* const __restrict__ losortPtr = lduAddr.losortAddr().begin();
+    const scalar* const __restrict__ lowerPtr = matrix.lower().begin();
+
+    const label nFaces = matrix.upper().size();
+    if (cellColourPtr_.valid())
+    {
+        const auto& cellColour = cellColourPtr_();
+        for (label face=0; face<nFaces; face++)
+        {
+            const label sface = losortPtr[face];
+            const label cell = lPtr[sface];
+            if (cellColour[cell] == colouri)
+            {
+                wAPtr[uPtr[sface]] -=
+                    rDPtr[uPtr[sface]]*lowerPtr[sface]*wAPtr[cell];
+            }
+        }
+    }
+    else
+    {
+        for (label face=0; face<nFaces; face++)
+        {
+            const label sface = losortPtr[face];
+            wAPtr[uPtr[sface]] -=
+                rDPtr[uPtr[sface]]*lowerPtr[sface]*wAPtr[lPtr[sface]];
+        }
+    }
+}
+
+
+void Foam::distributedDILUPreconditioner::backwardInternal
+(
+    solveScalarField& wA,
+    const label colouri
+) const
+{
+    const auto& matrix = solver_.matrix();
+    const auto& lduAddr = matrix.lduAddr();
+
+    solveScalar* __restrict__ wAPtr = wA.begin();
+    const solveScalar* __restrict__ rDPtr = rD_.begin();
+
+    const label* const __restrict__ uPtr = lduAddr.upperAddr().begin();
+    const label* const __restrict__ lPtr = lduAddr.lowerAddr().begin();
+    const scalar* const __restrict__ upperPtr = matrix.upper().begin();
+
+    const label nFacesM1 = matrix.upper().size() - 1;
+
+    if (cellColourPtr_.valid())
+    {
+        const auto& cellColour = cellColourPtr_();
+        for (label face=nFacesM1; face>=0; face--)
+        {
+            const label cell = uPtr[face];
+            if (cellColour[cell] == colouri)
+            {
+                // Note: lower cell guaranteed in same colour
+                wAPtr[lPtr[face]] -=
+                    rDPtr[lPtr[face]]*upperPtr[face]*wAPtr[cell];
+            }
+        }
+    }
+    else
+    {
+        for (label face=nFacesM1; face>=0; face--)
+        {
+            wAPtr[lPtr[face]] -=
+                rDPtr[lPtr[face]]*upperPtr[face]*wAPtr[uPtr[face]];
+        }
+    }
+}
+
+
+void Foam::distributedDILUPreconditioner::calcReciprocalD
+(
+    solveScalarField& rD
+) const
+{
+    const auto& matrix = solver_.matrix();
+
+    // Make sure no outstanding receives
+    wait(lowerRecvRequests_);
+
+    // Start reads (into recvBufs_)
+    receive(lowerNbrs_, lowerRecvRequests_);
+
+    // Start with diagonal
+    const scalarField& diag = matrix.diag();
+    rD.resize_nocopy(diag.size());
+    std::copy(diag.begin(), diag.end(), rD.begin());
+
+
+    // Subtract coupled contributions
+    {
+        // Wait for finish. Received result in recvBufs
+        wait(lowerRecvRequests_);
+
+        for (const label inti : lowerNbrs_)
+        {
+            addInterfaceDiag(rD, inti, recvBufs_[inti]);
+        }
+    }
+
+
+    // - divide the internal mesh into domains/colour, similar to domain
+    //   decomposition
+    // - we could use subsetted bits of the mesh or just loop over only
+    //   the cells of the current domain
+    // - a domain only uses boundary values of lower numbered domains
+    // - this also limits the interfaces that need to be evaluated since
+    //   we assume an interface only changes local faceCells and these
+    //   are all of the same colour
+
+    for (label colouri = 0; colouri < nColours_; colouri++)
+    {
+        if (cellColourPtr_.valid())// && colourBufs_.set(colouri))
+        {
+            for (const label inti : lowerGlobalRecv_[colouri])
+            {
+                addInterfaceDiag(rD, inti, colourBufs_[colouri][inti]);
+            }
+        }
+
+        forwardInternalDiag(rD, colouri);
+
+        // Store effect of exchanging rD to higher interfaces in colourBufs_
+        if (cellColourPtr_.valid())
+        {
+            sendGlobal(higherGlobalSend_[colouri], rD, higherColour_[colouri]);
+        }
+    }
+
+    // Make sure no outstanding sends
+    wait(higherSendRequests_);
+
+    // Start writes of rD (using sendBufs)
+    send(higherNbrs_, rD, higherSendRequests_);
+
+    // Calculate the reciprocal of the preconditioned diagonal
+    const label nCells = rD.size();
+
+    for (label cell=0; cell<nCells; cell++)
+    {
+        rD[cell] = 1.0/rD[cell];
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::distributedDILUPreconditioner::distributedDILUPreconditioner
+(
+    const lduMatrix::solver& sol,
+    const dictionary& dict
+)
+:
+    lduMatrix::preconditioner(sol),
+    coupled_(dict.getOrDefault<bool>("coupled", true, keyType::LITERAL)),
+    rD_(sol.matrix().diag().size())
+{
+    const lduMesh& mesh = sol.matrix().mesh();
+    const auto& interfaces = sol.interfaces();
+    const auto& interfaceBouCoeffs = sol.interfaceBouCoeffs();
+
+    // Allocate buffers
+    // ~~~~~~~~~~~~~~~~
+
+    sendBufs_.setSize(interfaces.size());
+    recvBufs_.setSize(interfaces.size());
+    forAll(interfaceBouCoeffs, inti)
+    {
+        if (interfaces.set(inti))
+        {
+            const auto& bc = interfaceBouCoeffs[inti];
+            sendBufs_.set(inti, new scalarField(bc.size(), Zero));
+            recvBufs_.set(inti, new scalarField(bc.size(), Zero));
+        }
+    }
+
+
+    // Determine processor colouring
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    //const processorColour& colours = processorColour::New(mesh);
+    //const label myColour = colours[Pstream::myProcNo()];
+    if (&mesh != meshPtr_)
+    {
+        procColoursPtr_.reset(new labelList(Pstream::nProcs(mesh.comm())));
+        const label nProcColours =
+            processorColour::colour(mesh, procColoursPtr_());
+
+        if (debug)
+        {
+            Info<< typeName << " : calculated processor colours:"
+                << nProcColours << endl;
+        }
+
+        meshPtr_ = &mesh;
+    }
+    const labelList& procColours = procColoursPtr_();
+
+    const label myColour = procColours[Pstream::myProcNo(mesh.comm())];
+
+    bool haveGlobalCoupled = false;
+    bool haveDistributedAMI = false;
+    forAll(interfaces, inti)
+    {
+        if (interfaces.set(inti))
+        {
+            const auto& intf = interfaces[inti].interface();
+            const auto* ppp = isA<const processorLduInterface>(intf);
+            if (ppp)
+            {
+                const label nbrColour = procColours[ppp->neighbProcNo()];
+                //const label nbrColour = ppp->neighbProcNo();
+                if (nbrColour < myColour)
+                {
+                    lowerNbrs_.append(inti);
+                }
+                else if (nbrColour > myColour)
+                {
+                    higherNbrs_.append(inti);
+                }
+                else
+                {
+                    WarningInFunction
+                        << "weird processorLduInterface"
+                        << " from " << ppp->myProcNo()
+                        << " to " << ppp->neighbProcNo()
+                        << endl;
+                }
+            }
+            else //if (isA<const cyclicAMILduInterface>(intf))
+            {
+                haveGlobalCoupled = true;
+
+                const auto* AMIpp = isA<const cyclicAMILduInterface>(intf);
+                if (AMIpp)
+                {
+                    //const auto& AMI =
+                    //(
+                    //    AMIpp->owner()
+                    //  ? AMIpp->AMI()
+                    //  : AMIpp->neighbPatch().AMI()
+                    //);
+                    //haveDistributedAMI = AMI.distributed();
+
+                    if (debug)
+                    {
+                        Info<< typeName << " : interface " << inti
+                            << " of type " << intf.type()
+                            << " is distributed:" << haveDistributedAMI << endl;
+                    }
+                }
+            }
+        }
+    }
+
+
+
+    if (debug)
+    {
+        Info<< typeName
+            << " : lower proc interfaces:" << flatOutput(lowerNbrs_)
+            << " higher proc interfaces:" << flatOutput(higherNbrs_)
+            << endl;
+    }
+
+
+    // Determine local colouring/zoneing
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Start off with all cells colour 0
+    nColours_ = 1;
+    cellColourPtr_.clear();
+    if (coupled_ && haveGlobalCoupled)
+    {
+        labelList patchToColour;
+        cellColourPtr_.reset(new labelList(mesh.lduAddr().size()));
+        //if (haveDistributedAMI)
+        //{
+        //    nColours_ = 1;
+        //    cellColourPtr_() = 0;
+        //    patchToColour = labelList(interfaces.size(), 0);
+        //}
+        //else
+        {
+            nColours_ = processorColour::cellColour
+            (
+                mesh,
+                cellColourPtr_(),
+                patchToColour
+            );
+        }
+
+        lowerGlobalRecv_.resize_nocopy(nColours_);
+        lowerGlobalSend_.resize_nocopy(nColours_);
+        lowerColour_.setSize(nColours_, -1);
+        higherGlobalRecv_.resize_nocopy(nColours_);
+        higherGlobalSend_.resize_nocopy(nColours_);
+        higherColour_.setSize(nColours_, -1);
+        colourBufs_.setSize(nColours_);
+
+        forAll(interfaces, inti)
+        {
+            if (interfaces.set(inti))
+            {
+                const auto& intf = interfaces[inti].interface();
+
+                label nbrInti = -1;
+                const auto* AMIpp = isA<const cyclicAMILduInterface>(intf);
+                if (AMIpp)
+                {
+                    nbrInti = AMIpp->neighbPatchID();
+                }
+                const auto* cycpp = isA<const cyclicLduInterface>(intf);
+                if (cycpp)
+                {
+                    nbrInti = cycpp->neighbPatchID();
+                }
+
+                if (nbrInti != -1)
+                {
+                    const label colouri = patchToColour[inti];
+                    const label nbrColouri = patchToColour[nbrInti];
+
+                    if
+                    (
+                        (haveDistributedAMI && (inti < nbrInti))
+                     || (!haveDistributedAMI && (colouri < nbrColouri))
+                    )
+                    {
+                        // Send to higher
+                        higherGlobalSend_[colouri].append(nbrInti);
+                        higherColour_[colouri] = nbrColouri;
+                        // Receive from higher
+                        higherGlobalRecv_[colouri].append(inti);
+                    }
+                    else
+                    {
+                        // Send to lower
+                        lowerGlobalSend_[colouri].append(nbrInti);
+                        lowerColour_[colouri] = nbrColouri;
+                        // Receive from lower
+                        lowerGlobalRecv_[colouri].append(inti);
+                    }
+                }
+            }
+        }
+    }
+
+
+    if (debug)
+    {
+        Info<< typeName << " : local colours:" << nColours_ << endl;
+        Info<< incrIndent;
+        forAll(lowerGlobalRecv_, colouri)
+        {
+            const auto& lowerRecv = lowerGlobalRecv_[colouri];
+            if (lowerRecv.size())
+            {
+                const auto& lowerSend = lowerGlobalSend_[colouri];
+                const auto& lowerCol = lowerColour_[colouri];
+                Info<< indent
+                    << "Lower global interfaces for colour:" << colouri
+                    << " receiving from:" << flatOutput(lowerRecv)
+                    << " sending to:" << flatOutput(lowerSend)
+                    << " with colour:" << lowerCol << endl;
+            }
+        }
+        forAll(higherGlobalRecv_, colouri)
+        {
+            const auto& higherRecv = higherGlobalRecv_[colouri];
+            if (higherRecv.size())
+            {
+                const auto& higherSend = higherGlobalSend_[colouri];
+                const auto& higherCol = higherColour_[colouri];
+                Info<< indent
+                    << "Higher global interfaces for colour:" << colouri
+                    << " receiving from:" << flatOutput(higherRecv)
+                    << " sending to:" << flatOutput(higherSend)
+                    << " with colour:" << higherCol << endl;
+            }
+        }
+    }
+
+    calcReciprocalD(rD_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::distributedDILUPreconditioner::~distributedDILUPreconditioner()
+{
+    // Wait on all requests before storage is being taken down
+    // (could rely on construction order?)
+    wait(lowerSendRequests_);
+    wait(lowerRecvRequests_);
+    wait(higherSendRequests_);
+    wait(higherRecvRequests_);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::distributedDILUPreconditioner::precondition
+(
+    solveScalarField& wA,
+    const solveScalarField& rA,
+    const direction
+) const
+{
+    solveScalar* __restrict__ wAPtr = wA.begin();
+    const solveScalar* __restrict__ rAPtr = rA.begin();
+    const solveScalar* __restrict__ rDPtr = rD_.begin();
+
+    const label nCells = wA.size();
+
+
+    // Forward sweep
+    // ~~~~~~~~~~~~~
+
+    // Make sure no receives are still in flight (should not happen)
+    wait(lowerRecvRequests_);
+
+    // Start reads (into recvBufs)
+    receive(lowerNbrs_, lowerRecvRequests_);
+
+    // Initialise 'internal' cells
+    for (label cell=0; cell<nCells; cell++)
+    {
+        wAPtr[cell] = rDPtr[cell]*rAPtr[cell];
+    }
+
+    // Do 'halo' contributions from lower numbered procs
+    {
+        // Wait for finish. Received result in recvBufs
+        wait(lowerRecvRequests_);
+
+        for (const label inti : lowerNbrs_)
+        {
+            addInterface(wA, inti, recvBufs_[inti]);
+        }
+
+        for (label colouri = 0; colouri < nColours_; colouri++)
+        {
+            // Do non-processor boundaries for this colour
+            if (cellColourPtr_.valid())// && colourBufs_.set(colouri))
+            {
+                for (const label inti : lowerGlobalRecv_[colouri])
+                {
+                    addInterface(wA, inti, colourBufs_[colouri][inti]);
+                }
+            }
+
+            forwardInternal(wA, colouri);
+
+            // Store effect of exchanging rD to higher interfaces
+            // in colourBufs_
+            if (cellColourPtr_.valid())
+            {
+                sendGlobal
+                (
+                    higherGlobalSend_[colouri],
+                    wA,
+                    higherColour_[colouri]
+                );
+            }
+        }
+    }
+
+    // Make sure no outstanding sends from previous iteration
+    wait(higherSendRequests_);
+
+    // Start writes of wA (using sendBufs)
+    send(higherNbrs_, wA, higherSendRequests_);
+
+
+    // Backward sweep
+    // ~~~~~~~~~~~~~~
+
+    // Make sure no outstanding receives
+    wait(higherRecvRequests_);
+
+    // Start receives
+    receive(higherNbrs_, higherRecvRequests_);
+
+    {
+        // Wait until receives finished
+        wait(higherRecvRequests_);
+
+        for (const label inti : higherNbrs_)
+        {
+            addInterface(wA, inti, recvBufs_[inti]);
+        }
+
+        for (label colouri = nColours_-1; colouri >= 0; colouri--)
+        {
+            // Do non-processor boundaries for this colour
+            if (cellColourPtr_.valid())// && colourBufs_.set(colouri))
+            {
+                for (const label inti : higherGlobalRecv_[colouri])
+                {
+                    addInterface(wA, inti, colourBufs_[colouri][inti]);
+                }
+            }
+
+            backwardInternal(wA, colouri);
+
+            // Store effect of exchanging rD to higher interfaces
+            // in colourBufs_
+            if (cellColourPtr_.valid())
+            {
+                sendGlobal
+                (
+                    lowerGlobalSend_[colouri],
+                    wA,
+                    lowerColour_[colouri]
+                );
+            }
+        }
+    }
+
+    // Make sure no outstanding sends
+    wait(lowerSendRequests_);
+
+    // Start writes of wA (using sendBufs)
+    send(lowerNbrs_, wA, lowerSendRequests_);
+}
+
+
+// ************************************************************************* //

src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDILUPreconditioner.H

@@ -0,0 +1,274 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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/>.
+
+Class
+    Foam::distributedDILUPreconditioner
+
+Group
+    grpLduMatrixPreconditioners
+
+Description
+    Version of DILUpreconditioner that uses preconditioning across processor
+    (and coupled) boundaries. Based on 'Parallel Preconditioners' chapter from
+    'Iterative Methods for Sparse Linear Systems' by Yousef Saad.
+
+    Leaves out the handling of boundary nodes after internal nodes since
+    probably not beneficial (so no overlap of comms and internal calculation)
+
+    By default preconditions across coupled boundaries (currently only
+    tested for cyclicAMI). This can be disabled with the 'coupled' setting
+
+        solver          PCG;
+        preconditioner
+        {
+            preconditioner  distributedDILU;
+            coupled         false;
+        }
+
+    The cyclicAMI boundary behaviour will only work if
+    - running non-parallel or
+    - different sides of cyclicAMI run on different processors i.e.
+      there is no processor which has cells on both owner and neighbour
+      of the patch pair.
+
+See Also
+    Foam::DILUPreconditioner
+    Foam::distributedDICPreconditioner
+
+SourceFiles
+    distributedDILUPreconditioner.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef distributedDILUPreconditioner_H
+#define distributedDILUPreconditioner_H
+
+#include "lduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                Class distributedDILUPreconditioner Declaration
+\*---------------------------------------------------------------------------*/
+
+class distributedDILUPreconditioner
+:
+    public lduMatrix::preconditioner
+{
+protected:
+
+    // Protected data
+
+        //- Precondition across global coupled bc
+        const bool coupled_;
+
+
+        //- Processor interface buffers and colouring
+
+            //- Previous mesh
+            static const lduMesh* meshPtr_;
+
+            //- Previous processor colours
+            static autoPtr<labelList> procColoursPtr_;
+
+            //- Buffers for sending and receiving data
+            mutable FieldField<Field, solveScalar> sendBufs_;
+            mutable FieldField<Field, solveScalar> recvBufs_;
+            mutable DynamicList<UPstream::Request> recvRequests_;
+
+            //- Interfaces to lower coloured processors
+            DynamicList<label> lowerNbrs_;
+            mutable DynamicList<UPstream::Request> lowerSendRequests_;
+            mutable DynamicList<UPstream::Request> lowerRecvRequests_;
+
+            //- Interfaces to higher coloured processors
+            DynamicList<label> higherNbrs_;
+            mutable DynamicList<UPstream::Request> higherSendRequests_;
+            mutable DynamicList<UPstream::Request> higherRecvRequests_;
+
+
+        //- Local (cell) colouring from global interfaces
+
+            //- Colour/zone per cell
+            autoPtr<labelList> cellColourPtr_;
+
+            //- Number of colours (in case of multiple disconnected regions
+            //  in single mesh)
+            label nColours_;
+
+
+        //- Global interfaces. Per colour the interfaces that (might)
+        //- influence it
+
+            mutable PtrList<FieldField<Field, solveScalar>> colourBufs_;
+
+            //- Interfaces to non-processor lower coupled interfaces
+            List<DynamicList<label>> lowerGlobalRecv_;
+
+            //- Interfaces to non-processor lower coupled interfaces
+            List<DynamicList<label>> lowerGlobalSend_;
+
+            //- Corresponding destination colour (for lowerGlobal)
+            List<label> lowerColour_;
+
+            //- Interfaces to non-processor higher coupled interfaces
+            List<DynamicList<label>> higherGlobalRecv_;
+
+            //- Interfaces to non-processor higher coupled interfaces
+            List<DynamicList<label>> higherGlobalSend_;
+
+            //- Corresponding destination colour (for higherGlobal)
+            List<label> higherColour_;
+
+
+        //- The reciprocal preconditioned diagonal
+        solveScalarField rD_;
+
+
+    // Private Member Functions
+
+        //- Variant of lduMatrix::updateMatrixInterfaces on selected interfaces
+        void updateMatrixInterfaces
+        (
+            const bool add,
+            const FieldField<Field, solveScalar>& coupleCoeffs,
+            const labelList& selectedInterfaces,
+            const solveScalarField& psiif,
+            solveScalarField& result,
+            const direction cmpt
+        ) const;
+
+        //- Send (and store in colourBufs_[colouri]) the effect of
+        //  doing selectedInterfaces
+        void sendGlobal
+        (
+            const labelList& selectedInterfaces,
+            solveScalarField& psi,
+            const label colouri
+        ) const;
+
+        //- Start receiving in recvBufs_
+        void receive
+        (
+            const labelList& selectedInterfaces,
+            DynamicList<UPstream::Request>& requests
+        ) const;
+
+        //- Start sending sendBufs_
+        void send
+        (
+            const labelList& selectedInterfaces,
+            const solveScalarField& psiInternal,
+            DynamicList<UPstream::Request>& requests
+        ) const;
+
+        //- Wait for explicit requests
+        void wait(DynamicList<UPstream::Request>& requests) const;
+
+        //- Update diagonal for interface
+        virtual void addInterfaceDiag
+        (
+            solveScalarField& rD,
+            const label inti,
+            const Field<solveScalar>& recvBuf
+        ) const;
+
+        //- Update diagonal for all faces of a certain colour
+        virtual void forwardInternalDiag
+        (
+            solveScalarField& rD,
+            const label colouri
+        ) const;
+
+        //- Update preconditioned variable from interface
+        virtual void addInterface
+        (
+            solveScalarField& wA,
+            const label inti,
+            const Field<solveScalar>& recvBuf
+        ) const;
+
+        //- Update preconditioned variable walking forward on internal faces
+        virtual void forwardInternal
+        (
+            solveScalarField& wA,
+            const label colouri
+        ) const;
+
+        //- Update preconditioned variable walking backward on internal faces
+        virtual void backwardInternal
+        (
+            solveScalarField& wA,
+            const label colouri
+        ) const;
+
+        //- Calculate reciprocal of diagonal
+        virtual void calcReciprocalD(solveScalarField& rD) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("distributedDILU");
+
+
+    // Constructors
+
+        //- Construct from matrix components and preconditioner solver controls
+        distributedDILUPreconditioner
+        (
+            const lduMatrix::solver&,
+            const dictionary& solverControlsUnused
+        );
+
+
+    //- Destructor
+    virtual ~distributedDILUPreconditioner();
+
+
+    // Member Functions
+
+        //- Return wA the preconditioned form of residual rA
+        virtual void precondition
+        (
+            solveScalarField& wA,
+            const solveScalarField& rA,
+            const direction cmpt=0
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/processorColour.C

@@ -0,0 +1,363 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 M. Janssens
+-------------------------------------------------------------------------------
+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 "processorColour.H"
+#include "processorLduInterface.H"
+#include "processorTopologyNew.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(processorColour, 0);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::label Foam::processorColour::colour
+(
+    const lduMesh& mesh,
+    labelList& procColour
+)
+{
+    procColour.resize_nocopy(Pstream::nProcs(mesh.comm()));
+    procColour = -1;
+
+    // Re-use processor-topology analysis
+
+    labelListList procNeighbours(Pstream::nProcs(mesh.comm()));
+
+    // Fill my entry
+    {
+        const lduInterfacePtrsList patches = mesh.interfaces();
+
+        auto& procToProcs = procNeighbours[Pstream::myProcNo(mesh.comm())];
+        label n = 0;
+        forAll(patches, patchi)
+        {
+            if (patches.set(patchi))
+            {
+                if (isA<processorLduInterface>(patches[patchi]))
+                {
+                    n++;
+                }
+            }
+        }
+
+        procToProcs.resize_nocopy(n);
+        n = 0;
+        forAll(patches, patchi)
+        {
+            if (patches.set(patchi))
+            {
+                const auto* ppPtr = isA<processorLduInterface>(patches[patchi]);
+                if (ppPtr)
+                {
+                    procToProcs[n++] = ppPtr->neighbProcNo();
+                }
+            }
+        }
+    }
+    // Send to master
+    Pstream::gatherList(procNeighbours, UPstream::msgType(), mesh.comm());
+
+
+    // Use greedy algorithm for now
+    // (see e.g. https://iq.opengenus.org/graph-colouring-greedy-algorithm/)
+
+    if (Pstream::master(mesh.comm()))
+    {
+        DynamicList<label> front;
+        front.append(0);    // start from processor 0
+
+        DynamicList<label> newFront;
+        while (front.size())
+        {
+            //Pout<< "Front:" << front << endl;
+
+            newFront.clear();
+            for (const label proci : front)
+            {
+                if (procColour[proci] == -1)
+                {
+                    const labelList& nbrs = procNeighbours[proci];
+                    const UIndirectList<label> nbrColour(procColour, nbrs);
+
+                    for
+                    (
+                        label colouri = 0;
+                        colouri < Pstream::nProcs();
+                        colouri++
+                    )
+                    {
+                        if (!nbrColour.found(colouri))
+                        {
+                            procColour[proci] = colouri;
+                            for (label nbrProci : nbrs)
+                            {
+                                if (procColour[nbrProci] == -1)
+                                {
+                                    newFront.append(nbrProci);
+                                }
+                            }
+                            break;
+                        }
+                    }
+                }
+            }
+
+            front = std::move(newFront);
+        }
+    }
+
+    Pstream::broadcast(procColour, mesh.comm());
+
+
+    //if (false)
+    //{
+    //    volScalarField volColour
+    //    (
+    //        IOobject
+    //        (
+    //            "colour",
+    //            mesh.time().timeName(),
+    //            mesh,
+    //            IOobject::NO_READ,
+    //            IOobject::AUTO_WRITE,
+    //            false
+    //        ),
+    //        mesh,
+    //        dimensionedScalar(dimless, procColour[Pstream::myProcNo()]),
+    //        zeroGradientFvPatchScalarField::typeName
+    //    );
+    //    volColour.write();
+    //}
+
+    const label nColours = max(procColour)+1;
+
+    if (debug)
+    {
+        Info<< typeName << " : coloured " << Pstream::nProcs(mesh.comm())
+            << " processors with in total " << nColours << " colours" << endl;
+    }
+
+    return nColours;
+}
+
+
+void Foam::processorColour::walkFront
+(
+    const lduMesh& mesh,
+    DynamicList<label>& front,
+    labelList& cellColour
+)
+{
+    // Colour with the (min) coupled (global) patch
+
+    const lduAddressing& addr = mesh.lduAddr();
+    const label* const __restrict__ uPtr = addr.upperAddr().begin();
+    const label* const __restrict__ lPtr = addr.lowerAddr().begin();
+    const label* const __restrict__ ownStartPtr = addr.ownerStartAddr().begin();
+    const label* const __restrict__ losortStartAddrPtr =
+        addr.losortStartAddr().begin();
+    const label* const __restrict__ losortAddrPtr = addr.losortAddr().begin();
+
+    DynamicList<label> newFront;
+    while (true)
+    {
+        newFront.clear();
+        for (const label celli : front)
+        {
+            const label colouri = cellColour[celli];
+
+            {
+                const label fStart = ownStartPtr[celli];
+                const label fEnd = ownStartPtr[celli + 1];
+
+                for (label facei=fStart; facei<fEnd; facei++)
+                {
+                    const label nbr =
+                    (
+                        lPtr[facei] == celli
+                      ? uPtr[facei]
+                      : lPtr[facei]
+                    );
+                    if (cellColour[nbr] == -1)
+                    {
+                        cellColour[nbr] = colouri;
+                        newFront.append(nbr);
+                    }
+                }
+            }
+            {
+                const label fStart = losortStartAddrPtr[celli];
+                const label fEnd = losortStartAddrPtr[celli + 1];
+
+                for (label i=fStart; i<fEnd; i++)
+                {
+                    label facei = losortAddrPtr[i];
+                    const label nbr =
+                    (
+                        lPtr[facei] == celli
+                      ? uPtr[facei]
+                      : lPtr[facei]
+                    );
+                    if (cellColour[nbr] == -1)
+                    {
+                        cellColour[nbr] = colouri;
+                        newFront.append(nbr);
+                    }
+                }
+            }
+
+        }
+
+        if (newFront.empty())
+        {
+            break;
+        }
+
+        front.transfer(newFront);
+    }
+}
+
+
+Foam::label Foam::processorColour::cellColour
+(
+    const lduMesh& mesh,
+    labelList& cellColour,
+    labelList& patchToColour
+)
+{
+    // Colour with the (min) coupled (global) patch. Compact to have colour 0
+    // if a single region. Problematic if same patch on multiple disconnected
+    // regions.
+
+    const lduAddressing& addr = mesh.lduAddr();
+    const lduInterfacePtrsList patches = mesh.interfaces();
+
+    const label nCells = addr.size();
+
+    patchToColour.resize_nocopy(patches.size());
+    patchToColour = -1;
+    cellColour.resize_nocopy(nCells);
+    cellColour = -1;
+
+
+    label colouri = 0;
+
+    // Starting front from patch faceCells
+    DynamicList<label> front;
+    forAll(patches, inti)
+    {
+        if
+        (
+            patches.set(inti)
+        && !isA<const processorLduInterface>(patches[inti])
+        )
+        {
+            // 'global' interface. Seed faceCells with patch index
+
+            if (patchToColour[inti] == -1)
+            {
+                patchToColour[inti] = colouri++;
+            }
+
+            const auto& fc = patches[inti].faceCells();
+            forAll(fc, face)
+            {
+                const label cell = fc[face];
+                if (cellColour[cell] != patchToColour[inti])
+                {
+                    cellColour[cell] = patchToColour[inti];
+                    front.append(cell);
+                }
+            }
+        }
+    }
+
+
+    // Walk current mesh
+    walkFront(mesh, front, cellColour);
+
+
+    // Any unvisited cell is still labelMax. Put into colour 0.
+    for (auto& colour : cellColour)
+    {
+        if (colour == -1)
+        {
+            colour = 0;
+        }
+    }
+
+    if (debug)
+    {
+        Info<< typeName << " : coloured " << nCells
+            << " cells with in total " << colouri << " colours" << endl;
+    }
+
+    return colouri;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::processorColour::processorColour(const lduMesh& mesh)
+:
+    MeshObject<lduMesh, Foam::MoveableMeshObject, processorColour>(mesh)
+{
+    nColours_ = colour(mesh, *this);
+}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+const Foam::processorColour& Foam::processorColour::New(const lduMesh& mesh)
+{
+    const processorColour* ptr =
+        mesh.thisDb().objectRegistry::template cfindObject<processorColour>
+        (
+            processorColour::typeName
+        );
+
+    if (ptr)
+    {
+        return *ptr;
+    }
+
+    processorColour* objectPtr = new processorColour(mesh);
+
+    //regIOobject::store(static_cast<MoveableMeshObject<lduMesh>*>(objectPtr));
+    regIOobject::store(objectPtr);
+
+    return *objectPtr;
+}
+
+
+// ************************************************************************* //

src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/processorColour.H

@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 M. Janssens
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::processorColour
+
+Description
+    Colouring processors such that no neighbours have the same colour
+
+SourceFiles
+    processorColour.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_processorColour_H
+#define Foam_processorColour_H
+
+#include "MeshObject.H"
+#include "lduMesh.H"
+#include "lduInterfacePtrsList.H"
+#include "primitiveFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class lduMatrix;
+
+/*---------------------------------------------------------------------------*\
+                      Class processorColour Declaration
+\*---------------------------------------------------------------------------*/
+
+class processorColour
+:
+    public MeshObject<lduMesh, MoveableMeshObject, processorColour>,
+    public labelList
+{
+protected:
+
+    // Protected data
+
+        //- Max number of colours
+        label nColours_;
+
+
+    // Protected Member Fucntions
+
+        static void walkFront
+        (
+            const lduMesh& mesh,
+            DynamicList<label>& front,
+            labelList& cellColour
+        );
+
+
+private:
+
+        //- No copy construct
+        processorColour(const processorColour&) = delete;
+
+        //- No copy assignment
+        void operator=(const processorColour&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("processorColour");
+
+
+    // Constructors
+
+        //- Construct given mesh
+        processorColour(const lduMesh& mesh);
+
+
+    // Selectors
+
+        //- Should use the MeshObject provided one but that needs a
+        //  mesh.name() ...
+        static const processorColour& New(const lduMesh& mesh);
+
+
+    //- Destructor
+    virtual ~processorColour() = default;
+
+
+    // Member Functions
+
+        // Access
+
+            label nColours() const
+            {
+                return nColours_;
+            }
+
+            label myColour() const
+            {
+                return operator[](Pstream::myProcNo());
+            }
+
+
+        //- Calculate processor colouring from processor connectivity. Sets
+        //- colour per processor such that no neighbouring processors have the
+        //- same colour. Returns number of colours used.
+        static label colour(const lduMesh& mesh, labelList& procColour);
+
+        //- Calculate (locally) per cell colour according to walking from
+        //- global patches. Returns number of colours used.
+        //- Note: asssumes all non-processor interfaces are global.
+        static label cellColour
+        (
+            const lduMesh& mesh,
+            labelList& cellColour,
+            labelList& patchToColour
+        );
+
+        virtual bool movePoints()
+        {
+            return true;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //