From c8e224b598bb8162b55fca4f9294b7ad20a1ff38 Mon Sep 17 00:00:00 2001 From: Mattijs Janssens Date: Thu, 1 Jun 2023 16:21:00 +0000 Subject: [PATCH] ENH: lduMatrix: parallel version of DIC,DILU --- src/meshTools/Make/files | 5 + .../distributedDICPreconditioner.C | 135 +++ .../distributedDICPreconditioner.H | 131 +++ .../distributedDILUPreconditioner.C | 897 ++++++++++++++++++ .../distributedDILUPreconditioner.H | 274 ++++++ .../processorColour.C | 363 +++++++ .../processorColour.H | 156 +++ 7 files changed, 1961 insertions(+) create mode 100644 src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDICPreconditioner.C create mode 100644 src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDICPreconditioner.H create mode 100644 src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDILUPreconditioner.C create mode 100644 src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDILUPreconditioner.H create mode 100644 src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/processorColour.C create mode 100644 src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/processorColour.H diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index 4f974d5d2f..bd0112da39 100644 --- a/src/meshTools/Make/files +++ b/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 diff --git a/src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDICPreconditioner.C b/src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDICPreconditioner.C new file mode 100644 index 0000000000..64a0bd6210 --- /dev/null +++ b/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 . + +\*---------------------------------------------------------------------------*/ + +#include "distributedDICPreconditioner.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(distributedDICPreconditioner, 0); + + lduMatrix::preconditioner:: + addsymMatrixConstructorToTable + 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. + +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 + +// ************************************************************************* // diff --git a/src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDILUPreconditioner.C b/src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDILUPreconditioner.C new file mode 100644 index 0000000000..1f92ab7dc4 --- /dev/null +++ b/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 . + +\*---------------------------------------------------------------------------*/ + +#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 + adddistributedDILUPreconditionerAsymMatrixConstructorToTable_; + + const lduMesh* distributedDILUPreconditioner::meshPtr_ = nullptr; + + autoPtr distributedDILUPreconditioner::procColoursPtr_(nullptr); +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::distributedDILUPreconditioner::updateMatrixInterfaces +( + const bool add, + const FieldField& 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 one(interfaces.size()); + FieldField 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 + ( + 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(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& 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(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& requests +) const +{ + const auto& interfaces = solver_.interfaces(); + + // Start writes + for (const label inti : selectedInterfaces) + { + const auto& intf = interfaces[inti].interface(); + const auto* ppp = isA(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& requests +) const +{ + UPstream::waitRequests(requests); + requests.clear(); +} + + +// Diagonal calculation +// ~~~~~~~~~~~~~~~~~~~~ + +void Foam::distributedDILUPreconditioner::addInterfaceDiag +( + solveScalarField& rD, + const label inti, + const Field& 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& 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=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("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(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(intf)) + { + haveGlobalCoupled = true; + + const auto* AMIpp = isA(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(intf); + if (AMIpp) + { + nbrInti = AMIpp->neighbPatchID(); + } + const auto* cycpp = isA(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= 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_); +} + + +// ************************************************************************* // diff --git a/src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDILUPreconditioner.H b/src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDILUPreconditioner.H new file mode 100644 index 0000000000..bda91b01ea --- /dev/null +++ b/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 . + +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 procColoursPtr_; + + //- Buffers for sending and receiving data + mutable FieldField sendBufs_; + mutable FieldField recvBufs_; + mutable DynamicList recvRequests_; + + //- Interfaces to lower coloured processors + DynamicList