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,
97 for (
auto& pfld : bfld)
99 if (typeOnly ==
bool(isA<PatchType>(pfld)))
101 pfld.initEvaluate(commsType);
108 for (
auto& pfld : bfld)
110 if (typeOnly ==
bool(isA<PatchType>(pfld)))
112 pfld.evaluate(commsType);
121 const fvMatrix<Type>& m
132 const FieldField<Field, Type>& internalCoeffs = m.internalCoeffs();
135 forAll(internalCoeffs, patchi)
137 const labelUList& fc = mesh_.lduAddr().patchAddr(patchi);
138 const Field<Type>& intCoeffs = internalCoeffs[patchi];
139 const scalarField cmptCoeffs(intCoeffs.component(cmpt));
142 norm[fc[i]] += cmptCoeffs[i];
151 const scalar&
n = norm[celli];
163 Pout<<
"For field " << m.psi().
name() <<
" have zero diagonals for " 171 const labelList& own = mesh_.faceOwner();
172 const labelList& nei = mesh_.faceNeighbour();
182 extrapolatedNorm[celli] = -GREAT;
187 bitSet isFront(mesh_.nFaces());
188 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
190 label ownType = types[own[facei]];
191 label neiType = types[nei[facei]];
205 label facei = mesh_.nInternalFaces();
206 facei < mesh_.nFaces();
210 const label ownType = types[own[facei]];
211 const label neiType = nbrTypes[facei-mesh_.nInternalFaces()];
228 bitSet newIsFront(mesh_.nFaces());
232 for (
const label facei : isFront)
234 if (extrapolatedNorm[own[facei]] == -GREAT)
237 newNorm[own[facei]] = cellAverage
250 mesh_.isInternalFace(facei)
251 && extrapolatedNorm[nei[facei]] == -GREAT
255 newNorm[nei[facei]] = cellAverage
274 extrapolatedNorm.transfer(newNorm);
275 isFront.transfer(newIsFront);
282 scalar&
n = norm[celli];
292 n = extrapolatedNorm[celli];
305 const bool setHoleCellValue,
306 const Type& holeCellValue
319 Field<Type>& source = m.source();
325 const lduAddressing& addr = lduAddr();
326 const labelUList& upperAddr = addr.upperAddr();
327 const labelUList& lowerAddr = addr.lowerAddr();
330 if (!isA<fvMeshPrimitiveLduAddressing>(addr))
333 <<
"Problem : addressing is not fvMeshPrimitiveLduAddressing" 339 upper.setSize(upperAddr.size(), 0.0);
341 lower.setSize(lowerAddr.size(), 0.0);
346 const label nOldInterfaces = mesh_.fvMesh::interfaces().size();
349 if (interfaces.size() > nOldInterfaces)
352 m.internalCoeffs().setSize(interfaces.size());
353 m.boundaryCoeffs().setSize(interfaces.size());
358 label patchi = nOldInterfaces;
359 patchi < interfaces.size();
363 const labelUList& fc = interfaces[patchi].faceCells();
364 m.internalCoeffs().set(patchi,
new Field<Type>(fc.size(),
Zero));
365 m.boundaryCoeffs().set(patchi,
new Field<Type>(fc.size(),
Zero));
370 typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
372 typename GeoField::Boundary& bfld =
375 bfld.setSize(interfaces.size());
382 for (label patchi = 0; patchi < nOldInterfaces; patchi++)
384 if (isA<processorFvPatch>(bfld[patchi].
patch()))
393 label patchi = nOldInterfaces;
394 patchi < interfaces.size();
401 new calculatedProcessorFvPatchField<Type>
404 bfld[addPatchi].
patch(),
425 const label l = lowerAddr[facei];
426 scaleConnection(
upper, types, factor, setHoleCellValue, l, facei);
427 const label u = upperAddr[facei];
428 scaleConnection(
lower, types, factor, setHoleCellValue, u, facei);
445 for (label patchi = 0; patchi < nOldInterfaces; ++patchi)
447 const labelUList& fc = addr.patchAddr(patchi);
448 Field<Type>& bCoeffs = m.boundaryCoeffs()[patchi];
449 Field<Type>& iCoeffs = m.internalCoeffs()[patchi];
452 scaleConnection(bCoeffs, types, factor, setHoleCellValue, fc[i], i);
454 scaleConnection(iCoeffs, types, factor, setHoleCellValue, fc[i], i);
494 const Type wantedValue
500 diag[celli] = normalisation[celli];
501 source[celli] = normalisation[celli]*wantedValue;
509 const scalar
f = factor[celli];
512 const labelList& nbrFaces = stencilFaces_[celli];
513 const labelList& nbrPatches = stencilPatches_[celli];
515 diag[celli] *= (1.0-
f);
516 source[celli] *= (1.0-
f);
520 const label patchi = nbrPatches[nbri];
521 const label facei = nbrFaces[nbri];
525 const label nbrCelli = nbrs[nbri];
527 const scalar
s = normalisation[celli]*
f*w[nbri];
529 scalar& u =
upper[facei];
530 scalar& l =
lower[facei];
531 if (celli < nbrCelli)
559 const scalar
s = normalisation[celli]*
f*w[nbri];
560 m.boundaryCoeffs()[patchi][facei] += pTraits<Type>::one*
s;
561 m.internalCoeffs()[patchi][facei] += pTraits<Type>::one*
s;
720 typename GeoField::Boundary& bpsi =
723 bool hasOverset =
false;
737 Pout<<
"oversetFvMeshBase::solveOverset() :" 738 <<
" bypassing matrix adjustment for field " << m.
psi().name()
742 return mesh_.fvMesh::solve(m,
dict);
747 Pout<<
"oversetFvMeshBase::solveOverset() :" 748 <<
" adjusting matrix for interpolation for field " 755 if (
debug && mesh_.time().writeTime())
761 m.
psi().name() +
"_normalisation",
762 mesh_.time().timeName(),
775 oversetFvPatchField<scalar>
776 >(scale.boundaryFieldRef(),
false);
781 Pout<<
"oversetFvMeshBase::solveOverset() :" 782 <<
" writing matrix normalisation for field " << m.
psi().name()
783 <<
" to " << scale.name() <<
endl;
798 Field<Type> oldSource(m.
source());
801 const label oldSize = bpsi.
size();
808 correctBoundaryConditions<GeoField, calculatedProcessorFvPatchField<Type>>
816 SolverPerformance<Type>
s(mesh_.fvMesh::solve(m,
dict));
819 bpsi.setSize(oldSize);
825 m.
source().transfer(oldSource);
834 const labelList& own = mesh_.faceOwner();
835 const labelList& nei = mesh_.faceNeighbour();
838 const_cast<GeometricField<Type, fvPatchField, volMesh>&
>(m.
psi());
841 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
843 const label ownType = types[own[facei]];
844 const label neiType = types[nei[facei]];
851 psi[own[facei]] =
psi[nei[facei]];
859 psi[nei[facei]] =
psi[own[facei]];
879 os <<
"** Matrix **" <<
endl;
895 forAll(interfaces, patchi)
897 if (interfaces.
set(patchi))
899 const labelUList& fc = interfaces[patchi].faceCells();
903 cellToPatch[fc[i]].append(patchi);
904 cellToPatchFace[fc[i]].append(i);
912 os <<
"cell:" << celli <<
" diag:" <<
diag[celli]
913 <<
" source:" << source[celli] <<
endl;
915 label startLabel = ownerStartAddr[celli];
916 label endLabel = ownerStartAddr[celli + 1];
918 for (label facei = startLabel; facei < endLabel; facei++)
920 os <<
" face:" << facei
921 <<
" nbr:" << upperAddr[facei] <<
" coeff:" <<
upper[facei]
928 for (label i = startLabel; i < endLabel; i++)
930 label facei = losortAddr[i];
931 os <<
" face:" << facei
932 <<
" nbr:" << lowerAddr[facei] <<
" coeff:" <<
lower[facei]
936 forAll(cellToPatch[celli], i)
938 label patchi = cellToPatch[celli][i];
939 label patchFacei = cellToPatchFace[celli][i];
941 os <<
" patch:" << patchi
942 <<
" patchface:" << patchFacei
954 os <<
"patch:" << patchi
956 <<
" size:" << fc.size() <<
endl;
959 interfaces.
set(patchi)
960 && isA<processorLduInterface>(interfaces[patchi])
963 const processorLduInterface& ppp =
964 refCast<const processorLduInterface>(interfaces[patchi]);
965 os <<
"(processor with my:" << ppp.myProcNo()
966 <<
" nbr:" << ppp.neighbProcNo()
972 os <<
" cell:" << fc[i]
982 m.
psi().boundaryField().scalarInterfaces();
983 forAll(interfaceFields, inti)
985 if (interfaceFields.set(inti))
987 os <<
"interface:" << inti
988 <<
" if type:" << interfaceFields[inti].interface().type()
989 <<
" fld type:" << interfaceFields[inti].type() <<
endl;
993 os <<
"** End of Matrix **" <<
endl;
997 template<
class GeoField>
1000 auto& bfld =
fld.boundaryFieldRef();
1006 for (
auto& pfld : bfld)
1011 pfld.initEvaluate(commsType);
1018 for (
auto& pfld : bfld)
1023 pfld.evaluate(commsType);
1029 template<
class GeoField>
1032 Pout<<
"** starting checking of " <<
fld.name() <<
endl;
1034 GeoField fldCorr(
fld.name()+
"_correct",
fld);
1035 correctCoupledBoundaryConditions(fldCorr);
1037 const typename GeoField::Boundary& bfld =
fld.boundaryField();
1038 const typename GeoField::Boundary& bfldCorr = fldCorr.boundaryField();
1042 const typename GeoField::Patch& pfld = bfld[patchi];
1043 const typename GeoField::Patch& pfldCorr = bfldCorr[patchi];
1049 if (pfld[i] != pfldCorr[i])
1051 Pout<<
" " << i <<
" orig:" << pfld[i]
1052 <<
" corrected:" << pfldCorr[i] <<
endl;
UPtrList< const lduInterfaceField > lduInterfaceFieldPtrsList
List of coupled interface fields to be used in coupling.
List< scalar > scalarList
List of scalar.
virtual const scalarList & cellInterpolationWeight() const
Per interpolated cell the interpolation factor. (0 = use.
const scalarField & diag() const
void size(const label n)
Older name for setAddressableSize.
errorManipArg< error, int > exit(error &err, const int errNo=1)
commsTypes
Communications types.
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 cellCellStencilObject & New(const fvMesh &mesh, Args &&... args)
Get existing or create MeshObject registered with typeName.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static label nRequests() noexcept
Number of outstanding requests (on the internal list of 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.
cellMask correctBoundaryConditions()
Ostream & endl(Ostream &os)
Add newline and flush stream.
A traits class, which is primarily used for primitives and vector-space.
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.
::Foam::direction nComponents(const expressions::valueTypeCode) noexcept
The number of components associated with given valueTypeCode.
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.
virtual const fileName & name() const override
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
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.
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.
static void waitRequests()
Wait for all requests to finish.
UList< label > labelUList
A UList of labels.
virtual const labelUList & cellTypes() const
Return the cell type list.
#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.
label size() const noexcept
Return number of equations.
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
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.
virtual bool write(const bool writeOnProc=true) const
Write using setting from DB.
const scalarField & lower() const
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.
Us boundaryFieldRef().evaluateCoupled< coupledFaPatch >()
static void correctBoundaryConditions(typename GeoField::Boundary &bfld, const bool typeOnly)
Correct boundary conditions of certain type (typeOnly = true) or explicitly not of the type (typeOnly...
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)
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.
const scalarField & upper() const
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.
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)
Attempt dynamic_cast to Type.
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.
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.
Do not request registration (bool: false)
const labelUList & ownerStartAddr() const
Return owner start addressing.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Boundary condition for use on overset patches. To be run in combination with special dynamicFvMesh ty...
const T * set(const label i) const
Return const pointer to element (can be nullptr), or nullptr for out-of-range access (ie...
static constexpr const zero Zero
Global zero (0)