44 const bool setHoleCellValue,
49 const label cType = types[celli];
50 const scalar
f = factor[celli];
54 coeffs[facei] *= 1.0-
f;
71 coeffs[facei] *= 1.0-
f;
76 template<
class GeoField,
class PatchType>
79 typename GeoField::Boundary& bfld,
88 if (typeOnly == (isA<PatchType>(bfld[patchi]) !=
nullptr))
90 bfld[patchi].initEvaluate(commsType);
102 if (typeOnly == (isA<PatchType>(bfld[patchi]) !=
nullptr))
104 bfld[patchi].evaluate(commsType);
113 const fvMatrix<Type>& m
124 const FieldField<Field, Type>& internalCoeffs = m.internalCoeffs();
125 for (
direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
127 forAll(internalCoeffs, patchi)
129 const labelUList& fc = mesh_.lduAddr().patchAddr(patchi);
130 const Field<Type>& intCoeffs = internalCoeffs[patchi];
131 const scalarField cmptCoeffs(intCoeffs.component(cmpt));
134 norm[fc[i]] += cmptCoeffs[i];
143 const scalar&
n = norm[celli];
155 Pout<<
"For field " << m.psi().
name() <<
" have zero diagonals for " 163 const labelList& own = mesh_.faceOwner();
164 const labelList& nei = mesh_.faceNeighbour();
174 extrapolatedNorm[celli] = -GREAT;
179 bitSet isFront(mesh_.nFaces());
180 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
182 label ownType = types[own[facei]];
183 label neiType = types[nei[facei]];
197 label facei = mesh_.nInternalFaces();
198 facei < mesh_.nFaces();
202 const label ownType = types[own[facei]];
203 const label neiType = nbrTypes[facei-mesh_.nInternalFaces()];
220 bitSet newIsFront(mesh_.nFaces());
224 for (
const label facei : isFront)
226 if (extrapolatedNorm[own[facei]] == -GREAT)
229 newNorm[own[facei]] = cellAverage
242 mesh_.isInternalFace(facei)
243 && extrapolatedNorm[nei[facei]] == -GREAT
247 newNorm[nei[facei]] = cellAverage
266 extrapolatedNorm.transfer(newNorm);
267 isFront.transfer(newIsFront);
274 scalar&
n = norm[celli];
284 n = extrapolatedNorm[celli];
297 const bool setHoleCellValue,
298 const Type& holeCellValue
311 Field<Type>& source = m.source();
317 const lduAddressing& addr = lduAddr();
318 const labelUList& upperAddr = addr.upperAddr();
319 const labelUList& lowerAddr = addr.lowerAddr();
322 if (!isA<fvMeshPrimitiveLduAddressing>(addr))
325 <<
"Problem : addressing is not fvMeshPrimitiveLduAddressing" 331 upper.setSize(upperAddr.size(), 0.0);
333 lower.setSize(lowerAddr.size(), 0.0);
338 const label nOldInterfaces = mesh_.fvMesh::interfaces().size();
341 if (interfaces.size() > nOldInterfaces)
344 m.internalCoeffs().setSize(interfaces.size());
345 m.boundaryCoeffs().setSize(interfaces.size());
350 label patchi = nOldInterfaces;
351 patchi < interfaces.size();
355 const labelUList& fc = interfaces[patchi].faceCells();
356 m.internalCoeffs().set(patchi,
new Field<Type>(fc.size(),
Zero));
357 m.boundaryCoeffs().set(patchi,
new Field<Type>(fc.size(),
Zero));
362 typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
364 typename GeoField::Boundary& bfld =
365 const_cast<GeoField&
>(m.psi()).boundaryFieldRef();
367 bfld.setSize(interfaces.size());
374 for (label patchi = 0; patchi < nOldInterfaces; patchi++)
376 if (isA<processorFvPatch>(bfld[patchi].
patch()))
385 label patchi = nOldInterfaces;
386 patchi < interfaces.size();
393 new calculatedProcessorFvPatchField<Type>
396 bfld[addPatchi].
patch(),
417 const label l = lowerAddr[facei];
418 scaleConnection(
upper, types, factor, setHoleCellValue, l, facei);
419 const label u = upperAddr[facei];
420 scaleConnection(
lower, types, factor, setHoleCellValue, u, facei);
437 for (label patchi = 0; patchi < nOldInterfaces; ++patchi)
439 const labelUList& fc = addr.patchAddr(patchi);
440 Field<Type>& bCoeffs = m.boundaryCoeffs()[patchi];
441 Field<Type>& iCoeffs = m.internalCoeffs()[patchi];
444 scaleConnection(bCoeffs, types, factor, setHoleCellValue, fc[i], i);
446 scaleConnection(iCoeffs, types, factor, setHoleCellValue, fc[i], i);
486 const Type wantedValue
492 diag[celli] = normalisation[celli];
493 source[celli] = normalisation[celli]*wantedValue;
501 const scalar
f = factor[celli];
504 const labelList& nbrFaces = stencilFaces_[celli];
505 const labelList& nbrPatches = stencilPatches_[celli];
507 diag[celli] *= (1.0-
f);
508 source[celli] *= (1.0-
f);
512 const label patchi = nbrPatches[nbri];
513 const label facei = nbrFaces[nbri];
517 const label nbrCelli = nbrs[nbri];
519 const scalar
s = normalisation[celli]*
f*w[nbri];
521 scalar& u =
upper[facei];
522 scalar& l =
lower[facei];
523 if (celli < nbrCelli)
551 const scalar
s = normalisation[celli]*
f*w[nbri];
552 m.boundaryCoeffs()[patchi][facei] += pTraits<Type>::one*
s;
553 m.internalCoeffs()[patchi][facei] += pTraits<Type>::one*
s;
712 typename GeoField::Boundary& bpsi =
713 const_cast<GeoField&
>(m.
psi()).boundaryFieldRef();
715 bool hasOverset =
false;
729 Pout<<
"oversetFvMeshBase::solveOverset() :" 730 <<
" bypassing matrix adjustment for field " << m.
psi().name()
734 return mesh_.fvMesh::solve(m,
dict);
739 Pout<<
"oversetFvMeshBase::solveOverset() :" 740 <<
" adjusting matrix for interpolation for field " 747 if (
debug && mesh_.time().outputTime())
753 m.
psi().name() +
"_normalisation",
754 mesh_.time().timeName(),
767 oversetFvPatchField<scalar>
768 >(scale.boundaryFieldRef(),
false);
773 Pout<<
"oversetFvMeshBase::solveOverset() :" 774 <<
" writing matrix normalisation for field " << m.
psi().name()
775 <<
" to " << scale.name() <<
endl;
790 Field<Type> oldSource(m.
source());
793 const label oldSize = bpsi.
size();
800 correctBoundaryConditions<GeoField, calculatedProcessorFvPatchField<Type>>
808 SolverPerformance<Type>
s(mesh_.fvMesh::solve(m,
dict));
811 bpsi.setSize(oldSize);
817 m.
source().transfer(oldSource);
826 const labelList& own = mesh_.faceOwner();
827 const labelList& nei = mesh_.faceNeighbour();
830 const_cast<GeometricField<Type, fvPatchField, volMesh>&
>(m.
psi());
833 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
835 const label ownType = types[own[facei]];
836 const label neiType = types[nei[facei]];
843 psi[own[facei]] =
psi[nei[facei]];
851 psi[nei[facei]] =
psi[own[facei]];
871 os <<
"** Matrix **" <<
endl;
887 forAll(interfaces, patchi)
889 if (interfaces.
set(patchi))
891 const labelUList& fc = interfaces[patchi].faceCells();
895 cellToPatch[fc[i]].append(patchi);
896 cellToPatchFace[fc[i]].append(i);
904 os <<
"cell:" << celli <<
" diag:" <<
diag[celli]
905 <<
" source:" << source[celli] <<
endl;
907 label startLabel = ownerStartAddr[celli];
908 label endLabel = ownerStartAddr[celli + 1];
910 for (label facei = startLabel; facei < endLabel; facei++)
912 os <<
" face:" << facei
913 <<
" nbr:" << upperAddr[facei] <<
" coeff:" <<
upper[facei]
920 for (label i = startLabel; i < endLabel; i++)
922 label facei = losortAddr[i];
923 os <<
" face:" << facei
924 <<
" nbr:" << lowerAddr[facei] <<
" coeff:" <<
lower[facei]
928 forAll(cellToPatch[celli], i)
930 label patchi = cellToPatch[celli][i];
931 label patchFacei = cellToPatchFace[celli][i];
933 os <<
" patch:" << patchi
934 <<
" patchface:" << patchFacei
946 os <<
"patch:" << patchi
948 <<
" size:" << fc.size() <<
endl;
951 interfaces.
set(patchi)
952 && isA<processorLduInterface>(interfaces[patchi])
955 const processorLduInterface& ppp =
956 refCast<const processorLduInterface>(interfaces[patchi]);
957 os <<
"(processor with my:" << ppp.myProcNo()
958 <<
" nbr:" << ppp.neighbProcNo()
964 os <<
" cell:" << fc[i]
974 m.
psi().boundaryField().scalarInterfaces();
975 forAll(interfaceFields, inti)
977 if (interfaceFields.set(inti))
979 os <<
"interface:" << inti
980 <<
" if type:" << interfaceFields[inti].interface().type()
981 <<
" fld type:" << interfaceFields[inti].type() <<
endl;
985 os <<
"** End of Matrix **" <<
endl;
989 template<
class GeoField>
992 typename GeoField::Boundary& bfld =
fld.boundaryFieldRef();
1002 bfld[patchi].initEvaluate(commsType);
1017 bfld[patchi].evaluate(commsType);
1023 template<
class GeoField>
1026 Pout<<
"** starting checking of " <<
fld.name() <<
endl;
1028 GeoField fldCorr(
fld.name()+
"_correct",
fld);
1029 correctCoupledBoundaryConditions(fldCorr);
1031 const typename GeoField::Boundary& bfld =
fld.boundaryField();
1032 const typename GeoField::Boundary& bfldCorr = fldCorr.boundaryField();
1036 const typename GeoField::Patch& pfld = bfld[patchi];
1037 const typename GeoField::Patch& pfldCorr = bfldCorr[patchi];
1043 if (pfld[i] != pfldCorr[i])
1045 Pout<<
" " << i <<
" orig:" << pfld[i]
1046 <<
" corrected:" << pfldCorr[i] <<
endl;
UPtrList< const lduInterfaceField > lduInterfaceFieldPtrsList
List of coupled interface fields to be used in coupling.
virtual const scalarList & cellInterpolationWeight() const
Per interpolated cell the interpolation factor. (0 = use.
virtual const fileName & name() const
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
List< labelList > labelListList
A List of labelList.
void size(const label n)
Older name for setAddressableSize.
errorManipArg< error, int > exit(error &err, const int errNo=1)
commsTypes
Types of communications.
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
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.
const labelUList & losortStartAddr() const
Return losort start addressing.
static const Type & New(const Mesh &mesh, Args &&... args)
Get existing or create a new MeshObject.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static label nRequests() noexcept
Number of outstanding requests.
static void checkCoupledBC(const GeoField &fld)
Debug: check halo swap is ok.
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
cellMask correctBoundaryConditions()
Ostream & endl(Ostream &os)
Add newline and flush stream.
A traits class, which is primarily used for primitives.
string upper(const std::string &s)
Return string copy transformed with std::toupper on each character.
void scaleConnection(Field< Type > &coeffs, const labelUList &types, const scalarList &factor, const bool setHoleCellValue, const label celli, const label facei) const
Freeze values at holes.
Generic GeometricField class.
tmp< scalarField > normalisation(const fvMatrix< Type > &m) const
Determine normalisation for interpolation. This equals the original diagonal except stabilised for ze...
Ignore writing from objectRegistry::writeObject()
const dimensionSet dimless
Dimensionless.
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.
virtual const labelUList & lowerAddr() const =0
Return lower addressing.
static void correctCoupledBoundaryConditions(GeoField &fld)
Debug: correct coupled bc.
SolverPerformance< Type > solveOverset(fvMatrix< Type > &m, const dictionary &) const
Solve given dictionary with settings.
UList< label > labelUList
A UList of labels.
virtual const labelUList & cellTypes() const
Return the cell type list.
const T * set(const label i) const
Return const pointer to element (can be nullptr), or nullptr for out-of-range access (ie...
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
const GeometricField< Type, fvPatchField, volMesh > & psi(const label i=0) const
Return psi.
Generic templated field type.
virtual const labelListList & cellStencil() const
Per interpolated cell the neighbour cells (in terms of slots as.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual const labelUList & upperAddr() const =0
Return upper addressing.
void addInterpolation(fvMatrix< Type > &m, const scalarField &normalisation, const bool setHoleCellValue, const Type &holeCellValue) const
Manipulate the matrix to add the interpolation/set hole.
const cellCellStencilObject & overlap
List< scalar > scalarList
A List of scalars.
virtual const List< scalarList > & cellInterpolationWeights() const
Weights for cellStencil.
UPtrList< const lduInterface > lduInterfacePtrsList
Store lists of lduInterface as a UPtrList.
virtual const labelUList & patchAddr(const label patchNo) const =0
Return patch to internal addressing given patch number.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
static void correctBoundaryConditions(typename GeoField::Boundary &bfld, const bool typeOnly)
Correct boundary conditions of certain type (typeOnly = true)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
int debug
Static debugging option.
OBJstream os(runTime.globalPath()/outputName)
static void waitRequests(const label start=0)
Wait until all requests (from start onwards) have finished.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< ' ';}gmvFile<< nl;for(const word &name :lagrangianScalarNames){ IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const Field< Type > & field() const noexcept
Return const-reference to the field values.
Field< Type > & source() noexcept
static commsTypes defaultCommsType
Default commsType.
string lower(const std::string &s)
Return string copy transformed with std::tolower on each character.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const std::string patch
OpenFOAM patch number as a std::string.
const labelUList & losortAddr() const
Return losort addressing.
"nonBlocking" : (MPI_Isend, MPI_Irecv)
Internal & ref(const bool updateAccessTime=true)
Same as internalFieldRef()
const FieldField< Field, Type > & internalCoeffs() const noexcept
fvBoundary scalar field containing pseudo-matrix coeffs for internal cells
const volScalarField & psi
The class contains the addressing required by the lduMatrix: upper, lower and losort.
const FieldField< Field, Type > & boundaryCoeffs() const noexcept
fvBoundary scalar field containing pseudo-matrix coeffs for boundary cells
void boundaryManipulate(typename GeometricField< Type, fvPatchField, volMesh >::Boundary &values)
Manipulate based on a boundary field.
const Type * isA(const U &obj)
Check if dynamic_cast to Type is possible.
List< label > labelList
A List of labels.
A class for managing temporary objects.
void write(Ostream &, const fvMatrix< Type > &, const lduAddressing &, const lduInterfacePtrsList &) const
Debug: print matrix.
virtual bool write(const bool valid=true) const
Write using setting from DB.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
const labelUList & ownerStartAddr() const
Return owner start addressing.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
label size() const
Return number of equations.
Boundary condition for use on overset patches. To be run in combination with special dynamicFvMesh ty...
static constexpr const zero Zero
Global zero (0)