diff --git a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H index 9f751ecd200ba0bbcdaab8a033fc6cffe94316cf..1df6e108f011f4227f2e642687f718d7b0051dd9 100644 --- a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H +++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H @@ -4,7 +4,7 @@ // Maximum length for dynamicList const label maxDynListLength ( - viewFactorDict.getOrDefault<label>("maxDynListLength", 100000) + viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000) ); for (const int proci : Pstream::allProcs()) @@ -47,7 +47,7 @@ for (const int proci : Pstream::allProcs()) const vector& remA = remoteArea[j]; const label& remAgg = remoteAgg[j]; - const vector& d = remFc - fc; + const vector d(remFc - fc); if (((d & fA) < 0.) && ((d & remA) > 0)) { diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C index a32b878ec90c095f9a4fd70963d28df714c49011..fafdfb811b9d3cf9e757efefdf897c8bacf9fbee 100644 --- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C +++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C @@ -408,7 +408,7 @@ int main(int argc, char *argv[]) List<point> availablePoints ( upp.faceCentres().size() - + upp.localPoints().size() + + upp.localPoints().size() ); SubList<point> @@ -863,21 +863,6 @@ int main(int argc, char *argv[]) IOglobalFaceFaces.write(); - labelListIOList IOvisibleFaceFaces - ( - IOobject - ( - "visibleFaceFaces", - mesh.facesInstance(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - std::move(visibleFaceFaces) - ); - IOvisibleFaceFaces.write(); - IOmapDistribute IOmapDist ( IOobject diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.C b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.C new file mode 100644 index 0000000000000000000000000000000000000000..1108e24a4c0b1b9e9b4760c13066222fbd786ff6 --- /dev/null +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.C @@ -0,0 +1,236 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 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 "lduCalculatedProcessorField.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template<class Type> +Foam::lduCalculatedProcessorField<Type>::lduCalculatedProcessorField +( + const lduInterface& interface, + const Field<Type>& iF +) +: + LduInterfaceField<Type>(interface), + procInterface_(refCast<const lduPrimitiveProcessorInterface>(interface)), + field_(iF), + sendBuf_(procInterface_.faceCells().size()), + receiveBuf_(procInterface_.faceCells().size()), + scalarSendBuf_(procInterface_.faceCells().size()), + scalarReceiveBuf_(procInterface_.faceCells().size()), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1) +{} + + +template<class Type> +Foam::lduCalculatedProcessorField<Type>::lduCalculatedProcessorField +( + const lduCalculatedProcessorField<Type>& ptf +) +: + LduInterfaceField<Type>(refCast<const lduInterface>(ptf)), + procInterface_(ptf.procInterface_), + field_(ptf.field_), + sendBuf_(procInterface_.faceCells().size()), + receiveBuf_(procInterface_.faceCells().size()), + scalarSendBuf_(procInterface_.faceCells().size()), + scalarReceiveBuf_(procInterface_.faceCells().size()), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class Type> +bool Foam::lduCalculatedProcessorField<Type>::ready() const +{ + if + ( + this->outstandingSendRequest_ >= 0 + && this->outstandingSendRequest_ < Pstream::nRequests() + ) + { + if (!UPstream::finishedRequest(this->outstandingSendRequest_)) + { + return false; + } + } + this->outstandingSendRequest_ = -1; + + if + ( + this->outstandingRecvRequest_ >= 0 + && this->outstandingRecvRequest_ < Pstream::nRequests() + ) + { + if (!UPstream::finishedRequest(this->outstandingRecvRequest_)) + { + return false; + } + } + this->outstandingRecvRequest_ = -1; + + return true; +} + + +template<class Type> +void Foam::lduCalculatedProcessorField<Type>::initInterfaceMatrixUpdate +( + solveScalarField& result, + const bool add, + const lduAddressing& lduAddr, + const label patchId, + const solveScalarField& psiInternal, + const scalarField& coeffs, + const direction cmpt, + const Pstream::commsTypes commsType +) const +{ + // Bypass patchInternalField since uses fvPatch addressing + const labelList& fc = lduAddr.patchAddr(patchId); + + scalarSendBuf_.setSize(fc.size()); + forAll(fc, i) + { + scalarSendBuf_[i] = psiInternal[fc[i]]; + } + + if (!this->ready()) + { + FatalErrorInFunction + << "On patch " + << " outstanding request." + << abort(FatalError); + } + + + scalarReceiveBuf_.setSize(scalarSendBuf_.size()); + outstandingRecvRequest_ = UPstream::nRequests(); + + UIPstream::read + ( + Pstream::commsTypes::nonBlocking, + procInterface_.neighbProcNo(), + scalarReceiveBuf_.data_bytes(), + scalarReceiveBuf_.size_bytes(), + procInterface_.tag(), + procInterface_.comm() + ); + + outstandingSendRequest_ = UPstream::nRequests(); + + UOPstream::write + ( + Pstream::commsTypes::nonBlocking, + procInterface_.neighbProcNo(), + scalarSendBuf_.cdata_bytes(), + scalarSendBuf_.size_bytes(), + procInterface_.tag(), + procInterface_.comm() + ); + + const_cast<lduInterfaceField&> + ( + static_cast<const lduInterfaceField&>(*this) + ).updatedMatrix() = false; +} + + +template<class Type> +void Foam::lduCalculatedProcessorField<Type>::addToInternalField +( + solveScalarField& result, + const bool add, + const scalarField& coeffs, + const solveScalarField& vals +) const +{ + const labelUList& faceCells = this->procInterface_.faceCells(); + + if (add) + { + forAll(faceCells, elemI) + { + result[faceCells[elemI]] += coeffs[elemI]*vals[elemI]; + } + } + else + { + forAll(faceCells, elemI) + { + result[faceCells[elemI]] -= coeffs[elemI]*vals[elemI]; + } + } +} + + +template<class Type> +void Foam::lduCalculatedProcessorField<Type>::updateInterfaceMatrix +( + solveScalarField& result, + const bool add, + const lduAddressing& lduAddr, + const label patchId, + const solveScalarField& psiInternal, + const scalarField& coeffs, + const direction cmpt, + const Pstream::commsTypes commsType +) const +{ + if (this->updatedMatrix()) + { + return; + } + + if + ( + outstandingRecvRequest_ >= 0 + && outstandingRecvRequest_ < Pstream::nRequests() + ) + { + UPstream::waitRequest(outstandingRecvRequest_); + } + // Recv finished so assume sending finished as well. + outstandingSendRequest_ = -1; + outstandingRecvRequest_ = -1; + + // Consume straight from scalarReceiveBuf_. Note use of our own + // helper to avoid using fvPatch addressing + addToInternalField(result, !add, coeffs, scalarReceiveBuf_); + + const_cast<lduInterfaceField&> + ( + static_cast<const lduInterfaceField&>(*this) + ).updatedMatrix() = true; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.H b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.H new file mode 100644 index 0000000000000000000000000000000000000000..f9a36b1bfb898a389d409e65d69704cbef7fcef9 --- /dev/null +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.H @@ -0,0 +1,254 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 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::lduCalculatedProcessorField + +Group + grpGenericBoundaryConditions + +Description + A lduProcessorField type bypassing coupledFvPatchField and holding + a reference to the Field<Type>. + + Used to add updateInterfaceMatrix capabilities to a lduMatrix + which is fully uncoupled from the fvMesh. + + Its functionality is purely to init and update the processor interfaces. + +SourceFiles + lduCalculatedProcessorField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef lduCalculatedProcessorField_H +#define lduCalculatedProcessorField_H + +#include "lduPrimitiveProcessorInterface.H" +#include "processorLduInterfaceField.H" +#include "LduInterfaceField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class lduCalculatedProcessorField Declaration +\*---------------------------------------------------------------------------*/ + +template<class Type> +class lduCalculatedProcessorField +: + public LduInterfaceField<Type>, + public processorLduInterfaceField +{ +protected: + + // Protected Data + + //- Local reference cast into the interface + const lduPrimitiveProcessorInterface& procInterface_; + + //- Local Field + const Field<Type>& field_; + + // Sending and receiving + + //- Send buffer + mutable Field<Type> sendBuf_; + + //- Receive buffer + mutable Field<Type> receiveBuf_; + + //- Scalar send buffer + mutable solveScalarField scalarSendBuf_; + + //- Scalar receive buffer + mutable solveScalarField scalarReceiveBuf_; + + //- Outstanding request + mutable label outstandingSendRequest_; + + //- Outstanding request + mutable label outstandingRecvRequest_; + + + // Protected Member Functions + + void addToInternalField + ( + solveScalarField& result, + const bool add, + const scalarField& coeffs, + const solveScalarField& vals + ) const; + + +public: + + //- Runtime type information + ClassName("lduCalculatedProcessorField"); + + + // Constructors + + //- Construct from patch and internal field + lduCalculatedProcessorField + ( + const lduInterface& interface, + const Field<Type>& + ); + + //- Construct as copy + lduCalculatedProcessorField + ( + const lduCalculatedProcessorField<Type>& + ); + + + //- Destructor + virtual ~lduCalculatedProcessorField() = default; + + + // Member Functions + + // Access + + //- Return communicator used for comms + virtual label comm() const + { + return procInterface_.comm(); + } + + //- Return processor number + virtual int myProcNo() const + { + return procInterface_.myProcNo(); + } + + //- Return neighbour processor number + virtual int neighbProcNo() const + { + return procInterface_.myProcNo(); + } + + //- Is the transform required + virtual bool doTransform() const + { + return false; + } + + //- Return face transformation tensor + virtual const tensorField& forwardT() const + { + return procInterface_.forwardT(); + } + + //- Return rank of component for transform + virtual int rank() const + { + return pTraits<Type>::rank; + } + + + // Evaluation + + //- Is all data available + virtual bool ready() const; + + //- Initialise neighbour matrix update + virtual void initInterfaceMatrixUpdate + ( + solveScalarField& result, + const bool add, + const lduAddressing& lduAddr, + const label patchId, + const solveScalarField& psiInternal, + const scalarField& coeffs, + const direction cmpt, + const Pstream::commsTypes commsType + ) const; + + //- Update result field based on interface functionality + virtual void updateInterfaceMatrix + ( + solveScalarField& result, + const bool add, + const lduAddressing& lduAddr, + const label patchId, + const solveScalarField& psiInternal, + const scalarField& coeffs, + const direction cmpt, + const Pstream::commsTypes commsType + ) const; + + //- Initialise neighbour matrix update + virtual void initInterfaceMatrixUpdate + ( + Field<scalar>& result, + const bool add, + const lduAddressing& lduAddr, + const label patchId, + const Field<scalar>& psiInternal, + const scalarField& coeffs, + const Pstream::commsTypes commsType + ) const + { + NotImplemented; + } + + //- Update result field based on interface functionality + virtual void updateInterfaceMatrix + ( + Field<scalar>& result, + const bool add, + const lduAddressing& lduAddr, + const label patchId, + const Field<scalar>& psiInternal, + const scalarField& coeffs, + const Pstream::commsTypes commsType + ) const + { + NotImplemented; + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "lduCalculatedProcessorField.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.C b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.C index 2df7eeec308264344e2f3ca114a5b496766e042a..48074e2e3ea64e69e924a9acfd707d20d28b1c78 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2016-2019 OpenCFD Ltd. + Copyright (C) 2016-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -38,6 +38,20 @@ namespace Foam // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // +Foam::lduPrimitiveProcessorInterface::lduPrimitiveProcessorInterface +( + const lduPrimitiveProcessorInterface& p +) +: + faceCells_(p.faceCells()), + myProcNo_(p.myProcNo()), + neighbProcNo_(p.neighbProcNo()), + forwardT_(p.forwardT()), + tag_(p.tag()), + comm_(p.comm()) +{} + + Foam::lduPrimitiveProcessorInterface::lduPrimitiveProcessorInterface ( const labelUList& faceCells, diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.H b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.H index 47217dc6a2ed4ce3ba9e9be816b7fc128e4261f9..f13a2b7f96eded6a40085e4784e2493f1389206c 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.H +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2016 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -78,11 +78,6 @@ class lduPrimitiveProcessorInterface // Private Member Functions - //- No copy construct - lduPrimitiveProcessorInterface - ( - const lduPrimitiveProcessorInterface& - ) = delete; //- No copy assignment void operator=(const lduPrimitiveProcessorInterface&) = delete; @@ -107,6 +102,11 @@ public: const label comm = UPstream::worldComm ); + //- Copy constructor + lduPrimitiveProcessorInterface + ( + const lduPrimitiveProcessorInterface& + ); //- Destructor virtual ~lduPrimitiveProcessorInterface() = default; diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C index 7042b98f6246d29eac052c05d6743ab30bdd85c7..cacc91b306c7b82a412083b8c4f63d2f878cd869 100644 --- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C +++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C @@ -33,6 +33,7 @@ License #include "typeInfo.H" #include "addToRunTimeSelectionTable.H" #include "boundaryRadiationProperties.H" +#include "lduCalculatedProcessorField.H" using namespace Foam::constant; @@ -77,6 +78,8 @@ void Foam::radiation::viewFactor::initialise() DebugInFunction << "Total number of clusters : " << totalNCoarseFaces_ << endl; + useDirect_ = coeffs_.getOrDefault<bool>("useDirectSolver", true); + map_.reset ( new IOmapDistribute @@ -93,94 +96,385 @@ void Foam::radiation::viewFactor::initialise() ) ); - scalarListIOList FmyProc + FmyProc_.reset ( - IOobject + new scalarListIOList ( - "F", - mesh_.facesInstance(), - mesh_, - IOobject::MUST_READ, - IOobject::NO_WRITE, - false + IOobject + ( + "F", + mesh_.facesInstance(), + mesh_, + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ) ) ); - labelListIOList globalFaceFaces + globalFaceFaces_.reset ( - IOobject + new labelListIOList ( - "globalFaceFaces", - mesh_.facesInstance(), - mesh_, - IOobject::MUST_READ, - IOobject::NO_WRITE, - false + IOobject + ( + "globalFaceFaces", + mesh_.facesInstance(), + mesh_, + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ) ) ); - List<labelListList> globalFaceFacesProc(Pstream::nProcs()); - globalFaceFacesProc[Pstream::myProcNo()] = globalFaceFaces; - Pstream::gatherList(globalFaceFacesProc); - - List<scalarListList> F(Pstream::nProcs()); - F[Pstream::myProcNo()] = FmyProc; - Pstream::gatherList(F); globalIndex globalNumbering(nLocalCoarseFaces_); - if (Pstream::master()) + // Size coarse qr + qrBandI_.setSize(nBands_); + forAll (qrBandI_, bandI) + { + qrBandI_[bandI].setSize(nLocalCoarseFaces_, 0.0); + } + + if (!useDirect_) { - Fmatrix_.reset + DynamicList<label> dfaceJ; + + // Per processor to owner (local)/neighbour (remote) + List<DynamicList<label>> procOwner(Pstream::nProcs()); + List<DynamicList<label>> dynProcNeighbour(Pstream::nProcs()); + + forAll (globalFaceFaces_(), iFace) + { + const labelList& visFaces = globalFaceFaces_()[iFace]; + forAll (visFaces, j) + { + label gFacej = visFaces[j]; + label proci = globalNumbering.whichProcID(gFacej); + label faceJ = globalNumbering.toLocal(proci, gFacej); + + if (Pstream::myProcNo() == proci) + { + edge e(iFace, faceJ); + if (rays_.insert(e)) + { + dfaceJ.append(j); + } + } + else + { + label gFaceI = + globalNumbering.toGlobal(Pstream::myProcNo(), iFace); + + label proci = globalNumbering.whichProcID(gFacej); + + label facei = + globalNumbering.toLocal(Pstream::myProcNo(), gFaceI); + + label remoteFacei = globalNumbering.toLocal(proci, gFacej); + + procOwner[proci].append(facei); + dynProcNeighbour[proci].append(remoteFacei); + } + } + } + + mapRayToFmy_.transfer(dfaceJ); + + labelList upper(rays_.size(), -1); + labelList lower(rays_.size(), -1); + + const edgeList& raysLst = rays_.sortedToc(); + label rayI = 0; + for (const auto& e : raysLst) + { + label faceI = e.start(); + label faceJ = e.end(); + upper[rayI] = max(faceI, faceJ); + lower[rayI] = min(faceI, faceJ); + rayI++; + } + + labelListList procNeighbour(dynProcNeighbour.size()); + forAll(procNeighbour, i) + { + procNeighbour[i] = std::move(dynProcNeighbour[i]); + } + labelListList mySendCells; + Pstream::exchange<labelList, label>(procNeighbour, mySendCells); + + label nbri = 0; + forAll(procOwner, proci) + { + if (procOwner[proci].size()) + { + nbri++; + } + if (mySendCells[proci].size()) + { + nbri++; + } + } + + DebugInFunction<< "Number of procBound : " << nbri << endl; + + PtrList<const lduPrimitiveProcessorInterface> primitiveInterfaces; + + primitiveInterfaces.setSize(nbri); + internalCoeffs_.setSize(0); + boundaryCoeffs_.setSize(nbri); + + nbri = 0; + + procToInterface_.setSize(Pstream::nProcs(), -1); + + forAll(procOwner, proci) + { + if (proci < Pstream::myProcNo() && procOwner[proci].size()) + { + if (debug) + { + Pout<< "Adding interface " << nbri + << " to receive my " << procOwner[proci].size() + << " from " << proci << endl; + } + procToInterface_[proci] = nbri; + primitiveInterfaces.set + ( + nbri++, + new lduPrimitiveProcessorInterface + ( + procOwner[proci], + Pstream::myProcNo(), + proci, + tensorField(0), + Pstream::msgType()+2 + ) + ); + } + else if (proci > Pstream::myProcNo() && mySendCells[proci].size()) + { + if (debug) + { + Pout<< "Adding interface " << nbri + << " to send my " << mySendCells[proci].size() + << " to " << proci << endl; + } + primitiveInterfaces.set + ( + nbri++, + new lduPrimitiveProcessorInterface + ( + mySendCells[proci], + Pstream::myProcNo(), + proci, + tensorField(0), + Pstream::msgType()+2 + ) + ); + } + } + forAll(procOwner, proci) + { + if (proci > Pstream::myProcNo() && procOwner[proci].size()) + { + if (debug) + { + Pout<< "Adding interface " << nbri + << " to receive my " << procOwner[proci].size() + << " from " << proci << endl; + } + procToInterface_[proci] = nbri; + primitiveInterfaces.set + ( + nbri++, + new lduPrimitiveProcessorInterface + ( + procOwner[proci], + Pstream::myProcNo(), + proci, + tensorField(0), + Pstream::msgType()+3 + ) + ); + } + else if (proci < Pstream::myProcNo() && mySendCells[proci].size()) + { + if (debug) + { + Pout<< "Adding interface " << nbri + << " to send my " << mySendCells[proci].size() + << " to " << proci << endl; + } + primitiveInterfaces.set + ( + nbri++, + new lduPrimitiveProcessorInterface + ( + mySendCells[proci], + Pstream::myProcNo(), + proci, + tensorField(0), + Pstream::msgType()+3 + ) + ); + } + } + + forAll (boundaryCoeffs_, proci) + { + boundaryCoeffs_.set + ( + proci, + new Field<scalar> + ( + primitiveInterfaces[proci].faceCells().size(), + Zero + ) + ); + } + + lduInterfacePtrsList allInterfaces; + allInterfaces.setSize(primitiveInterfaces.size()); + + forAll(primitiveInterfaces, i) + { + const lduPrimitiveProcessorInterface& pp = primitiveInterfaces[i]; + + allInterfaces.set(i, &pp); + } + + const lduSchedule ps ( - new scalarSquareMatrix(totalNCoarseFaces_, Zero) + lduPrimitiveMesh::nonBlockingSchedule + <lduPrimitiveProcessorInterface>(allInterfaces) ); - DebugInFunction - << "Insert elements in the matrix..." << endl; - - for (const int procI : Pstream::allProcs()) + PtrList<const lduInterface> allInterfacesPtr(allInterfaces.size()); + forAll (allInterfacesPtr, i) { - insertMatrixElements + const lduPrimitiveProcessorInterface& pp = primitiveInterfaces[i]; + + allInterfacesPtr.set ( - globalNumbering, - procI, - globalFaceFacesProc[procI], - F[procI], - Fmatrix_() + i, + new lduPrimitiveProcessorInterface(pp) ); } + lduPtr_.reset + ( + new lduPrimitiveMesh + ( + rays_.size(), + lower, + upper, + allInterfacesPtr, + ps, + UPstream::worldComm + ) + ); + + // Set size for local lduMatrix + matrixPtr_.reset(new lduMatrix(lduPtr_())); + + scalarListList& myF = FmyProc_(); if (coeffs_.get<bool>("smoothing")) { - DebugInFunction << "Smoothing the matrix..." << endl; - - for (label i=0; i<totalNCoarseFaces_; i++) + scalar maxDelta = 0; + scalar totalDelta = 0; + forAll (myF, i) { scalar sumF = 0.0; - for (label j=0; j<totalNCoarseFaces_; j++) + scalarList& myFij = myF[i]; + forAll (myFij, j) { - sumF += Fmatrix_()(i, j); + sumF += myFij[j]; } - const scalar delta = sumF - 1.0; - for (label j=0; j<totalNCoarseFaces_; j++) + forAll (myFij, j) { - Fmatrix_()(i, j) *= (1.0 - delta/(sumF + 0.001)); + myFij[j] *= (1.0 - delta/(sumF + 0.001)); + } + totalDelta += delta; + if (delta > maxDelta) + { + maxDelta = delta; } } + totalDelta /= myF.size(); + reduce(totalDelta, sumOp<scalar>()); + reduce(maxDelta, maxOp<scalar>()); + Info << "Smoothng average delta : " << totalDelta << endl; + Info << "Smoothng maximum delta : " << maxDelta << nl << endl; } + } + + if (useDirect_) + { + List<labelListList> globalFaceFacesProc(Pstream::nProcs()); + globalFaceFacesProc[Pstream::myProcNo()] = globalFaceFaces_(); + Pstream::gatherList(globalFaceFacesProc); + + List<scalarListList> F(Pstream::nProcs()); + F[Pstream::myProcNo()] = FmyProc_(); + Pstream::gatherList(F); - coeffs_.readEntry("constantEmissivity", constEmissivity_); - if (constEmissivity_) + if (Pstream::master()) { - CLU_.reset + Fmatrix_.reset ( new scalarSquareMatrix(totalNCoarseFaces_, Zero) ); - pivotIndices_.setSize(CLU_().m()); + DebugInFunction + << "Insert elements in the matrix..." << endl; + + for (const int procI : Pstream::allProcs()) + { + insertMatrixElements + ( + globalNumbering, + procI, + globalFaceFacesProc[procI], + F[procI], + Fmatrix_() + ); + } + + if (coeffs_.get<bool>("smoothing")) + { + DebugInFunction << "Smoothing the matrix..." << endl; + + for (label i=0; i<totalNCoarseFaces_; i++) + { + scalar sumF = 0.0; + for (label j=0; j<totalNCoarseFaces_; j++) + { + sumF += Fmatrix_()(i, j); + } + + const scalar delta = sumF - 1.0; + for (label j=0; j<totalNCoarseFaces_; j++) + { + Fmatrix_()(i, j) *= (1.0 - delta/(sumF + 0.001)); + } + } + } + + coeffs_.readEntry("constantEmissivity", constEmissivity_); + if (constEmissivity_) + { + CLU_.reset + ( + new scalarSquareMatrix(totalNCoarseFaces_, Zero) + ); + + pivotIndices_.setSize(CLU_().m()); + } } } @@ -258,7 +552,8 @@ Foam::radiation::viewFactor::viewFactor(const volScalarField& T) pivotIndices_(0), useSolarLoad_(false), solarLoad_(), - nBands_(coeffs_.getOrDefault<label>("nBands", 1)) + nBands_(coeffs_.getOrDefault<label>("nBands", 1)), + FmyProc_() { initialise(); } @@ -320,7 +615,8 @@ Foam::radiation::viewFactor::viewFactor pivotIndices_(0), useSolarLoad_(false), solarLoad_(), - nBands_(coeffs_.getOrDefault<label>("nBands", 1)) + nBands_(coeffs_.getOrDefault<label>("nBands", 1)), + FmyProc_() { initialise(); } @@ -372,8 +668,13 @@ void Foam::radiation::viewFactor::calculate() solarLoad_->calculate(); } - // Net radiation - scalarField q(totalNCoarseFaces_, 0.0); + // Global net radiation + scalarField qNet(totalNCoarseFaces_, 0.0); + + // Local net radiation + scalarField qTotalCoarse(nLocalCoarseFaces_, 0.0); + + // Referen to fvMesh qr field volScalarField::Boundary& qrBf = qr_.boundaryFieldRef(); globalIndex globalNumbering(nLocalCoarseFaces_); @@ -383,6 +684,9 @@ void Foam::radiation::viewFactor::calculate() for (label bandI = 0; bandI < nBands_; bandI++) { + // Global bandI radiation + scalarField qBandI(totalNCoarseFaces_, 0.0); + scalarField compactCoarseT4(map_->constructSize(), 0.0); scalarField compactCoarseE(map_->constructSize(), 0.0); scalarField compactCoarseHo(map_->constructSize(), 0.0); @@ -467,6 +771,7 @@ void Foam::radiation::viewFactor::calculate() SubList<scalar>(compactCoarseHo, nLocalCoarseFaces_) = localCoarseHoave; + // Distribute data map_->distribute(compactCoarseT4); map_->distribute(compactCoarseE); @@ -487,123 +792,260 @@ void Foam::radiation::viewFactor::calculate() map_->distribute(compactGlobalIds); - // Create global size vectors - scalarField T4(totalNCoarseFaces_, 0.0); - scalarField E(totalNCoarseFaces_, 0.0); - scalarField qrExt(totalNCoarseFaces_, 0.0); - - // Fill lists from compact to global indexes. - forAll(compactCoarseT4, i) + if (!useDirect_) { - T4[compactGlobalIds[i]] = compactCoarseT4[i]; - E[compactGlobalIds[i]] = compactCoarseE[i]; - qrExt[compactGlobalIds[i]] = compactCoarseHo[i]; - } + const labelList globalToCompact + ( + invert(totalNCoarseFaces_, compactGlobalIds) + ); - Pstream::listCombineAllGather(T4, maxEqOp<scalar>()); - Pstream::listCombineAllGather(E, maxEqOp<scalar>()); - Pstream::listCombineAllGather(qrExt, maxEqOp<scalar>()); + scalarField& diag = matrixPtr_->diag(localCoarseEave.size()); + scalarField& upper = matrixPtr_->upper(rays_.size()); + scalarField& lower = matrixPtr_->lower(rays_.size()); - if (Pstream::master()) - { - // Variable emissivity - if (!constEmissivity_) + scalarField source(nLocalCoarseFaces_, 0); + + + // Local diag and source + forAll(source, i) { - scalarSquareMatrix C(totalNCoarseFaces_, 0.0); + const scalar sigmaT4 = + physicoChemical::sigma.value()*localCoarseT4ave[i]; - for (label i=0; i<totalNCoarseFaces_; i++) + diag[i] = 1/localCoarseEave[i]; + + source[i] += -sigmaT4 + localCoarseHoave[i]; + } + + // Local matrix coefficients + if (!constEmissivity_ || iterCounter_ == 0) + { + const edgeList& raysLst = rays_.sortedToc(); + + label rayI = 0; + for (const auto& e : raysLst) { - for (label j=0; j<totalNCoarseFaces_; j++) + label facelJ = e.end(); + label faceI = e.start(); + + label faceJ = mapRayToFmy_[rayI]; + + lower[rayI] = + (1-1/localCoarseEave[faceI])*FmyProc_()[faceI][faceJ]; + + upper[rayI] = + (1-1/localCoarseEave[facelJ])*FmyProc_()[faceI][faceJ]; + + rayI++; + } + iterCounter_++; + } + + // Extra local contribution to the source and extra matrix coefs + // from other procs (boundaryCoeffs) + label nInterfaces = lduPtr_().interfaces().size(); + labelList boundCoeffI(nInterfaces, Zero); + forAll (globalFaceFaces_(), iFace) + { + const labelList& visFaces = globalFaceFaces_()[iFace]; + forAll (visFaces, jFace) + { + label gFacej = visFaces[jFace]; + label proci = globalNumbering.whichProcID(gFacej); + if (Pstream::myProcNo() == proci) { - const scalar invEj = 1.0/E[j]; + label lFacej = + globalNumbering.toLocal + ( + Pstream::myProcNo(), + gFacej + ); + const scalar sigmaT4 = - physicoChemical::sigma.value()*T4[j]; + physicoChemical::sigma.value() + *localCoarseT4ave[lFacej]; - if (i==j) - { - C(i, j) = invEj - (invEj - 1.0)*Fmatrix_()(i, j); - q[i] += - (Fmatrix_()(i, j) - 1.0)*sigmaT4 + qrExt[j]; - } - else - { - C(i, j) = (1.0 - invEj)*Fmatrix_()(i, j); - q[i] += Fmatrix_()(i, j)*sigmaT4; - } + source[iFace] += FmyProc_()[iFace][jFace]*sigmaT4; + } + else + { + label compactFaceJ = globalToCompact[gFacej]; + const scalar sigmaT4 = + physicoChemical::sigma.value() + *compactCoarseT4[compactFaceJ]; + + source[iFace] += FmyProc_()[iFace][jFace]*sigmaT4; + + label interfaceI = procToInterface_[proci]; + boundaryCoeffs_ + [interfaceI][boundCoeffI[interfaceI]++] = + -(1-1/compactCoarseE[compactFaceJ]) + *FmyProc_()[iFace][jFace]; } } + } - Info<< "Solving view factor equations for band :" - << bandI << endl; + PtrList<const lduInterfaceField> interfaces(nInterfaces); + for(label i = 0; i < interfaces.size(); i++) + { + interfaces.set + ( + i, + new lduCalculatedProcessorField<scalar> + ( + lduPtr_().interfaces()[i], + qrBandI_[bandI] + ) + ); + } - // Negative coming into the fluid - LUsolve(C, q); + const dictionary& solverControls = + qr_.mesh().solverDict + ( + qr_.select + ( + qr_.mesh().data::template getOrDefault<bool> + ("finalIteration", false) + ) + ); + + // Solver call + solverPerformance solverPerf = lduMatrix::solver::New + ( + "qr", + matrixPtr_(), + boundaryCoeffs_, + internalCoeffs_, + interfaces, + solverControls + )->solve(qrBandI_[bandI], source); + + solverPerf.print(Info.masterStream(qr_.mesh().comm())); + + qTotalCoarse += qrBandI_[bandI]; + } + + if (useDirect_) + { + // Create global size vectors + scalarField T4(totalNCoarseFaces_, 0.0); + scalarField E(totalNCoarseFaces_, 0.0); + scalarField qrExt(totalNCoarseFaces_, 0.0); + + // Fill lists from compact to global indexes. + forAll(compactCoarseT4, i) + { + T4[compactGlobalIds[i]] = compactCoarseT4[i]; + E[compactGlobalIds[i]] = compactCoarseE[i]; + qrExt[compactGlobalIds[i]] = compactCoarseHo[i]; } - else //Constant emissivity + + Pstream::listCombineGather(T4, maxEqOp<scalar>()); + Pstream::listCombineGather(E, maxEqOp<scalar>()); + Pstream::listCombineGather(qrExt, maxEqOp<scalar>()); + + Pstream::listCombineScatter(T4); + Pstream::listCombineScatter(E); + Pstream::listCombineScatter(qrExt); + + if (Pstream::master()) { - // Initial iter calculates CLU and caches it - if (iterCounter_ == 0) + // Variable emissivity + if (!constEmissivity_) { + scalarSquareMatrix C(totalNCoarseFaces_, 0.0); + for (label i=0; i<totalNCoarseFaces_; i++) { for (label j=0; j<totalNCoarseFaces_; j++) { const scalar invEj = 1.0/E[j]; + const scalar sigmaT4 = + physicoChemical::sigma.value()*T4[j]; + if (i==j) { - CLU_()(i, j) = - invEj-(invEj-1.0)*Fmatrix_()(i, j); + C(i, j) = invEj; + qBandI[i] += -sigmaT4 + qrExt[j]; } else { - CLU_()(i, j) = (1.0 - invEj)*Fmatrix_()(i, j); + C(i, j) = (1.0 - invEj)*Fmatrix_()(i, j); + qBandI[i] += Fmatrix_()(i, j)*sigmaT4; } } } - DebugInFunction << "\nDecomposing C matrix..." << endl; + Info<< "Solving view factor equations for band :" + << bandI << endl; - LUDecompose(CLU_(), pivotIndices_); + // Negative coming into the fluid + LUsolve(C, qBandI); } - - for (label i=0; i<totalNCoarseFaces_; i++) + else //Constant emissivity { - for (label j=0; j<totalNCoarseFaces_; j++) + // Initial iter calculates CLU and caches it + if (iterCounter_ == 0) { - const scalar sigmaT4 = - constant::physicoChemical::sigma.value()*T4[j]; - - if (i==j) - { - q[i] += - (Fmatrix_()(i, j) - 1.0)*sigmaT4 + qrExt[j]; - } - else + for (label i=0; i<totalNCoarseFaces_; i++) { - q[i] += Fmatrix_()(i, j)*sigmaT4; + for (label j=0; j<totalNCoarseFaces_; j++) + { + const scalar invEj = 1.0/E[j]; + if (i==j) + { + CLU_()(i, j) = invEj; + } + else + { + CLU_()(i, j) = + (1.0-invEj)*Fmatrix_()(i, j); + } + } } + + DebugInFunction << "\nDecomposing C matrix..." << endl; + + LUDecompose(CLU_(), pivotIndices_); } - } + for (label i=0; i<totalNCoarseFaces_; i++) + { + for (label j=0; j<totalNCoarseFaces_; j++) + { + const scalar sigmaT4 = + constant::physicoChemical::sigma.value()*T4[j]; - Info<< "Solving view factor equations for band : " - << bandI << endl; + if (i==j) + { + qBandI[i] += -sigmaT4 + qrExt[j]; + } + else + { + qBandI[i] += Fmatrix_()(i, j)*sigmaT4; + } + } + } + Info<< "Solving view factor equations for band : " + << bandI << endl; - LUBacksubstitute(CLU_(), pivotIndices_, q); - iterCounter_ ++; + LUBacksubstitute(CLU_(), pivotIndices_, qBandI); + iterCounter_ ++; + } } - } + // Broadcast qBandI and fill qr + Pstream::broadcast(qBandI); + qNet += qBandI; + } } - // Broadcast q and fill qr - Pstream::broadcast(q); label globCoarseId = 0; for (const label patchID : selectedPatches_) { - const polyPatch& pp = mesh_.boundaryMesh()[patchID]; + const polyPatch& pp = coarseMesh_.boundaryMesh()[patchID]; if (pp.size() > 0) { @@ -617,22 +1059,31 @@ void Foam::radiation::viewFactor::calculate() const labelList& coarsePatchFace = coarseMesh_.patchFaceMap()[patchID]; - /// scalar heatFlux = 0.0; + //scalar heatFlux = 0.0; forAll(coarseToFine, coarseI) { - label globalCoarse = - globalNumbering.toGlobal - (Pstream::myProcNo(), globCoarseId); + label globalCoarse = globalNumbering.toGlobal + ( + Pstream::myProcNo(), + globCoarseId + ); const label coarseFaceID = coarsePatchFace[coarseI]; const labelList& fineFaces = coarseToFine[coarseFaceID]; for (const label facei : fineFaces) { - qrp[facei] = q[globalCoarse]; - /// heatFlux += qrp[facei]*sf[facei]; + if (useDirect_) + { + qrp[facei] = qNet[globalCoarse]; + } + else + { + qrp[facei] = qTotalCoarse[globCoarseId]; + } + //heatFlux += qrp[facei]*sf[facei]; } - globCoarseId ++; + globCoarseId++; } } } @@ -645,8 +1096,7 @@ void Foam::radiation::viewFactor::calculate() const scalarField& magSf = mesh_.magSf().boundaryField()[patchID]; const scalar heatFlux = gSum(qrp*magSf); - InfoInFunction - << "Total heat transfer rate at patch: " + Info<< "Total heat transfer rate at patch: " << patchID << " " << heatFlux << endl; } diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H index 7592cfdf50541ae18aede651218f204c69935d57..662c6586beeedf10bd3537041f4c3a580328833b 100644 --- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H +++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H @@ -60,6 +60,10 @@ SourceFiles #include "IOmapDistribute.H" #include "solarLoad.H" +#include "lduPrimitiveMesh.H" +#include "lduPrimitiveProcessorInterface.H" +#include "lduMatrix.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -108,12 +112,18 @@ protected: //- Net radiative heat flux [W/m2] volScalarField qr_; + //- Coarse radiative heat flux + List<scalarField> qrBandI_; + //- View factor matrix autoPtr<scalarSquareMatrix> Fmatrix_; //- Inverse of C matrix autoPtr<scalarSquareMatrix> CLU_; + //- Visible global faces + autoPtr<labelListIOList> globalFaceFaces_; + //- Selected patches labelList selectedPatches_; @@ -141,6 +151,35 @@ protected: //-Number of bands label nBands_; + //- Primitive addressing for lduMatrix + autoPtr<lduPrimitiveMesh> lduPtr_; + + //- Matrix formed from view factors + autoPtr<lduMatrix> matrixPtr_; + + //- Boundary scalar field containing pseudo-matrix coeffs + //- for internal cells + FieldField<Field, scalar> internalCoeffs_; + + //- Boundary scalar field containing pseudo-matrix coeffs + //- for boundary cells + FieldField<Field, scalar> boundaryCoeffs_; + + //- Rays on local proc + edgeHashSet rays_; + + //- Map local-ray to j-column for F + List<label> mapRayToFmy_; + + //- Local view factors + autoPtr<scalarListIOList> FmyProc_; + + //- Map from proc to interafce + labelList procToInterface_; + + //- Use direct or iterative solver + bool useDirect_; + // Private Member Functions