64 finalAgglom_.setSize(
patches.size());
65 for (label patchi=0; patchi <
patches.size(); patchi++)
77 "coarse:" + mesh_.name(),
78 mesh_.polyMesh::instance(),
91 selectedPatches_ = mesh_.boundaryMesh().
indices(viewFactorWalls);
93 for (
const label patchi : selectedPatches_)
95 nLocalCoarseFaces_ += coarsePatches[patchi].
size();
100 Pout<<
"radiation::viewFactor::initialise() Selected patches:" 101 << selectedPatches_ <<
endl;
102 Pout<<
"radiation::viewFactor::initialise() Number of coarse faces:" 103 << nLocalCoarseFaces_ <<
endl;
106 totalNCoarseFaces_ =
returnReduce(nLocalCoarseFaces_, sumOp<label>());
109 <<
"Total number of clusters : " << totalNCoarseFaces_ <<
endl;
111 useDirect_ = coeffs_.getOrDefault<
bool>(
"useDirectSolver",
true);
122 mesh_.facesInstance(),
138 mesh_.facesInstance(),
147 globalFaceFaces_.reset
154 mesh_.facesInstance(),
164 globalIndex globalNumbering(nLocalCoarseFaces_);
167 qrBandI_.setSize(nBands_);
170 qrBandI_[bandI].setSize(nLocalCoarseFaces_, 0.0);
175 DynamicList<label> dfaceJ;
181 forAll (globalFaceFaces_(), iFace)
183 const labelList& visFaces = globalFaceFaces_()[iFace];
186 label gFacej = visFaces[j];
188 label faceJ = globalNumbering.
toLocal(proci, gFacej);
192 edge
e(iFace, faceJ);
208 label remoteFacei = globalNumbering.
toLocal(proci, gFacej);
210 procOwner[proci].append(facei);
211 dynProcNeighbour[proci].append(remoteFacei);
216 mapRayToFmy_.transfer(dfaceJ);
221 const edgeList raysLst(rays_.sortedToc());
223 for (
const auto&
e : raysLst)
225 label faceI =
e.start();
226 label faceJ =
e.end();
235 procNeighbour[i] = std::move(dynProcNeighbour[i]);
238 Pstream::exchange<labelList, label>(procNeighbour, mySendCells);
243 if (procOwner[proci].size())
247 if (mySendCells[proci].size())
255 PtrList<const lduPrimitiveProcessorInterface> primitiveInterfaces(nbri);
256 internalCoeffs_.clear();
257 boundaryCoeffs_.resize_null(nbri);
269 Pout<<
"Adding interface " << nbri
270 <<
" to receive my " << procOwner[proci].size()
271 <<
" from " << proci <<
endl;
273 procToInterface_[proci] = nbri;
274 primitiveInterfaces.set
277 new lduPrimitiveProcessorInterface
291 Pout<<
"Adding interface " << nbri
292 <<
" to send my " << mySendCells[proci].size()
293 <<
" to " << proci <<
endl;
295 primitiveInterfaces.set
298 new lduPrimitiveProcessorInterface
315 Pout<<
"Adding interface " << nbri
316 <<
" to receive my " << procOwner[proci].size()
317 <<
" from " << proci <<
endl;
319 procToInterface_[proci] = nbri;
320 primitiveInterfaces.set
323 new lduPrimitiveProcessorInterface
337 Pout<<
"Adding interface " << nbri
338 <<
" to send my " << mySendCells[proci].size()
339 <<
" to " << proci <<
endl;
341 primitiveInterfaces.set
344 new lduPrimitiveProcessorInterface
356 forAll (boundaryCoeffs_, proci)
363 primitiveInterfaces[proci].faceCells().size(),
371 forAll(primitiveInterfaces, i)
373 const lduPrimitiveProcessorInterface&
pp = primitiveInterfaces[i];
375 allInterfaces.set(i, &
pp);
381 <lduPrimitiveProcessorInterface>(allInterfaces)
384 PtrList<const lduInterface> allInterfacesPtr(allInterfaces.size());
385 forAll (allInterfacesPtr, i)
387 const lduPrimitiveProcessorInterface&
pp = primitiveInterfaces[i];
392 new lduPrimitiveProcessorInterface(
pp)
400 globalFaceFaces_().size(),
410 matrixPtr_.reset(
new lduMatrix(lduPtr_()));
414 if (coeffs_.get<
bool>(
"smoothing"))
417 scalar totalDelta = 0;
429 const scalar
delta = sumF - 1.0;
432 myFij[j] *= (1.0 -
delta/(sumF + 0.001));
435 if (
delta > maxDelta)
440 totalDelta /= myF.size();
442 reduce(totalDelta, sumOp<scalar>());
443 reduce(maxDelta, maxOp<scalar>());
444 Info <<
"Smoothing average delta : " << totalDelta <<
endl;
445 Info <<
"Smoothing maximum delta : " << maxDelta <<
nl <<
endl;
467 <<
"Insert elements in the matrix..." <<
endl;
475 globalFaceFacesProc[procI],
481 if (coeffs_.get<
bool>(
"smoothing"))
485 for (label i=0; i<totalNCoarseFaces_; i++)
488 for (label j=0; j<totalNCoarseFaces_; j++)
490 sumF += Fmatrix_()(i, j);
493 const scalar
delta = sumF - 1.0;
494 for (label j=0; j<totalNCoarseFaces_; j++)
496 Fmatrix_()(i, j) *= (1.0 -
delta/(sumF + 0.001));
501 coeffs_.readEntry(
"constantEmissivity", constEmissivity_);
502 if (constEmissivity_)
509 pivotIndices_.setSize(CLU_().m());
514 coeffs_.readIfPresent(
"useSolarLoad", useSolarLoad_);
518 const dictionary& solarDict = this->subDict(
"solarLoadCoeffs");
519 solarLoad_.reset(
new solarLoad(solarDict, T_));
521 if (solarLoad_->nBands() != nBands_)
524 <<
"Solar radiation and view factor band numbers " 529 Info<<
"Creating Solar Load Model " <<
nl;
544 mesh_.facesInstance(),
580 selectedPatches_(mesh_.
boundary().size(), -1),
581 totalNCoarseFaces_(0),
582 nLocalCoarseFaces_(0),
583 constEmissivity_(false),
586 useSolarLoad_(false),
588 nBands_(coeffs_.getOrDefault<label>(
"nBands", 1)),
595 Foam::radiation::viewFactor::viewFactor
607 mesh_.facesInstance(),
643 selectedPatches_(mesh_.
boundary().size(), -1),
644 totalNCoarseFaces_(0),
645 nLocalCoarseFaces_(0),
646 constEmissivity_(false),
649 useSolarLoad_(false),
651 nBands_(coeffs_.getOrDefault<label>(
"nBands", 1)),
680 forAll(viewFactors, faceI)
683 const labelList& globalFaces = globalFaceFaces[faceI];
685 label globalI = globalNumbering.
toGlobal(procI, faceI);
688 Fmatrix[globalI][globalFaces[i]] = vf[i];
701 solarLoad_->calculate();
713 globalIndex globalNumbering(nLocalCoarseFaces_);
715 const boundaryRadiationProperties& boundaryRadiation =
718 for (label bandI = 0; bandI < nBands_; bandI++)
723 scalarField compactCoarseT4(map_->constructSize(), 0.0);
724 scalarField compactCoarseE(map_->constructSize(), 0.0);
725 scalarField compactCoarseHo(map_->constructSize(), 0.0);
728 DynamicList<scalar> localCoarseT4ave(nLocalCoarseFaces_);
729 DynamicList<scalar> localCoarseEave(nLocalCoarseFaces_);
730 DynamicList<scalar> localCoarseHoave(nLocalCoarseFaces_);
732 forAll(selectedPatches_, i)
734 label
patchID = selectedPatches_[i];
744 greyDiffusiveViewFactorFixedValueFvPatchScalarField
747 const tmp<scalarField> teb =
748 boundaryRadiation.emissivity(
patchID, bandI,
nullptr, &Tp);
752 const tmp<scalarField> tHoi = qrp.qro(bandI);
755 const polyPatch&
pp = coarseMesh_->boundaryMesh()[
patchID];
758 coarseMesh_->patchFaceMap()[
patchID];
767 label nAgglom =
max(agglom) + 1;
771 forAll(coarseToFine, coarseI)
773 const label coarseFaceID = coarsePatchFace[coarseI];
774 const labelList& fineFaces = coarseToFine[coarseFaceID];
775 UIndirectList<scalar> fineSf
781 const scalar
area =
sum(fineSf());
786 label facei = fineFaces[j];
787 T4ave[coarseI] += (
pow4(Tp[facei])*sf[facei])/
area;
788 Eave[coarseI] += (eb[facei]*sf[facei])/
area;
789 Hoiave[coarseI] += (Hoi[facei]*sf[facei])/
area;
794 localCoarseT4ave.
append(T4ave);
795 localCoarseEave.append(Eave);
796 localCoarseHoave.append(Hoiave);
801 SubList<scalar>(compactCoarseT4, nLocalCoarseFaces_) =
803 SubList<scalar>(compactCoarseE, nLocalCoarseFaces_) = localCoarseEave;
804 SubList<scalar>(compactCoarseHo, nLocalCoarseFaces_) =
809 map_->distribute(compactCoarseT4);
810 map_->distribute(compactCoarseE);
811 map_->distribute(compactCoarseHo);
826 map_->distribute(compactGlobalIds);
832 invert(totalNCoarseFaces_, compactGlobalIds)
845 const scalar sigmaT4 =
848 diag[i] = 1/localCoarseEave[i];
850 source[i] += -sigmaT4 + localCoarseHoave[i];
854 if (!constEmissivity_ || iterCounter_ == 0)
856 const edgeList raysLst(rays_.sortedToc());
859 for (
const auto&
e : raysLst)
861 label facelJ =
e.end();
862 label faceI =
e.start();
864 label faceJ = mapRayToFmy_[rayI];
867 (1-1/localCoarseEave[faceI])*FmyProc_()[faceI][faceJ];
870 (1-1/localCoarseEave[facelJ])*FmyProc_()[faceI][faceJ];
879 label nInterfaces = lduPtr_().interfaces().size();
881 forAll (globalFaceFaces_(), iFace)
883 const labelList& visFaces = globalFaceFaces_()[iFace];
886 label gFacej = visFaces[jFace];
897 const scalar sigmaT4 =
899 *localCoarseT4ave[lFacej];
901 source[iFace] += FmyProc_()[iFace][jFace]*sigmaT4;
905 label compactFaceJ = globalToCompact[gFacej];
906 const scalar sigmaT4 =
908 *compactCoarseT4[compactFaceJ];
910 source[iFace] += FmyProc_()[iFace][jFace]*sigmaT4;
912 label interfaceI = procToInterface_[proci];
915 [interfaceI][boundCoeffI[interfaceI]++] =
916 -(1-1/compactCoarseE[compactFaceJ])
917 *FmyProc_()[iFace][jFace];
922 PtrList<const lduInterfaceField> interfaces(nInterfaces);
923 for(label i = 0; i < interfaces.size(); i++)
928 new lduCalculatedProcessorField<scalar>
930 lduPtr_().interfaces()[i]
936 qr_.mesh().solverDict
938 qr_.select(qr_.mesh().data().isFinalIteration())
950 )->solve(qrBandI_[bandI], source);
954 qTotalCoarse += qrBandI_[bandI];
965 forAll(compactCoarseT4, i)
967 T4[compactGlobalIds[i]] = compactCoarseT4[i];
968 E[compactGlobalIds[i]] = compactCoarseE[i];
969 qrExt[compactGlobalIds[i]] = compactCoarseHo[i];
979 if (!constEmissivity_)
983 for (label i=0; i<totalNCoarseFaces_; i++)
985 for (label j=0; j<totalNCoarseFaces_; j++)
987 const scalar invEj = 1.0/E[j];
988 const scalar sigmaT4 =
994 qBandI[i] += -sigmaT4 + qrExt[j];
998 C(i, j) = (1.0 - invEj)*Fmatrix_()(i, j);
999 qBandI[i] += Fmatrix_()(i, j)*sigmaT4;
1004 Info<<
"Solving view factor equations for band :" 1013 if (iterCounter_ == 0)
1015 for (label i=0; i<totalNCoarseFaces_; i++)
1017 for (label j=0; j<totalNCoarseFaces_; j++)
1019 const scalar invEj = 1.0/E[j];
1022 CLU_()(i, j) = invEj;
1027 (1.0-invEj)*Fmatrix_()(i, j);
1037 for (label i=0; i<totalNCoarseFaces_; i++)
1039 for (label j=0; j<totalNCoarseFaces_; j++)
1041 const scalar sigmaT4 =
1046 qBandI[i] += -sigmaT4 + qrExt[j];
1050 qBandI[i] += Fmatrix_()(i, j)*sigmaT4;
1055 Info<<
"Solving view factor equations for band : " 1069 label globCoarseId = 0;
1070 for (
const label
patchID : selectedPatches_)
1072 const polyPatch&
pp = coarseMesh_->boundaryMesh()[
patchID];
1079 label nAgglom =
max(agglom)+1;
1084 coarseMesh_->patchFaceMap()[
patchID];
1087 forAll(coarseToFine, coarseI)
1089 label globalCoarse = globalNumbering.
toGlobal 1095 const label coarseFaceID = coarsePatchFace[coarseI];
1096 const labelList& fineFaces = coarseToFine[coarseFaceID];
1098 for (
const label facei : fineFaces)
1102 qrp[facei] = qNet[globalCoarse];
1106 qrp[facei] = qTotalCoarse[globCoarseId];
1117 for (
const label
patchID : selectedPatches_)
1121 const scalar heatFlux =
gSum(qrp*magSf);
1123 Info<<
"Total heat transfer rate at patch: " 1125 << heatFlux <<
endl;
label toLocal(const label proci, const label i) const
From global to local on proci.
Different types of constants.
List< scalar > scalarList
List of scalar.
virtual tmp< volScalarField > Rp() const
Source term component (for power of T^4)
const Type & value() const noexcept
Return const reference to value.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
fvMesh as subset of other mesh. Consists of one cell and all original boundary faces. Useful when manipulating boundary data. Single internal cell only needed to be able to manipulate in a standard way.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void append(const T &val)
Append an element at the end of the list.
List< edge > edgeList
List of edge.
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
static const boundaryRadiationProperties & New(const fvMesh &mesh, Args &&... args)
Get existing or create MeshObject registered with typeName.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
static rangeType allProcs(const label communicator=worldComm)
Range of process indices for all processes.
constexpr char nl
The newline '\n' character (0x0a)
Type & refCast(U &obj)
A dynamic_cast (for references) to Type reference.
label whichProcID(const label proci, const label i) const
Which processor does global id come from? Checks proci first (assumed to occur reasonably frequently)...
Ostream & endl(Ostream &os)
Add newline and flush stream.
static void gatherList(const UList< commsStruct > &comms, UList< T > &values, const int tag, const label comm)
Gather data, but keep individual values separate. Uses the specified communication schedule...
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
string upper(const std::string &s)
Return string copy transformed with std::toupper on each character.
static int & msgType() noexcept
Message tag of standard messages.
List< lduScheduleEntry > lduSchedule
A List of lduSchedule entries.
Ignore writing from objectRegistry::writeObject()
static int myProcNo(const label communicator=worldComm)
Rank of this process in the communicator (starting from masterNo()). Can be negative if the process i...
static label worldComm
Communicator for all ranks. May differ from commGlobal() if local worlds are in use.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
List< labelList > labelListList
List of labelList.
Macros for easy insertion into run-time selection tables.
static void broadcast(Type &value, const label comm=UPstream::worldComm)
Broadcast content (contiguous or non-contiguous) to all communicator ranks. Does nothing in non-paral...
#define forAll(list, i)
Loop across all elements in list.
labelList indices(const wordRe &matcher, const bool useGroups=true) const
Return (sorted) patch indices for all matches.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run. ...
const dimensionedScalar e
Elementary charge.
Type gSum(const FieldField< Field, Type > &f)
virtual bool read()=0
Read radiationProperties dictionary.
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
fvPatchField< scalar > fvPatchScalarField
const wordList area
Standard area field types (scalar, vector, tensor, etc)
List< scalarList > scalarListList
List of scalarList.
virtual tmp< volScalarField::Internal > Ru() const
Source term component (constant)
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const dimensionSet dimTemperature(0, 0, 0, 1, 0, 0, 0)
#define DebugInFunction
Report an information message using Foam::Info.
label size() const noexcept
The number of entries in the list.
static lduSchedule nonBlockingSchedule(const lduInterfacePtrsList &)
Get non-scheduled send/receive schedule.
UPtrList< const lduInterface > lduInterfacePtrsList
Store lists of lduInterface as a UPtrList.
void insertMatrixElements(const globalIndex &index, const label fromProci, const labelListList &globalFaceFaces, const scalarListList &viewFactors, scalarSquareMatrix &matrix)
Insert view factors into main matrix.
label localSize(const label proci) const
Size of proci data.
Top level model for radiation modelling.
volVectorField F(fluid.F())
void initialise()
Initialise.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
void calculate()
Solve system of equation(s)
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO...
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, IOobjectOption::registerOption regOpt, const Mesh &mesh, const dimensionSet &dims, const word &patchFieldType=fvPatchField< scalar >::calculatedType())
Return tmp field (NO_READ, NO_WRITE) from name, mesh, dimensions and patch type. [Takes current timeN...
static const word viewFactorWalls
Static name for view factor walls.
static autoPtr< solver > New(const word &solverName, const word &fieldName, const lduMatrix &matrix, const FieldField< Field, scalar > &interfaceBouCoeffs, const FieldField< Field, scalar > &interfaceIntCoeffs, const lduInterfaceFieldPtrsList &interfaces, const dictionary &solverControls)
Return a new solver of given type.
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
bool read()
Read radiation properties dictionary.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
Basic run-time type information using word as the type's name. Used to enhance the standard RTTI to c...
OSstream & masterStream(const label communicator)
Return OSstream for output operations on the master process only, Snull on other processes.
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
dimensionedScalar pow3(const dimensionedScalar &ds)
IOList< scalarList > scalarListIOList
IO for a List of scalarList.
label toGlobal(const label proci, const label i) const
From local to global on proci.
string lower(const std::string &s)
Return string copy transformed with std::tolower on each character.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
const polyBoundaryMesh & patches
static tmp< DimensionedField< Type, GeoMesh > > New(const word &name, IOobjectOption::registerOption regOpt, const Mesh &mesh, const dimensionSet &dims, const Field< Type > &iField)
Return tmp field (NO_READ, NO_WRITE) from name, mesh, dimensions, copy of internal field...
dimensionedScalar pow4(const dimensionedScalar &ds)
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
messageStream Info
Information stream (stdout output on master, null elsewhere)
label localStart(const label proci) const
Start of proci data.
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
List< label > labelList
A List of labels.
A class for managing temporary objects.
Defines the attributes of an object for which implicit objectRegistry management is supported...
#define addToRadiationRunTimeSelectionTables(model)
SquareMatrix< scalar > scalarSquareMatrix
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Do not request registration (bool: false)
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
autoPtr< radiation::radiationModel > radiation(radiation::radiationModel::New(T))
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
SolverPerformance< scalar > solverPerformance
SolverPerformance instantiated for a scalar.
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
Solve the matrix using LU decomposition with pivoting returning the LU form of the matrix and the sol...
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution in the source.
static void listCombineReduce(List< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Combines List elements. After completion all processors have the same data.
static constexpr const zero Zero
Global zero (0)
IOList< labelList > labelListIOList
IO for a List of labelList.